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