1 /* ----------------------------------------------------------------------
2     This is the
3 
4     ██╗     ██╗ ██████╗  ██████╗  ██████╗ ██╗  ██╗████████╗███████╗
5     ██║     ██║██╔════╝ ██╔════╝ ██╔════╝ ██║  ██║╚══██╔══╝██╔════╝
6     ██║     ██║██║  ███╗██║  ███╗██║  ███╗███████║   ██║   ███████╗
7     ██║     ██║██║   ██║██║   ██║██║   ██║██╔══██║   ██║   ╚════██║
8     ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║  ██║   ██║   ███████║
9     ╚══════╝╚═╝ ╚═════╝  ╚═════╝  ╚═════╝ ╚═╝  ╚═╝   ╚═╝   ╚══════╝®
10 
11     DEM simulation engine, released by
12     DCS Computing Gmbh, Linz, Austria
13     http://www.dcs-computing.com, office@dcs-computing.com
14 
15     LIGGGHTS® is part of CFDEM®project:
16     http://www.liggghts.com | http://www.cfdem.com
17 
18     Core developer and main author:
19     Christoph Kloss, christoph.kloss@dcs-computing.com
20 
21     LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22     License, version 2 or later. It is distributed in the hope that it will
23     be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25     received a copy of the GNU General Public License along with LIGGGHTS®.
26     If not, see http://www.gnu.org/licenses . See also top-level README
27     and LICENSE files.
28 
29     LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30     the producer of the LIGGGHTS® software and the CFDEM®coupling software
31     See http://www.cfdem.com/terms-trademark-policy for details.
32 
33 -------------------------------------------------------------------------
34     Contributing author and copyright for this file:
35 
36     Christoph Kloss (DCS Computing GmbH, Linz)
37     Arno Mayrhofer (CFDEMresearch GmbH, Linz)
38 
39     Copyright 2012-     DCS Computing GmbH, Linz
40     Copyright 2009-2012 JKU Linz
41     Copyright 2016-     CFDEMresearch GmbH, Linz
42 ------------------------------------------------------------------------- */
43 
44 #ifndef LMP_CONTACT_HISTORY_MESH_I_H
45 #define LMP_CONTACT_HISTORY_MESH_I_H
46 
47   /* ---------------------------------------------------------------------- */
48 
49   /* ---------------------------------------------------------------------- */
50 
handleContact(int iP,int idTri,double * & history,bool intersect,bool faceflag)51   inline bool FixContactHistoryMesh::handleContact(int iP, int idTri, double *&history,bool intersect,bool faceflag)
52   {
53 
54     // check if contact with iTri was there before
55     // if so, set history to correct location and return
56     if(haveContact(iP,idTri,history,intersect))
57       return true;
58 
59     // else new contact - add contact if did not calculate contact with coplanar neighbor already
60 
61     if(faceflag && coplanarContactAlready(iP,idTri))
62     {
63         // did not add new contact
64         return false;
65     }
66     else
67     {
68         /*if (coplanarContactAlready(iP,idTri) && !faceflag)
69             {
70                 fprintf(screen,"WEIRD: atom ID %d, ts " BIGINT_FORMAT "\n",atom->tag[iP],update->ntimestep);
71             }
72         */
73         addNewTriContactToExistingParticle(iP,idTri,history,intersect);
74 
75         // check if one of the contacts of previous steps is coplanar with iTri
76 
77         // if so, copy history
78         if(faceflag)
79             checkCoplanarContactHistory(iP,idTri,history);
80         return true;
81     }
82   }
83 
84   /* ---------------------------------------------------------------------- */
85 
swap(int ilocal,int ineigh,int jneigh,bool keepflag_swap)86   inline void FixContactHistoryMesh::swap(int ilocal,int ineigh, int jneigh, bool keepflag_swap)
87   {
88 
89       int id_temp;
90 
91       id_temp                  = partner_[ilocal][ineigh];
92       partner_[ilocal][ineigh] = partner_[ilocal][jneigh];
93       partner_[ilocal][jneigh] = id_temp;
94 
95       vectorCopyN(&(contacthistory_[ilocal][ineigh*dnum_]),swap_,                                   dnum_);
96       vectorCopyN(&(contacthistory_[ilocal][jneigh*dnum_]),&(contacthistory_[ilocal][ineigh*dnum_]),dnum_);
97       vectorCopyN(swap_,                                   &(contacthistory_[ilocal][jneigh*dnum_]),dnum_);
98 
99       if(keepflag_swap)
100       {
101           const bool keepflag_temp  = keepflag_[ilocal][ineigh];
102           keepflag_[ilocal][ineigh] = keepflag_[ilocal][jneigh];
103           keepflag_[ilocal][jneigh] = keepflag_temp;
104 
105           const bool intersectflag_temp = intersectflag_[ilocal][ineigh];
106           intersectflag_[ilocal][ineigh] = intersectflag_[ilocal][jneigh];
107           intersectflag_[ilocal][jneigh] = intersectflag_temp;
108       }
109   }
110 
111   /* ---------------------------------------------------------------------- */
112 
haveContact(int iP,int idTri,double * & history,bool intersect)113   inline bool FixContactHistoryMesh::haveContact(int iP, int idTri, double *&history,bool intersect)
114   {
115     int *tri = partner_[iP];
116     const int nneighs = fix_nneighs_->get_vector_atom_int(iP);
117 
118     for(int i = 0; i < nneighs; i++)
119     {
120         if(tri[i] == idTri)
121         {
122             if(dnum_ > 0) history = &(contacthistory_[iP][i*dnum_]);
123             keepflag_[iP][i] = true;
124             intersectflag_[iP][i] = intersect;
125             return true;
126         }
127     }
128     return false;
129   }
130 
131   /* ---------------------------------------------------------------------- */
132 
coplanarContactAlready(int iP,int idTri)133   inline bool FixContactHistoryMesh::coplanarContactAlready(int iP, int idTri)
134   {
135     const int nneighs = fix_nneighs_->get_vector_atom_int(iP);
136     for(int i = 0; i < nneighs; i++)
137     {
138 
139       int idPartnerTri = partner_[iP][i];
140 
141       if(idPartnerTri >= 0 && idPartnerTri != idTri && mesh_->map(idPartnerTri, 0) >= 0 && mesh_->areCoplanarNodeNeighs(idPartnerTri,idTri))
142       {
143 
144         // other coplanar contact handled already - do not handle this contact
145         if(keepflag_[iP][i]) return true;
146       }
147     }
148 
149     // no coplanar contact found - handle this contact
150     return false;
151   }
152 
153   /* ---------------------------------------------------------------------- */
154 
checkCoplanarContactHistory(int iP,int idTri,double * & history)155   inline void FixContactHistoryMesh::checkCoplanarContactHistory(int iP, int idTri, double *&history)
156   {
157     int *tri = partner_[iP];
158     const int nneighs = fix_nneighs_->get_vector_atom_int(iP);
159 
160     for(int i = 0; i < nneighs; i++)
161     {
162 
163       if(tri[i] >= 0 && tri[i] != idTri && mesh_->map(tri[i], 0) >= 0 && mesh_->areCoplanarNodeNeighs(tri[i],idTri))
164       {
165 
166           // copy contact history
167           if(dnum_ > 0) vectorCopyN(&(contacthistory_[iP][i*dnum_]),history,dnum_);
168 
169       }
170     }
171   }
172 
173   /* ---------------------------------------------------------------------- */
174 
addNewTriContactToExistingParticle(int iP,int idTri,double * & history,bool intersect)175   inline void FixContactHistoryMesh::addNewTriContactToExistingParticle(int iP, int idTri, double *&history, bool intersect)
176   {
177 
178       const int nneighs = fix_nneighs_->get_vector_atom_int(iP);
179       int iContact = -1;
180 
181       if(-1 == idTri)
182         error->one(FLERR,"internal error");
183 
184       if(npartner_[iP] >= nneighs)
185       {
186 
187         error->one(FLERR,"internal error");
188       }
189 
190       for(int ineigh = 0; ineigh < nneighs; ineigh++)
191       {
192           if(-1 == partner_[iP][ineigh])
193           {
194               iContact = ineigh;
195               break;
196           }
197       }
198 
199       if(iContact >= nneighs)
200         error->one(FLERR,"internal error");
201 
202       partner_[iP][iContact] = idTri;
203       keepflag_[iP][iContact] = true;
204       intersectflag_[iP][iContact] = intersect;
205 
206       if(dnum_ > 0)
207       {
208           history = &(contacthistory_[iP][iContact*dnum_]);
209           vectorZeroizeN(history,dnum_);
210       }
211       else
212           history = 0;
213 
214       npartner_[iP]++;
215 
216   }
217 
218   /* ---------------------------------------------------------------------- */
219 
n_contacts(int & nIntersect)220   inline int FixContactHistoryMesh::n_contacts(int & nIntersect)
221   {
222     int ncontacts = 0, nlocal = atom->nlocal;
223 
224     for(int i = 0; i < nlocal; i++)
225     {
226         for(int ipartner = 0; ipartner < npartner_[i]; ipartner++)
227         {
228             if(intersectflag_[i][ipartner])
229             {
230                 nIntersect++;
231             }
232             ncontacts++;
233         }
234     }
235     return ncontacts;
236   }
237 
238   /* ---------------------------------------------------------------------- */
239 
n_contacts(int contact_groupbit,int & nIntersect)240   inline int FixContactHistoryMesh::n_contacts(int contact_groupbit, int & nIntersect)
241   {
242     int ncontacts = 0, nlocal = atom->nlocal;
243     int *mask = atom->mask;
244 
245     for(int i = 0; i < nlocal; i++)
246     {
247         if(mask[i] & contact_groupbit)
248         {
249             for(int ipartner = 0; ipartner < npartner_[i]; ipartner++)
250             {
251                 if(intersectflag_[i][ipartner])
252                 {
253                     nIntersect++;
254                 }
255                 ncontacts++;
256             }
257         }
258     }
259     return ncontacts;
260   }
261 
262 #endif
263