ContactList.h

Go to the documentation of this file.

 
#include <vector>
#include <array>
#include <list>
#include "Typedefs.h"
#include "Parameters.h"
#include "Tools.h"
 
 
#ifndef CONTACTLIST
#define CONTACTLIST
 
// ------------------------------------ Action class -------------------------------------------
template <int d>
class Action {
public :
    Action() : Fn(d,0), Ft(d,0), Torquei(d*(d-1)/2,0), Torquej(d*(d-1)/2,0), vn(d,0), vt(d,0) {} 
    vector <double> Fn, Ft, Torquei, Torquej ;
    vector <double> vn, vt ;
    vector<double> Fn_el, Fn_visc, Ft_el, Ft_visc, Ft_fric ; 
    bool Ft_isfric ;  
    
    template <class Archive>
    void serialize(Archive &ar) {
        ar(Fn, Ft, Torquei, Torquej, vn, vt, Fn_el, Fn_visc, Ft_el, Ft_visc, Ft_fric, Ft_isfric) ; 
    }
    
    void set (cv1d &a, cv1d &b, cv1d &c, cv1d &dd) {Fn=a ; Ft=b ; Torquei=c ; Torquej=dd ; }
    void setvel(cv1d &vvn, cv1d &vvt) {vn=vvn; vt=vvt ;}
    void setzero () 
    {
        for (int dd=0 ; dd<d ; dd++) Fn[dd]=Ft[dd]=vn[dd]=vt[dd]=0 ;    
        for (int dd=0 ; dd<d*(d-1)/2 ; dd++) Torquei[dd]=Torquej[dd]=0 ;    
    }
    
    //void set_fwd (v1d a, v1d b, v1d c) {Fni=a; Fti=b; Torquei=c ;}
    //  void set_rev (v1d a, v1d b, v1d c) {Fnj=a; Ftj=b; Torquej=c ;}
} ;
template <int d>
class SpecificAction: public Action<d> {
public :
    int id ;
    int duration ;
    template <class Archive>
    void serialize(Archive &ar) {
        ar(cereal::base_class<Action<d>>(this), id, duration) ; 
    }
} ;
 
// ------------------------------------ Contact properties class -------------------------------------------
template <int d>
class cp
{
public:
 cp () {} 
 cp (int ii, int jj, double ctlength, Action<d> * default_action) : i(ii), j(jj), contactlength(ctlength), tspr (vector <double> (d, 0)), infos(default_action), owninfos(false), persisting(true) {} 
 cp(const cp& v) { *this=v ; }
 ~cp () { if (owninfos) delete(infos) ; } 
 cp & operator= (const cp & c)
 {
     i=c.i ;
     if (c.i<0) return *this ; //is a deleted element, just keep moving
     j=c.j ; // copy everything else if it is not a deleted element
     ghost=c.ghost ;
     ghostdir=c.ghostdir ;
     tspr=c.tspr ;
     contactlength=c.contactlength ;
     owninfos = c.owninfos ; 
     infos=c.infos ;
     persisting = c.persisting ; 
     return *this ;
 } 
 
 Action<d> & getinfo () {return *infos ; } 
 //void setinfo (Action & a) {if (!infos) infos = new Action ; *infos=a ; }
 void setinfo (Action<d> * a) {infos=a ; } 
 void saveinfo (Action<d> & a) {
     if (!owninfos)
     {
         infos = new Action<d> ;
         owninfos = true ;
     }
     infos->Fn = a.Fn ;
     //printf("%g %g %g\n",a.Fn[0], a.Fn[1], a.Fn[2]) ;
     infos->Ft = a.Ft ;
     infos->vn = a.vn ;
     infos->vt = a.vt ;
     infos->Fn_el = a.Fn_el ;
     infos->Fn_visc = a.Fn_visc ;
     infos->Ft_el = a.Ft_el ;
     infos->Ft_visc = a.Ft_visc ;
     infos->Ft_fric = a.Ft_fric ; 
     infos->Torquei = a.Torquei ;
     infos->Torquej = a.Torquej ;
 } 
 
 int i ;  
 int j ; 
 double contactlength ; 
 uint32_t ghost ; 
 uint32_t ghostdir ; 
 std::vector<double> tspr ; 
 Action<d> * infos ; 
 bool owninfos ; 
 bool persisting ; 
 
 template <class Archive>
 void save (Archive &ar) const {
        ar(i, j, contactlength, ghost, ghostdir, tspr, owninfos, persisting) ;
        if (owninfos)
            ar(*infos) ; 
    }
 template <class Archive>
 void load (Archive &ar) {
        ar(i, j, contactlength, ghost, ghostdir, tspr, owninfos, persisting) ;
        if (owninfos)
        {
            infos = new Action<d> ; 
            ar(*infos) ; 
        }
    }
 
 std::pair<vector<double>,vector<double>> compute_branchvector (cv2d &X, Parameters<d> &P)
 {
    vector <double> loc (d, 0), branch (d, 0)  ;
    if ( P.Boundaries[0].Type != WallType::PBC_LE || (ghost & 1)==0)
    {
        loc=X[j] ;
        uint32_t gh=ghost, ghd=ghostdir ;
        for (int n=0 ; gh>0 ; gh>>=1, ghd>>=1, n++)
            if (gh&1)
                loc[n] += P.Boundaries[n].delta * ((ghd&1)?-1:1) ;
    }
    else
    {
        uint32_t gh=ghost, ghd=ghostdir ; // Handle pbc in first dim
        loc=X[j] ;
        loc[0] += P.Boundaries[0].delta * ((ghd&1)?-1:1) ;
        loc[1] += (ghd&1?-1:1)*P.Boundaries[0].displacement ;
        if (P.contact_strategy != ContactStrategies::CELLS)
        {
            double additionaldelta = 0 ;
            if (loc[1] > P.Boundaries[1].xmax) {additionaldelta = -P.Boundaries[1].delta ;}
            if (loc[1] < P.Boundaries[1].xmin) {additionaldelta =  P.Boundaries[1].delta ;}
            loc[1] += additionaldelta ;
        }
 
        gh>>=1 ; ghd>>=1 ;
        for (int n=1 ; gh>0 ; gh>>=1, ghd>>=1, n++)
            if (gh&1)
                loc[n] += P.Boundaries[n].delta * ((ghd&1)?-1:1) ;
    }
    for (int dd = 0 ; dd<d ; dd++) branch[dd] = X[i][dd]-loc[dd] ; //j->i
    for (int dd = 0 ; dd<d ; dd++) loc[dd] = (X[i][dd]+loc[dd]*(P.r[i]/(P.r[i]+P.r[j]) )) ; // This location is wrong for sphere with different radius.
    return {loc, branch} ; 
 } 
} ;
 
template <int d>
class cpm : public cp<d> {
public: 
    cpm () {} 
    cpm (int ii, int jj, int sid, double ctlength, Action<d> * default_action) : cp<d>(ii,jj,ctlength,default_action), contactpoint(std::vector<double>(d,0)), submeshid(sid){} 
    cpm(const cpm& v): cp<d>(v) {*this=v ;}
    ~cpm () {} 
    cpm & operator= (const cpm & c)
    {
     cp<d>::operator=(c);
     contactpoint = c.contactpoint ;
     submeshid = c.submeshid ; 
     return *this ;
    } 
    vector <double> contactpoint; 
    int submeshid ; 
    
    template <class Archive>
    void load (Archive &ar) {
        ar( cereal::base_class<cp<d>>( this ), contactpoint, submeshid );
    }
    template <class Archive>
    void save (Archive &ar) const {
        ar( cereal::base_class<cp<d>>( this ), contactpoint, submeshid );
    }
} ; 
 
// ------------------------------------ Contact List class -------------------------------------------
template <int d>
class ContactList
{
public:
 ContactList () {check_ghost=&ContactList::check_ghost_LE;}
 void reset() {it = v.begin() ;}  
 int insert(const cp<d>& a) ; 
 void finalise () { while (it!=v.end()) it=v.erase(it) ; } 
 
 // Functions for Cell contact detections
 /*int init_for_cells(int N)
 {
     it_array_beg.resize(N,null_list.begin()) ;
     it_array_end.resize(N,null_list.begin()) ;
     return 0 ; 
 }
 int make_iterator_array(int N) 
 {
     for (int i=0 ; i<N ; i++) it_array_beg[i] = null_list.begin() ;
     while (v.begin()!=v.end() && !(v.begin()->persisting)) v.erase(v.begin()) ;
     
     int curi = -1 ; int n=0 ; 
     auto tmpit = v.begin() ; 
     for (auto it = v.begin() ; it!=v.end() ; it++, n++)
     {
         if (it->i != curi)
         {
             it_array_beg[it->i] = it ; 
             if (curi!= -1) it_array_end[curi] = it ; 
             curi = it->i ; 
         }
         tmpit = it ; 
         while (++tmpit != v.end() && !tmpit->persisting) {v.erase(tmpit) ; tmpit = it ;  }
     }
     if (curi!= -1) it_array_end[curi] = v.end() ; 
     return 0 ; 
 }
 std::pair<typename list<cp<d>>::iterator, typename list<cp<d>>::iterator> it_bounds (int nmin, int nmax, int N)
 {   
     for ( ; nmin < N && it_array_beg[nmin] == null_list.begin() ; nmin++) ;
     for (--nmax ; nmax >=0 && it_array_beg[nmax] == null_list.begin() ; nmax--) ;
     
     if (nmax<nmin) return {v.end(), v.end()} ; 
     return {it_array_beg[nmin], it_array_end[nmax]} ; 
 }
 void coalesce_list ()
 {
     std::vector<cp<d>> vtmp ; 
     vtmp.reserve(v.size()) ; 
     printf("%ld ", v.size()) ; fflush(stdout) ; 
     int i = 0 ; 
     for (auto& element : v) { printf("\r%d", i++) ; fflush(stdout) ; vtmp.push_back(std::move(element)); }
     v.clear();
     printf("\n") ; fflush(stdout) ;
     i=0 ; 
     for (auto& element : vtmp) {printf("\r%d", i++) ; fflush(stdout) ; v.push_back(element); }
     return ;      
 }*/
 
 //void check_ghost    (uint32_t gst, double partialsum, const Parameters & P, cv1d &X1, cv1d &X2, double R, cp & tmpcp) ;
 void check_ghost_dst(uint32_t gst, int n, double partialsum, uint32_t mask, const Parameters<d> & P, cv1d &X1, cv1d &X2, cp<d> & contact) ; 
 bool check_ghost_regular (bitdim gst, const Parameters<d> & P, cv1d &X1, cv1d &X2, double r1, double r2, cp<d> & tmpcp,
                   int startd=0, double partialsum=0, bitdim mask=0) ; 
 bool check_ghost_LE (bitdim gst, const Parameters<d> & P, cv1d &X1, cv1d &X2, double r1, double r2, cp<d> & tmpcp, int startd=0, double partialsum=0, bitdim mask=0) ;
 bool (ContactList::*check_ghost) (bitdim , const Parameters<d> & , cv1d &, cv1d &, double, double, cp<d> &,int startd, double partialsum, bitdim mask) ;
 void coordinance (v1d &Z) ; 
 Action<d> * default_action () {return (&def) ; } 
  
 template <class Archive>
 void serialize(Archive &ar) {ar(v, cid) ; }
 
 list <cp<d>> v ; 
 int cid=0 ; 
 //std::vector<typename list<cp<d>>::iterator> it_array_beg, it_array_end ;
 
private:
 typename list<cp<d>>::iterator it ; 
 Action<d> def ; 
 
 list<cp<d>> null_list ; 
};
 
template <int d>
inline bool operator< (const cp<d> &a, const cp<d> &b) {if (a.i==b.i) return (a.j<b.j) ; return a.i<b.i ; } 
template <int d>
inline bool operator== (const cp<d> &a, const cp<d> &b) {return (a.i==b.i && a.j==b.j) ; } 
 
//-------------------------------------------------
template <int d>
class ContactListMesh
{
public:
 list <cpm<d>> v ;    
 Action<d> * default_action () {return (&def) ; } 
 void reset() {it = v.begin() ;} 
 void finalise () { while (it!=v.end()) it=v.erase(it) ; }
 int insert(const cpm<d>& a) ; 
 int cid=0 ; 
 
 template <class Archive>
 void serialize(Archive &ar) {ar(v, cid) ; }
 
 bool check_mesh_dst_contact (Mesh<d> &mesh, cv1d & Xo, double r, cpm<d> & c)
 {     
  std::vector<double> dotproducts (d) ;
  auto X = Xo-mesh.origin ;
  
  // Treating the special case of having a single point
  if (mesh.dimensionality == 0)
  {
      double dst=Tools<d>::norm(X) ; 
      if (dst<r) // We got a contact
      {
          c.contactlength = dst ; 
          c.contactpoint=mesh.origin ; 
          c.submeshid=mesh.dimensionality ; 
          return true ; 
      }
      return false ; // There can't be any submesh for a 0-dimensionality mesh
  }
  
  // All other cases. 
  double dstsqr=0 ;
  std::vector<double> coefficient (d, 0) ; 
  for (int i=0 ; i<d-mesh.dimensionality ; i++)
  {
      for (int dd=0 ; dd<d ; dd++)
          coefficient[i] += mesh.invertbase[i*d+dd]*X[dd] ; 
      dstsqr += coefficient[i]*coefficient[i] ; 
      if (dstsqr > r*r) goto submeshesprocessing ;       
  }
  if (dstsqr<r*r)
  {
    double coefficient_sum = 0 ; 
    for (int i=d-mesh.dimensionality ; i<d ; i++)
    {
      for (int dd=0 ; dd<d ; dd++)
          coefficient[i] += mesh.invertbase[i*d+dd]*X[dd] ;
      if (coefficient[i]<0 || coefficient[i]>1) goto submeshesprocessing ; 
      coefficient_sum+=coefficient[i] ; 
      if (coefficient_sum > 1) goto submeshesprocessing ; 
    }
    if (coefficient_sum < 1)
    {
      c.contactlength = sqrt(dstsqr) ;
      c.contactpoint=mesh.origin ; 
      c.submeshid=mesh.dimensionality ; 
      for (int i=d-mesh.dimensionality ; i<d ; i++)
        for (int dd=0 ; dd<d ; dd++)
          c.contactpoint[dd] += mesh.mixedbase[i*d+dd]*coefficient[i] ;
      //for (int i=0 ; i<d-mesh.dimensionality ; i++) printf("%g ", coefficient[i]) ;
      //printf("#") ; 
      //for (int i=d-mesh.dimensionality ; i<d; i++) printf("%g ", coefficient[i]) ; 
      
      //Tools<d>::print(mesh.mixedbase) ; 
      
      return true ;
    }
  }
  
  
  
  
  
  /*for (int i=0; i<d-mesh.dimensionality ; i++)
  {
    dotproducts[i]=Tools<d>::dot(mesh.mixedbase[i], X) ; 
    dstsqr+= dotproducts[i]*dotproducts[i] ; 
    //if (mesh.dimensionality==2)  printf("%g %g %g %g|", mesh.mixedbase[i][0], mesh.mixedbase[i][1], mesh.mixedbase[i][2], dotproducts[i]) ; 
  }
  if (dstsqr<r*r) //Potential contact
  {
    std::vector<double> coefficient(d) ;
    double coefficient_sum = 0 ;
    //printf("-> %g | [%g %g %g | %g %g %g | %g %g %g] ", sqrt(dstsqr), mesh.mixedbase[0][0], mesh.mixedbase[0][1], mesh.mixedbase[0][2], mesh.mixedbase[1][0], mesh.mixedbase[1][1], mesh.mixedbase[1][2],mesh.mixedbase[2][0], mesh.mixedbase[2][1], mesh.mixedbase[2][2]) ; 
    for (int i=d-mesh.dimensionality ; i<d ; i++)
    {
      dotproducts[i] = Tools<d>::dot(mesh.mixedbase[i], X) ; 
      coefficient[i] = dotproducts[i]/mesh.norms[i]/mesh.norms[i] ; 
      //printf("/%d %g ", i, coefficient[i]) ; fflush(stdout) ; 
      if (coefficient[i]<0 || coefficient[i]>1) goto submeshesprocessing ;  // We are outside        
      coefficient_sum += coefficient[i] ; 
      if (coefficient_sum>1) goto submeshesprocessing ;  
    }
    if (coefficient_sum>=0 && coefficient_sum<=1) // Seems like we got a contact!
    {
      c.contactlength = sqrt(dstsqr) ;
      c.contactpoint=mesh.origin ; 
      for (int i=d-mesh.dimensionality ; i<d ; i++)
        c.contactpoint += mesh.mixedbase[i]*coefficient[i] ; 
      return true ; 
    }
  }*/
  
  submeshesprocessing:                  // yeah yeah, using goto ... sue me
  if (mesh.submeshes.size()>0)
      for (int i=mesh.dimensionality-1 ; i>=0 ; i--)
          for (size_t j=0 ; j<mesh.submeshes[i].size() ; j++)
          {
              bool res = check_mesh_dst_contact(mesh.submeshes[i][j], Xo, r, c) ; 
              if (res) return true ; 
          }
  return false ; 
 }
  
private: 
 typename list<cpm<d>>::iterator it ;  
 Action<d> def ; 
} ; 
 
 
/*****************************************************************************************************
 *                                                                                                   *
 *                                                                                                   *
 *                                                                                                   *
 * IMPLEMENTATIONS                                                                                   *
 *                                                                                                   *
 *                                                                                                   *
 *                                                                                                   *
 * ***************************************************************************************************/
 
template <int d>
int ContactList<d>::insert(const cp<d> &a)
{
    while (it!=v.end() && (*it) < a)
    { it= v.erase(it) ; }
    if (it==v.end()) {v.push_back(a) ; it=v.end() ; }
    else
    {
        if ((*it)==a)
        {
            it->contactlength=a.contactlength ;
            it->ghost=a.ghost ;
            it->ghostdir=a.ghostdir ;
            it->persisting = true ; 
            it++ ;
        }
        else {it=v.insert(it,a) ; it++ ; }
    }
    return (cid++) ;
}
//----------------------------------------------
template <int d>
int ContactListMesh<d>::insert(const cpm<d> &a)
{
    while (it!=v.end() && (*it) < a)
    { it= v.erase(it) ; }
    if (it==v.end()) {v.push_back(a) ; it=v.end() ; }
    else
    {
        if ((*it)==a)
        {
            it->contactlength=a.contactlength ;
            it->ghost=a.ghost ;
            it->ghostdir=a.ghostdir ;
            it->contactpoint = a.contactpoint ; 
            it->submeshid = a.submeshid ; 
            it++ ;
        }
        else {it=v.insert(it,a) ; it++ ; }
    }
    return (cid++) ;
}
 
//-----------------------------------Fastest version so far ...
template <int d>
bool ContactList<d>::check_ghost_regular (bitdim gst, const Parameters<d> & P, 
                                          cv1d &X1, cv1d &X2, double r1, double r2, 
                                          cp<d> & tmpcp, int startd, double partialsum, bitdim mask)
{
    double sum=partialsum ;
    double rr = (r1+r2)*(r1+r2) ; 
    for (int dd=startd ; sum<rr && dd<d ; dd++, gst>>=1)
    {
        sum += (X1[dd]-X2[dd]) * (X1[dd]-X2[dd]) ;
        if (gst & 1)
        {
            double Delta= (tmpcp.ghostdir&(1<<dd)?-1:1) * P.Boundaries[dd].delta ;
            double sumspawn = partialsum + (X1[dd]-X2[dd]-Delta) * (X1[dd]-X2[dd]-Delta) ;
            //printf("/%g %g %g %g %g %g %g/", partialsum, sumspawn, X1[0], X1[1], X2[0], X2[1], Delta ) ;
            if (sumspawn<rr)
            {
                bool res = check_ghost_regular (gst>>1, P, X1, X2, r1, r2, tmpcp, dd+1, sumspawn, mask | (1<<dd)) ;
                if (res) return true ; 
            }
        }
        partialsum = sum ;
    }
    if (sum<rr)
    {
        tmpcp.contactlength=sqrt(sum) ;
        tmpcp.ghost=mask ;
        insert(tmpcp) ;
        return true ; 
        //printf("[%d %d %X %X %g]\n", tmpcp.i, tmpcp.j, tmpcp.ghost, tmpcp.ghostdir, tmpcp.contactlength) ;
    }
    return false ; 
}
 
//--------------------------------------------------------------------
template <int d>
bool ContactList<d>::check_ghost_LE (bitdim gst, const Parameters<d> & P, cv1d &X1, cv1d &X2, double r1, double r2, cp<d> & tmpcp,
                                     [[maybe_unused]] int startd, [[maybe_unused]]double partialsum, [[maybe_unused]]bitdim mask)
{
    if (P.Boundaries[0].Type != WallType::PBC_LE)
            return check_ghost_regular(gst, P, X1, X2, r1, r2, tmpcp) ;
    else
    {
     if ((gst & 1)==0)
     {
         double partialsum = (X1[0]-X2[0]) * (X1[0]-X2[0]) ;
         assert (((gst>>30 & 1) ==0)) ;
         check_ghost_regular(gst>>1, P, X1, X2, r1, r2, tmpcp, 1, partialsum, 0) ;
     }
     else //There is an image through the LE
     {
         //1 : case without taking that image
         double partialsum = (X1[0]-X2[0]) * (X1[0]-X2[0]) ;
         bitdim newgst = gst ; newgst &= (~(1<<30)) ;
         check_ghost_regular(newgst>>1, P, X1, X2, r1, r2, tmpcp, 1, partialsum, 0) ;
 
         //2: now is the hard case: we take the image path
         double Delta= (tmpcp.ghostdir&1?-1:1) * P.Boundaries[0].delta ;
         partialsum = (X1[0]-X2[0]-Delta) * (X1[0]-X2[0]-Delta) ;
         newgst = gst ;
         newgst &= (~(1<<1)) ;      // Clearing bit 1
         newgst |= ((gst>>30)<<1) ; // Setting bit 1 to the value of bit 30
         newgst &= (~(1<<30)) ;     // Clearing bit 30
         tmpcp.ghostdir &= (~(1<<1)) ; //Clearing bit 1
         tmpcp.ghostdir |= ((tmpcp.ghostdir>>30)<<1) ; // Setting the value of bit 1 to the value of bit 30. No need for clearing in the ghostdir.
         auto tmpX2 = X2 ;
         tmpX2[1] += (tmpcp.ghostdir&1?-1:1)*P.Boundaries[0].displacement ;
         if (tmpX2[1] > P.Boundaries[1].xmax) tmpX2[1] -= P.Boundaries[1].delta ;
         if (tmpX2[1] < P.Boundaries[1].xmin) tmpX2[1] += P.Boundaries[1].delta ;
         //printf("{%g %g %X %X", tmpX2[0], tmpX2[1], newgst>>1, tmpcp.ghostdir) ;
         check_ghost_regular(newgst>>1, P, X1, tmpX2, r1, r2, tmpcp, 1, partialsum, 1) ;
         //printf("}") ;
     }
    }
    
    return false ; 
}
 
//----------------------------------------
template <int d>
void ContactList<d>::check_ghost_dst(uint32_t gst, int n, double partialsum, uint32_t mask, const Parameters<d> & P, cv1d &X1, cv1d &X2, cp<d> & contact)
{
  if (gst==0) return ;
  else
  {
    for ( ;(gst&1)==0; gst>>=1,n++) ;
    check_ghost_dst(gst-1, n, partialsum, mask, P, X1, X2, contact) ;
    double Delta= (contact.ghostdir&(1<<n)?-1:1) * P.Boundaries[n].delta ;
    partialsum = partialsum + Delta*(2*(X2[n]-X1[n]) + Delta) ;
    if (partialsum<contact.contactlength) // Found a lower distance with this ghost
    {
      contact.contactlength=partialsum ;
      contact.ghost = mask|(1<<n) ;
    }
    check_ghost_dst(gst-1, n, partialsum, mask|(1<<n), P, X1, X2, contact) ;
  }
}
 
//-----------------------------------
template <int d>
void ContactList<d>::coordinance (v1d &Z)
{
  for (auto & w : v)
  {
      Z[w.i]++ ; Z[w.j] ++ ;
  }
}
 
#endif