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     Christoph Kloss (JKU Linz)
38     Philippe Seil (JKU Linz)
39 
40     Copyright 2012-     DCS Computing GmbH, Linz
41     Copyright 2009-2012 JKU Linz
42 ------------------------------------------------------------------------- */
43 
44 #ifndef LMP_SURFACE_MESH_I_H
45 #define LMP_SURFACE_MESH_I_H
46 
47 #define NTRY_MC_SURFACE_MESH_I_H 30000
48 #define NITER_MC_SURFACE_MESH_I_H 5
49 #define TOLERANCE_MC_SURFACE_MESH_I_H 0.05
50 
51 /* ----------------------------------------------------------------------
52    constructors, destructor
53 ------------------------------------------------------------------------- */
54 
55 template<int NUM_NODES, int NUM_NEIGH_MAX>
SurfaceMesh(LAMMPS * lmp)56 SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::SurfaceMesh(LAMMPS *lmp)
57 :   TrackingMesh<NUM_NODES>(lmp),
58     curvature_(1.-EPSILON_CURVATURE),
59     curvature_tolerant_(false),
60 
61     minAngle_softLimit_(cos(0.5*M_PI/180.)),
62 
63     minAngle_hardLimit_(0.9999995),
64 
65     // TODO should keep areaMeshSubdomain up-to-date more often for insertion faces
66     areaMesh_     (*this->prop().template addGlobalProperty   < ScalarContainer<double> >                 ("areaMesh",     "comm_none","frame_trans_rot_invariant","restart_no",2)),
67 
68     nBelowAngle_softLimit_(0),
69     nTooManyNeighs_(0),
70     nOverlapping_(0),
71 
72     area_         (*this->prop().template addElementProperty< ScalarContainer<double> >                   ("area",         "comm_none","frame_trans_rot_invariant", "restart_no",2)),
73     areaAcc_      (*this->prop().template addElementProperty< ScalarContainer<double> >                   ("areaAcc",      "comm_none","frame_trans_rot_invariant", "restart_no",2)),
74     edgeLen_      (*this->prop().template addElementProperty< VectorContainer<double,NUM_NODES> >         ("edgeLen",      "comm_none","frame_trans_rot_invariant", "restart_no")),
75     edgeVec_      (*this->prop().template addElementProperty< MultiVectorContainer<double,NUM_NODES,3> >  ("edgeVec",      "comm_none","frame_scale_trans_invariant","restart_no")),
76     edgeNorm_     (*this->prop().template addElementProperty< MultiVectorContainer<double,NUM_NODES,3> >  ("edgeNorm",     "comm_none","frame_scale_trans_invariant","restart_no")),
77     surfaceNorm_  (*this->prop().template addElementProperty< VectorContainer<double,3> >                 ("surfaceNorm",  "comm_none","frame_scale_trans_invariant","restart_no")),
78     obtuseAngleIndex_   (*this->prop().template addElementProperty< ScalarContainer<int> >                ("obtuseAngleIndex","comm_exchange_borders","frame_invariant","restart_no")),
79     nNeighs_      (*this->prop().template addElementProperty< ScalarContainer<int> >                      ("nNeighs",      "comm_exchange_borders","frame_invariant","restart_no")),
80     neighFaces_   (*this->prop().template addElementProperty< VectorContainer<int,NUM_NEIGH_MAX> >        ("neighFaces",   "comm_exchange_borders","frame_invariant","restart_no")),
81     hasNonCoplanarSharedNode_(*this->prop().template addElementProperty< VectorContainer<bool,NUM_NODES> >("hasNonCoplanarSharedNode","comm_exchange_borders","frame_invariant", "restart_no")),
82     edgeActive_   (*this->prop().template addElementProperty< VectorContainer<bool,NUM_NODES> >           ("edgeActive",   "comm_exchange_borders","frame_invariant","restart_no")),
83     cornerActive_ (*this->prop().template addElementProperty< VectorContainer<bool,NUM_NODES> >           ("cornerActive", "comm_exchange_borders","frame_invariant","restart_no")),
84     neighList_(*new RegionNeighborList<interpolate_no>(lmp))
85 {
86 
87     areaMesh_.add(0.);
88     areaMesh_.add(0.);
89     areaMesh_.add(0.);
90 
91 }
92 
93 template<int NUM_NODES, int NUM_NEIGH_MAX>
~SurfaceMesh()94 SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::~SurfaceMesh()
95 {
96     delete &neighList_;
97 }
98 
99 /* ----------------------------------------------------------------------
100    set mesh curvature, used for mesh topology
101 ------------------------------------------------------------------------- */
102 
103 template<int NUM_NODES, int NUM_NEIGH_MAX>
setCurvature(double _curvature)104 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::setCurvature(double _curvature)
105 {
106     curvature_ = _curvature;
107 }
108 
109 /* ----------------------------------------------------------------------
110    set mesh curvature tolerance
111 ------------------------------------------------------------------------- */
112 
113 template<int NUM_NODES, int NUM_NEIGH_MAX>
setCurvatureTolerant(bool _tol)114 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::setCurvatureTolerant(bool _tol)
115 {
116     curvature_tolerant_ = _tol;
117 }
118 
119 /* ----------------------------------------------------------------------
120    add and delete an element
121 ------------------------------------------------------------------------- */
122 
123 template<int NUM_NODES, int NUM_NEIGH_MAX>
addElement(double ** nodeToAdd,int lineNumb)124 bool SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::addElement(double **nodeToAdd,int lineNumb)
125 {
126     if(TrackingMesh<NUM_NODES>::addElement(nodeToAdd,lineNumb))
127     {
128 
129         calcSurfPropertiesOfNewElement();
130         return true;
131     }
132     return false;
133 }
134 
135 template<int NUM_NODES, int NUM_NEIGH_MAX>
deleteElement(int n)136 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::deleteElement(int n)
137 {
138     TrackingMesh<NUM_NODES>::deleteElement(n);
139 }
140 
141 /* ----------------------------------------------------------------------
142    recalculate properties on setup (on start and during simulation)
143 ------------------------------------------------------------------------- */
144 
145 template<int NUM_NODES, int NUM_NEIGH_MAX>
refreshOwned(int setupFlag)146 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::refreshOwned(int setupFlag)
147 {
148     TrackingMesh<NUM_NODES>::refreshOwned(setupFlag);
149     // (re)calculate all properties for owned elements
150 
151     recalcLocalSurfProperties();
152 }
153 
154 template<int NUM_NODES, int NUM_NEIGH_MAX>
refreshGhosts(int setupFlag)155 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::refreshGhosts(int setupFlag)
156 {
157     TrackingMesh<NUM_NODES>::refreshGhosts(setupFlag);
158 
159     recalcGhostSurfProperties();
160 }
161 
162 /* ----------------------------------------------------------------------
163    recalculate properties of local elements
164 ------------------------------------------------------------------------- */
165 
166 template<int NUM_NODES, int NUM_NEIGH_MAX>
recalcLocalSurfProperties()167 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::recalcLocalSurfProperties()
168 {
169 
170     // areaMeshGlobal [areaMesh_(0)] and areaMeshOwned [areaMesh_(1)]
171     // calculated here
172 
173     areaMesh_(0) = 0.;
174     areaMesh_(1) = 0.;
175 
176     int nlocal = this->sizeLocal();
177 
178     for(int i = 0; i < nlocal; i++)
179     {
180       calcEdgeVecLen(i, edgeLen(i), edgeVec(i));
181       calcSurfaceNorm(i, surfaceNorm(i));
182       calcEdgeNormals(i, edgeNorm(i));
183       for(int j=0;j<NUM_NODES;j++)
184       {
185           double dot;
186           calcObtuseAngleIndex(i,j,dot);
187       }
188 
189       area(i) = calcArea(i);
190       areaAcc(i) = area(i);
191       if(i > 0) areaAcc(i) += areaAcc(i-1);
192 
193       // add to local area
194       areaMesh_(1) += area(i);
195 
196     }
197 
198     // mesh area must be summed up
199     MPI_Sum_Scalar(areaMesh_(1),areaMesh_(0),this->world);
200 
201 }
202 
203 /* ----------------------------------------------------------------------
204    recalculate properties of ghost elements
205 ------------------------------------------------------------------------- */
206 
207 template<int NUM_NODES, int NUM_NEIGH_MAX>
recalcGhostSurfProperties()208 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::recalcGhostSurfProperties()
209 {
210     int nlocal = this->sizeLocal();
211     int nall = this->sizeLocal()+this->sizeGhost();
212 
213     // areaMeshGhost [areaMesh_(2)] calculated here
214 
215     // accumulated area includes owned and ghosts
216     areaMesh_(2) = 0.;
217     for(int i = nlocal; i < nall; i++)
218     {
219       calcEdgeVecLen(i, edgeLen(i), edgeVec(i));
220       calcSurfaceNorm(i, surfaceNorm(i));
221       calcEdgeNormals(i, edgeNorm(i));
222 
223       for(int j=0;j<NUM_NODES;j++)
224       {
225           double dot;
226           calcObtuseAngleIndex(i,j,dot);
227       }
228 
229       area(i) = calcArea(i);
230       areaAcc(i) = area(i);
231       if(i > 0) areaAcc(i) += areaAcc(i-1);
232 
233       // add to ghost area
234       areaMesh_(2) += area(i);
235     }
236 
237 }
238 
239 /* ----------------------------------------------------------------------
240    generate a random Element by areaAcc
241 ------------------------------------------------------------------------- */
242 
243 template<int NUM_NODES, int NUM_NEIGH_MAX>
randomOwnedGhostElement()244 inline int SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::randomOwnedGhostElement()
245 {
246 
247     if(!this->isInsertionMesh())
248         this->error->one(FLERR,"Illegal call for non-insertion mesh");
249 
250     double area = areaMeshOwned()+areaMeshGhost();
251 
252     double r = this->random_->uniform() * area;
253 
254     int first = 0;
255     int last = this->sizeLocal()+this->sizeGhost()-1;
256 
257     return searchElementByAreaAcc(r,first,last);
258 }
259 
260 template<int NUM_NODES, int NUM_NEIGH_MAX>
searchElementByAreaAcc(double area,int lo,int hi)261 inline int SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::searchElementByAreaAcc(double area,int lo, int hi)
262 {
263 
264     if( (lo < 1 || area > areaAcc(lo-1)) && (area <= areaAcc(lo)) )
265         return lo;
266     if( (hi < 1 || area > areaAcc(hi-1)) && (area <= areaAcc(hi)) )
267         return hi;
268 
269     int mid = static_cast<int>((lo+hi)/2);
270     if(area > areaAcc(mid))
271         return searchElementByAreaAcc(area,mid,hi);
272     else
273         return searchElementByAreaAcc(area,lo,mid);
274 }
275 
276 /* ----------------------------------------------------------------------
277    calculate surface properties of new element
278    only called once on import
279 ------------------------------------------------------------------------- */
280 
281 template<int NUM_NODES, int NUM_NEIGH_MAX>
calcSurfPropertiesOfNewElement()282 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::calcSurfPropertiesOfNewElement()
283 {
284 
285     int n = this->sizeLocal()-1;
286 
287     double *vecTmp3,*vecTmpNumNodes,**nodeTmp;
288     create<double>(vecTmp3,3);
289     create<double>(vecTmpNumNodes,NUM_NODES);
290     create<double>(nodeTmp,NUM_NODES,3);
291 
292     // calculate edge vectors and lengths
293     calcEdgeVecLen(n,vecTmpNumNodes,nodeTmp);
294     edgeLen_.set(n,vecTmpNumNodes);
295     edgeVec_.set(n,nodeTmp);
296 
297     // calc surface normal
298     calcSurfaceNorm(n,vecTmp3);
299     surfaceNorm_.set(n,vecTmp3);
300 
301     // calc edge normal in plane pointing outwards of area_
302     // should be (edgeVec_ cross surfaceNormal)
303     calcEdgeNormals(n,nodeTmp);
304     edgeNorm_.set(n,nodeTmp);
305 
306     obtuseAngleIndex_.set(n,NO_OBTUSE_ANGLE);
307 
308     bool hasSmallAngleSoftLimit = false;
309     bool hasSmallAngleHardLimit = false;
310 
311     for(int i=0;i<NUM_NODES;i++){
312       double dot;
313       calcObtuseAngleIndex(n,i,dot);
314       if(-dot > minAngle_softLimit_)
315         hasSmallAngleSoftLimit = true;
316       if(-dot >= minAngle_hardLimit_)
317         hasSmallAngleHardLimit = true;
318     }
319 
320     if(hasSmallAngleSoftLimit)
321     {
322         if(TrackingMesh<NUM_NODES>::verbose() && 0 == this->comm->me)
323             fprintf(this->screen,"Mesh %s: element %d (line %d) has high aspect ratio (soft limit: smallest angle must be > %f °) \n",
324                     this->mesh_id_,n,TrackingMesh<NUM_NODES>::lineNo(n),this->angleSoftLimit());
325         nBelowAngle_softLimit_++;
326     }
327 
328     if(hasSmallAngleHardLimit && MultiNodeMesh<NUM_NODES>::elementExclusionList())
329     {
330 
331         if(0 == this->comm->me)
332         {
333 
334             fprintf(MultiNodeMesh<NUM_NODES>::elementExclusionList(),"%d\n",TrackingMesh<NUM_NODES>::lineNo(n));
335         }
336     }
337 
338     // calc area_ from previously obtained values and add to container
339     // calcArea is pure virtual and implemented in derived class(es)
340 
341     double area_elem = calcArea(n);
342     areaMesh_(0) += area_elem;
343     area_(n) = area_elem;
344     areaAcc_(n) = area_elem;
345     if(n > 0) areaAcc_(n) += areaAcc_(n-1);
346 
347     // cannot calc areaMesh_(1), areaMesh_(2), areaMesh_(3) here since
348     // not parallelized at this point
349 
350     destroy<double>(nodeTmp);
351     destroy<double>(vecTmpNumNodes);
352     destroy<double>(vecTmp3);
353 
354 }
355 
356 /* ----------------------------------------------------------------------
357    sub-functions needed to calculate mesh properties
358 ------------------------------------------------------------------------- */
359 
360 template<int NUM_NODES, int NUM_NEIGH_MAX>
calcEdgeVecLen(int nElem,double * len,double ** vec)361 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::calcEdgeVecLen(int nElem, double *len, double **vec)
362 {
363     for(int i=0;i<NUM_NODES;i++)
364     {
365       vectorSubtract3D(
366         MultiNodeMesh<NUM_NODES>::node_(nElem)[(i+1)%NUM_NODES],
367         MultiNodeMesh<NUM_NODES>::node_(nElem)[i],vec[i]);
368       len[i] = vectorMag3D(vec[i]);
369       vectorScalarDiv3D(vec[i],len[i]);
370     }
371 }
372 
373 template<int NUM_NODES, int NUM_NEIGH_MAX>
calcSurfaceNorm(int nElem,double * surfNorm)374 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::calcSurfaceNorm(int nElem, double *surfNorm)
375 {
376     vectorCross3D(edgeVec(nElem)[0],edgeVec(nElem)[1],surfNorm);
377 
378     if(vectorMag3D(surfNorm) <1e-15)
379     {
380 
381         vectorCross3D(edgeVec(nElem)[1],edgeVec(nElem)[2],surfNorm);
382 
383         if(vectorMag3D(surfNorm) <1e-15)
384         {
385             vectorCross3D(edgeVec(nElem)[2],edgeVec(nElem)[0],surfNorm);
386 
387             if(vectorMag3D(surfNorm) <1e-15)
388             {
389                 double tmpvec[] = {1.1233,2.123231,-3.3343434};
390                 vectorCross3D(edgeVec(nElem)[0],tmpvec,surfNorm);
391 
392                 if(vectorMag3D(surfNorm) <1e-15)
393                 {
394                     vectorCross3D(edgeVec(nElem)[1],tmpvec,surfNorm);
395 
396                     if(vectorMag3D(surfNorm) <1e-15)
397                     {
398                         double tmpvec2[] = {1.1233,-2.123231,3.3343434};
399                         vectorCross3D(edgeVec(nElem)[0],tmpvec2,surfNorm);
400 
401                         if(vectorMag3D(surfNorm) <1e-15)
402                         {
403                             vectorCross3D(edgeVec(nElem)[1],tmpvec2,surfNorm);
404 
405                         }
406                     }
407                 }
408             }
409         }
410 
411     }
412 
413     vectorScalarDiv3D(surfNorm, vectorMag3D(surfNorm));
414 }
415 
416 template<int NUM_NODES, int NUM_NEIGH_MAX>
calcEdgeNormals(int nElem,double ** edgeNorm)417 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::calcEdgeNormals(int nElem, double **edgeNorm)
418 {
419     for(int i=0;i<NUM_NODES;i++)
420     {
421       vectorCross3D(edgeVec(nElem)[i],surfaceNorm(nElem),edgeNorm[i]);
422 
423       if(vectorMag3D(edgeNorm[i]) <1e-15)
424       {
425           int otherIndex = (i+1)%3;
426           vectorCopy3D(edgeVec(nElem)[otherIndex],edgeNorm[i]);
427 
428       }
429       else
430         vectorScalarDiv3D(edgeNorm[i],vectorMag3D(edgeNorm[i]));
431     }
432 }
433 
434 template<int NUM_NODES, int NUM_NEIGH_MAX>
calcObtuseAngleIndex(int nElem,int iNode,double & dot)435 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::calcObtuseAngleIndex(int nElem, int iNode, double &dot)
436 {
437 
438     dot = vectorDot3D(edgeVec_(nElem)[iNode],edgeVec_(nElem)[(iNode-1+NUM_NODES)%NUM_NODES]);
439 
440     if(dot > 0.)
441     {
442 
443         obtuseAngleIndex_.set(nElem,iNode);
444     }
445     else
446         obtuseAngleIndex_.set(nElem,NO_OBTUSE_ANGLE);
447 }
448 
449 /* ----------------------------------------------------------------------
450    build neighlist, generate mesh topology, check (in)active edges and nodes
451 ------------------------------------------------------------------------- */
452 
453 template<int NUM_NODES, int NUM_NEIGH_MAX>
buildNeighbours()454 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::buildNeighbours()
455 {
456 
457     int nall = this->sizeLocal()+this->sizeGhost();
458 
459     if (this->lmp->wb && this->comm->me == 0)
460         fprintf(this->screen,"\nBuilding mesh topology (mesh processing step 2/3) \n");
461 
462     bool t[NUM_NODES], f[NUM_NODES];
463     int neighs[NUM_NEIGH_MAX];
464 
465     for(int i = 0; i < NUM_NODES; i++)
466     {
467         t[i] = true;
468         f[i] = false;
469     }
470     for(int i = 0; i < NUM_NEIGH_MAX; i++)
471         neighs[i] = -1;
472 
473     for(int i = 0; i < nall; i++)
474     {
475         nNeighs_.set(i,0);
476         neighFaces_.set(i,neighs);
477         edgeActive_.set(i,t);
478         cornerActive_.set(i,t);
479         hasNonCoplanarSharedNode_.set(i,f);
480     }
481 
482     // build neigh topology and edge activity
483 
484     BoundingBox bb(this->domain->boxlo[0], this->domain->boxhi[0],
485                    this->domain->boxlo[1], this->domain->boxhi[1],
486                    this->domain->boxlo[2], this->domain->boxhi[2]);
487 
488     neighList_.clear();
489 
490     double rBound_max = this->rBound_.max();
491     double binsize_factor = rBound_max ;
492 
493     if(nall > 100000 && binsize_factor > cbrt(bb.getVolume())/(4*20.))
494     {
495         binsize_factor = cbrt(bb.getVolume())/(4*20.);
496     }
497 
498     if(nall > 0 &&  neighList_.setBoundingBox(bb, binsize_factor, true, true))
499     {
500         std::vector<int> overlaps;
501 
502         for (int i = 0; i < nall; ++i)
503         {
504             //useless since would need allreduce to work
505             //if (this->lmp->wb && this->comm->me == 0 && 0 == i % 100000)
506             //    fprintf(this->screen,"   successfully built for a chunk of 100000 mesh elements\n");
507 
508             overlaps.clear();
509             neighList_.hasOverlapWith(this->center_(i), this->rBound_(i),overlaps);
510 
511             for(size_t iOverlap = 0; iOverlap < overlaps.size(); iOverlap++)
512             {
513                 int j = overlaps[iOverlap];
514                 if(j < 0 || j >= nall)
515                     this->error->one(FLERR,"Mesh error: internal error");
516 
517                 int iEdge(0), jEdge(0);
518 
519                 if(shareEdge(i,j,iEdge,jEdge))
520                   handleSharedEdge(i,iEdge,j,jEdge, areCoplanar(TrackingMesh<NUM_NODES>::id(i),TrackingMesh<NUM_NODES>::id(j)));
521             }
522 
523             neighList_.insert(this->center_(i), this->rBound_(i),i);
524         }
525     }
526     else if(nall > 0)
527         this->error->one(FLERR,"Mesh error: bounding box for neigh topology not set sucessfully");
528 
529     int *idListVisited = new int[nall];
530     int *idListHasNode = new int[nall];
531     double **edgeList,**edgeEndPoint;
532     this->memory->create(edgeList,2*nall,3,"SurfaceMesh:edgeList");
533     this->memory->create(edgeEndPoint,2*nall,3,"SurfaceMesh:edgeEndPoint");
534 
535     // recursively handle corner activity, ~n
536     for(int i = 0; i < nall; i++)
537     {
538         for(int iNode = 0; iNode < NUM_NODES; iNode++)
539             handleCorner(i,iNode,idListVisited,idListHasNode,edgeList,edgeEndPoint);
540     }
541 
542     if(MultiNodeMesh<NUM_NODES>::minFeatureLength() > 0. && MultiNodeMesh<NUM_NODES>::elementExclusionList())
543         handleExclusion(idListVisited);
544 
545     delete []idListVisited;
546     delete []idListHasNode;
547     this->memory->destroy(edgeList);
548     this->memory->destroy(edgeEndPoint);
549 
550     fflush(MultiNodeMesh<NUM_NODES>::elementExclusionList());
551 
552     // correct edge and corner activation/deactivation and neighs in parallel
553 
554     parallelCorrectionActiveInactive();
555     parallelCorrectionNeighs();
556 
557 }
558 
559 /* ----------------------------------------------------------------------
560    quality check for surface mesh
561 ------------------------------------------------------------------------- */
562 
563 template<int NUM_NODES, int NUM_NEIGH_MAX>
qualityCheck()564 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::qualityCheck()
565 {
566     if (this->lmp->wb && this->comm->me == 0)
567         fprintf(this->screen,"\nChecking quality of mesh (mesh processing step 3/3) \n");
568 
569     // iterate over surfaces
570 
571     int nlocal = this->sizeLocal();
572     int nall = this->sizeLocal()+this->sizeGhost();
573     int me = this->comm->me;
574 
575     // check duplicate elements, O(n) operation
576 
577     BoundingBox bb(this->domain->boxlo[0], this->domain->boxhi[0],
578                    this->domain->boxlo[1], this->domain->boxhi[1],
579                    this->domain->boxlo[2], this->domain->boxhi[2]);
580 
581     neighList_.clear();
582 
583     double rBound_max = this->rBound_.max();
584     double binsize_factor = rBound_max;
585 
586     if(nall > 100000 && binsize_factor > cbrt(bb.getVolume())/(4*20.))
587     {
588         binsize_factor = cbrt(bb.getVolume())/(4*20.);
589     }
590 
591     if(nall > 0 &&  neighList_.setBoundingBox(bb, binsize_factor, true,true))
592     {
593         std::vector<int> overlaps;
594 
595         for (int i = 0; i < nall; ++i)
596         {
597             //useless since would need allreduce to work
598             //if (this->lmp->wb && this->comm->me == 0 && 0 == i % 100000)
599             //    fprintf(this->screen,"   successfully checked a chunk of 100000 mesh elements\n");
600 
601             overlaps.clear();
602             neighList_.hasOverlapWith(this->center_(i), this->rBound_(i),overlaps);
603 
604             for(size_t iOverlap = 0; iOverlap < overlaps.size(); iOverlap++)
605             {
606                 int j = overlaps[iOverlap];
607                 if(j < 0 || j >= nall)
608                     this->error->one(FLERR,"Mesh error: internal error");
609 
610                 if(this->nSharedNodes(i,j) == NUM_NODES)
611                 {
612                     fprintf(this->screen,"ERROR: Mesh %s: elements %d and %d (lines %d and %d) are duplicate\n",
613                             this->mesh_id_,TrackingMesh<NUM_NODES>::id(i),TrackingMesh<NUM_NODES>::id(j),
614                             TrackingMesh<NUM_NODES>::lineNo(i),TrackingMesh<NUM_NODES>::lineNo(j));
615 
616                     if(MultiNodeMesh<NUM_NODES>::elementExclusionList())
617                         fprintf(MultiNodeMesh<NUM_NODES>::elementExclusionList(),"%d\n",TrackingMesh<NUM_NODES>::lineNo(j));
618                     else if(!this->removeDuplicates())
619                         this->error->one(FLERR,"Fix mesh: Bad mesh, cannot continue. You can try re-running with 'heal auto_remove_duplicates'");
620                     else
621                         this->error->one(FLERR,"Fix mesh: Bad mesh, cannot continue. The mesh probably reached the precision you defined. "
622                                                "You can try re-running with a lower value for 'precision'");
623                 }
624             }
625 
626             neighList_.insert(this->center_(i), this->rBound_(i),i);
627         }
628     }
629     else if(nall > 0)
630         this->error->one(FLERR,"Mesh error: bounding box for neigh topology not set sucessfully");
631 
632     fflush(MultiNodeMesh<NUM_NODES>::elementExclusionList());
633 
634     for(int i = 0; i < nlocal; i++)
635     {
636       for(int iNode = 0; iNode < NUM_NODES; iNode++)
637       {
638         double dot;
639         calcObtuseAngleIndex(i,iNode,dot);
640         if(-dot > curvature_)
641         {
642             if(TrackingMesh<NUM_NODES>::verbose() || !curvature_tolerant_)
643                 fprintf(this->screen,"%s: Mesh %s: The minumum angle of mesh element %d (line %d) is lower than the specified curvature. "
644                                   "Increase mesh quality or decrease curvature (currently %f°)\n",
645                                   curvature_tolerant_?"WARNING:":"ERROR",this->mesh_id_,TrackingMesh<NUM_NODES>::id(i),
646                                   TrackingMesh<NUM_NODES>::lineNo(i),acos(curvature_)*180./M_PI);
647             if(!curvature_tolerant_)
648                 this->error->one(FLERR,"Fix mesh: Bad mesh, cannot continue. You can try setting 'curvature' to 1e-5 or lower or use 'curvature_tolerant yes'");
649         }
650       }
651     }
652 
653     if(this->nBelowAngleSoftLimit() > 0 && 0 == me)
654     {
655         fprintf(this->screen,"Mesh %s: %d elements have high aspect ratio (soft limit: smallest angle > %f ° required)\n",
656                 this->mesh_id_,this->nBelowAngleSoftLimit(),this->angleSoftLimit());
657         this->error->warning(FLERR,"Fix mesh: Mesh contains highly skewed element, moving mesh (if used) will not parallelize well");
658     }
659 
660     int nBelowAngle_hardLimit = 0;
661     for(int i = 0; i < nlocal; i++)
662     {
663       for(int iNode = 0; iNode < NUM_NODES; iNode++)
664       {
665         double dot;
666         calcObtuseAngleIndex(i,iNode,dot);
667 
668         if(-dot > minAngle_hardLimit_)
669         {
670              if(TrackingMesh<NUM_NODES>::verbose())
671                 fprintf(this->screen,"Mesh %s: element %d (line %d) has a really unreasonably high aspect ratio (hard limit: smallest angle must be > %f °) \n",
672                         this->mesh_id_,TrackingMesh<NUM_NODES>::id(i),TrackingMesh<NUM_NODES>::lineNo(i),this->angleHardLimit());
673              nBelowAngle_hardLimit++;
674         }
675       }
676     }
677 
678     MPI_Sum_Scalar(nBelowAngle_hardLimit,this->world);
679     if(nBelowAngle_hardLimit > 0 && 0 == me)
680     {
681         fprintf(this->screen,"Mesh %s: %d mesh elements have a really unreasonably high aspect ratio (hard limit: smallest angle must be > %f °)  \n",
682                 this->mesh_id_,nBelowAngle_hardLimit,this->angleHardLimit());
683         this->error->one(FLERR,"Fix mesh: Bad mesh, cannot continue. You will need to fix or remove these element. Remedies:\n"
684                                 " - You can use the 'element_exclusion_list' feature to remove elements with too many neighbors and elements with angles below the hard limit\n");
685     }
686 
687     if(this->nTooManyNeighs() > 0 && 0 == me)
688     {
689 
690         fprintf(this->screen,"Mesh %s: %d mesh elements have more than %d neighbors \n",
691                 this->mesh_id_,this->nTooManyNeighs(),NUM_NEIGH_MAX);
692         this->error->one(FLERR,"Fix mesh: Bad mesh, cannot continue. Possibly corrupt elements with too many neighbors. Remedies:\n"
693                                " - You can use the 'element_exclusion_list' feature to remove elements with too many neighbors and elements with angles below the hard limit\n"
694                                " -  If you know what you're doing, you can alternatively also try to change the definition of SurfaceMeshBase in tri_mesh.h and recompile");
695     }
696 
697     if(nOverlapping() > 0)
698     {
699         fprintf(this->screen,"WARNING: Mesh %s: proc %d has %d element pairs that are coplanar, "
700                 "share an edge and overlap (but are not duplicate)\n",
701                 this->mesh_id_,me,nOverlapping());
702     }
703 }
704 
705 /* ----------------------------------------------------------------------
706    correct edge and corner activation/deactivation in parallel
707 ------------------------------------------------------------------------- */
708 
709 template<int NUM_NODES, int NUM_NEIGH_MAX>
parallelCorrectionActiveInactive()710 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::parallelCorrectionActiveInactive()
711 {
712 
713     int iGlobal;
714     int mine = this->sizeLocal()+this->sizeGhost();
715     int sizeGlob = this->sizeGlobal();
716     int len = NUM_NODES*sizeGlob;
717 
718     int *edgea = new int[len];
719     int *cornera = new int[len];
720     vectorInitializeN(edgea,len,2);
721     vectorInitializeN(cornera,len,2);
722 
723     for(int i = 0; i < mine; i++)
724     {
725         iGlobal = TrackingMesh<NUM_NODES>::id(i);
726 
727         for(int j = 0; j < NUM_NODES; j++)
728         {
729             edgea[iGlobal*NUM_NODES+j] = (edgeActive(i)[j] && edgea[iGlobal*NUM_NODES+j] > 0) ? 1 : 0;
730             cornera[iGlobal*NUM_NODES+j] = (cornerActive(i)[j] && cornera[iGlobal*NUM_NODES+j] > 0) ? 1 : 0;
731         }
732     }
733 
734     MPI_Min_Vector(edgea,len,this->world);
735     MPI_Min_Vector(cornera,len,this->world);
736 
737     for(int i = 0; i < sizeGlob; i++)
738     {
739         const int nTri_j = this->map_size(i);
740         for (int j = 0; j < nTri_j; j++)
741         {
742             const int iLocal = this->map(i, j);
743             if(iLocal >= 0)
744             {
745                 for(int j = 0; j < NUM_NODES; j++)
746                 {
747                     if(edgea[i*NUM_NODES+j] == 0)
748                         edgeActive(iLocal)[j] = false;
749                     else if(edgea[i*NUM_NODES+j] == 1)
750                         edgeActive(iLocal)[j] = true;
751                     else
752                         this->error->one(FLERR,"Illegal situation in SurfaceMesh::parallelCorrection()");
753                     if(cornera[i*NUM_NODES+j] == 0)
754                         cornerActive(iLocal)[j] = false;
755                     else if(cornera[i*NUM_NODES+j] == 1)
756                         cornerActive(iLocal)[j] = true;
757                     else
758                         this->error->one(FLERR,"Illegal situation in SurfaceMesh::parallelCorrection()");
759                 }
760             }
761         }
762     }
763 
764     delete []edgea;
765     delete []cornera;
766 }
767 
768 /* ----------------------------------------------------------------------
769    correct neighs in parallel
770 ------------------------------------------------------------------------- */
771 
772 template<int NUM_NODES, int NUM_NEIGH_MAX>
parallelCorrectionNeighs()773 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::parallelCorrectionNeighs()
774 {
775 
776     int iGlobal;
777     int nlocal = this->sizeLocal();
778     int nghost = this->sizeGhost();
779     int sizeGlob = this->sizeGlobal();
780     int len1 = NUM_NEIGH_MAX*sizeGlob;
781     int len2 = sizeGlob;
782 
783     int *neighs_found_by_owned_id = new int[len1];
784     vectorInitializeN(neighs_found_by_owned_id,len1,-1);
785 
786     int *additionalNeigh_id = new int[len2];
787     vectorInitializeN(additionalNeigh_id,len2,-1);
788 
789     for(int i = 0; i < nlocal; i++)
790     {
791         iGlobal = TrackingMesh<NUM_NODES>::id(i);
792 
793         for(int iNeigh = 0; iNeigh< nNeighs_(i); iNeigh++)
794         {
795             neighs_found_by_owned_id[iGlobal*NUM_NEIGH_MAX+iNeigh] = neighFaces_(i)[iNeigh];
796         }
797     }
798 
799     MPI_Max_Vector(neighs_found_by_owned_id,len1,this->world);
800 
801     for(int i = nlocal; i < nlocal+nghost; i++)
802     {
803         iGlobal = TrackingMesh<NUM_NODES>::id(i);
804 
805         for(int iNeigh = 0; iNeigh < nNeighs_(i); iNeigh++)
806         {
807             bool already_found = false;
808 
809             for(int iFound = 0; iFound < NUM_NEIGH_MAX; iFound++)
810             {
811                 if(neighs_found_by_owned_id[iGlobal*NUM_NEIGH_MAX+iFound] == neighFaces_(i)[iNeigh])
812                     already_found = true;
813             }
814 
815             if(!already_found)
816                 additionalNeigh_id[iGlobal] = neighFaces_(i)[iNeigh];
817         }
818     }
819 
820     MPI_Max_Vector(additionalNeigh_id,len2,this->world);
821 
822     for(int i = 0; i < nlocal; i++)
823     {
824         iGlobal = TrackingMesh<NUM_NODES>::id(i);
825 
826         if(additionalNeigh_id[iGlobal] > -1)
827         {
828             // set neighbor topology
829             if(nNeighs_(i) < NUM_NEIGH_MAX)
830                 neighFaces_(i)[nNeighs_(i)] = additionalNeigh_id[iGlobal] ;
831             nNeighs_(i)++;
832         }
833     }
834 
835     delete []neighs_found_by_owned_id;
836     delete []additionalNeigh_id;
837 }
838 
839 /* ----------------------------------------------------------------------
840    functions to generate mesh topology
841 ------------------------------------------------------------------------- */
842 
843 template<int NUM_NODES, int NUM_NEIGH_MAX>
areCoplanar(int tag_a,int tag_b)844 bool SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::areCoplanar(int tag_a, int tag_b)
845 {
846     int a = this->map(tag_a, 0);
847     int b = this->map(tag_b, 0);
848 
849     if(a < 0 || b < 0)
850         this->error->one(FLERR,"Internal error: Illegal call to SurfaceMesh::areCoplanar()");
851 
852     // check if two faces are coplanar
853 
854     double dot = vectorDot3D(surfaceNorm(a),surfaceNorm(b));
855 
856     // need fabs in case surface normal is other direction
857     if(fabs(dot) >= curvature_) return true;
858     else return false;
859 }
860 
861 /* ---------------------------------------------------------------------- */
862 
863 template<int NUM_NODES, int NUM_NEIGH_MAX>
areCoplanarNeighs(int tag_a,int tag_b)864 bool SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::areCoplanarNeighs(int tag_a, int tag_b)
865 {
866     bool areNeighs = false;
867     int a = this->map(tag_a, 0);
868     int b = this->map(tag_b, 0);
869 
870     if(a < 0 || b < 0)
871         this->error->one(FLERR,"Internal error: Illegal call to SurfaceMesh::areCoplanarNeighs()");
872 
873     // check if two faces are coplanar
874 
875     // must be neighs, otherwise not considered coplanar
876     for(int i = 0; i < nNeighs_(a); i++)
877         if(neighFaces_(a)[i] == tag_b)
878             areNeighs = true;
879 
880     if(!areNeighs) return false;
881 
882     double dot = vectorDot3D(surfaceNorm(a),surfaceNorm(b));
883 
884     // need fabs in case surface normal is other direction
885     if(fabs(dot) > curvature_) return true;
886     else return false;
887 }
888 
889 /* ---------------------------------------------------------------------- */
890 
891 template<int NUM_NODES, int NUM_NEIGH_MAX>
areCoplanarNodeNeighs(int tag_a,int tag_b)892 bool SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::areCoplanarNodeNeighs(int tag_a, int tag_b)
893 {
894     bool areNeighs = false;
895     int a = this->map(tag_a, 0);
896     int b = this->map(tag_b, 0);
897 
898     if(a < 0 || b < 0)
899         this->error->one(FLERR,"Internal error: Illegal call to SurfaceMesh::areCoplanarNeighs()");
900 
901     // check if two faces are coplanar
902 
903     // must be neighs, otherwise not considered coplanar
904     for(int i = 0; i < nNeighs_(a); i++)
905         if(neighFaces_(a)[i] == tag_b)
906             areNeighs = true;
907 
908     const int nTri_j = this->map_size(tag_b);
909     bool found = false;
910     // only check if normals are equal if they are not listed as neigbors
911     if (!areNeighs)
912     {
913         for (int j = 0; j < nTri_j; j++)
914         {
915             const int b_tmp = this->map(tag_b, j);
916             if (MultiNodeMesh<NUM_NODES>::nSharedNodes(a,b_tmp) != 0)
917             {
918                 found = true;
919                 break;
920             }
921         }
922     }
923     if(!areNeighs && !found) return false;
924 
925     double dot = vectorDot3D(surfaceNorm(a),surfaceNorm(b));
926 
927     // need fabs in case surface normal is other direction
928     if(fabs(dot) > curvature_) return true;
929     else return false;
930 }
931 
932 /* ---------------------------------------------------------------------- */
933 
934 template<int NUM_NODES, int NUM_NEIGH_MAX>
coplanarNeighsOverlap(int iSrf,int iEdge,int jSrf,int jEdge)935 bool SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::coplanarNeighsOverlap(int iSrf,int iEdge,int jSrf,int jEdge)
936 {
937 
938     double vecI[3],vecJ[3], pRef[3], edgeN[3], dot1, dot2;
939 
940     vectorCopy3D(MultiNodeMesh<NUM_NODES>::node_(iSrf)[iEdge],pRef);
941     vectorCopy3D(edgeNorm(iSrf)[iEdge],edgeN);
942 
943     vectorSubtract3D(MultiNodeMesh<NUM_NODES>::node_(iSrf)[(iEdge+2)%NUM_NODES],pRef,vecI);
944     vectorSubtract3D(MultiNodeMesh<NUM_NODES>::node_(jSrf)[(jEdge+2)%NUM_NODES],pRef,vecJ);
945 
946     dot1 = vectorDot3D(vecI,edgeN);
947     dot2 = vectorDot3D(vecJ,edgeN);
948 
949     if(dot1*dot2 > 0.)
950     {
951         if(TrackingMesh<NUM_NODES>::verbose())
952         {
953 
954             int nlocal = this->sizeLocal();
955             fprintf(this->screen,"WARNING: Mesh %s: elements %d and %d are coplanar, "
956                     "share an edge and overlap (but are not duplicate)\n",
957                     this->mesh_id_,TrackingMesh<NUM_NODES>::id(iSrf),TrackingMesh<NUM_NODES>::id(jSrf));
958             if(iSrf < nlocal)
959                 fprintf(this->screen,"INFO: Mesh %s: element %d corresponds to line # %d\n",
960                     this->mesh_id_,TrackingMesh<NUM_NODES>::id(iSrf),TrackingMesh<NUM_NODES>::lineNo(iSrf));
961             if(jSrf < nlocal)
962                 fprintf(this->screen,"INFO: Mesh %s: element %d corresponds to line # %d\n",
963                     this->mesh_id_,TrackingMesh<NUM_NODES>::id(jSrf),TrackingMesh<NUM_NODES>::lineNo(jSrf));
964         }
965 
966         nOverlapping_++;
967 
968         //this->error->warning(FLERR,"Fix mesh: Check overlapping mesh elements");
969         return true;
970     }
971     else return false;
972 }
973 
974 /* ---------------------------------------------------------------------- */
975 
976 template<int NUM_NODES, int NUM_NEIGH_MAX>
edgeVecsColinear(double * v,double * w)977 bool SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::edgeVecsColinear(double *v,double *w)
978 {
979     // need normalized vectors
980     double dot = vectorDot3D(v,w);
981     // need fabs in case vectors are in different direction
982     if(fabs(dot) > curvature_) return true;
983     else return false;
984 }
985 
986 /* ---------------------------------------------------------------------- */
987 
988 template<int NUM_NODES, int NUM_NEIGH_MAX>
growSurface(int iSrf,double by)989 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::growSurface(int iSrf, double by)
990 {
991     double *tmp = new double[3];
992     for(int i=0;i<NUM_NODES;i++)
993     {
994       vectorSubtract3D(MultiNodeMesh<NUM_NODES>::node(iSrf)[i],this->center_(iSrf),tmp);
995       vectorScalarMult3D(tmp,by);
996       vectorAdd3D(MultiNodeMesh<NUM_NODES>::node(iSrf)[i],
997                       tmp,MultiNodeMesh<NUM_NODES>::node(iSrf)[i]);
998     }
999     delete[] tmp;
1000     return;
1001 }
1002 
1003 /* ---------------------------------------------------------------------- */
1004 
1005 template<int NUM_NODES, int NUM_NEIGH_MAX>
shareEdge(int iSrf,int jSrf,int & iEdge,int & jEdge)1006 bool SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::shareEdge(int iSrf, int jSrf, int &iEdge, int &jEdge)
1007 {
1008     int iNode1=0,jNode1=0,iNode2,jNode2;
1009     if(this->share2Nodes(iSrf,jSrf,iNode1,jNode1,iNode2,jNode2)){
1010       // following implementation of shareNode(), the only remaining option to
1011       // share an edge is that the next node of iSrf is equal to the next or previous
1012       // node if jSrf
1013 
1014       if(2 == iNode1+iNode2)
1015         iEdge = 2;
1016 
1017       else
1018         iEdge = std::min(iNode1,iNode2);
1019 
1020       if(2 == jNode1+jNode2)
1021         jEdge = 2;
1022       else
1023         jEdge = std::min(jNode1,jNode2);
1024 
1025       return true;
1026     }
1027 
1028     iEdge = -1; jEdge = -1;
1029     return false;
1030 }
1031 
1032 /* ---------------------------------------------------------------------- */
1033 
1034 template<int NUM_NODES, int NUM_NEIGH_MAX>
handleSharedEdge(int iSrf,int iEdge,int jSrf,int jEdge,bool coplanar,bool neighflag)1035 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::handleSharedEdge(int iSrf, int iEdge, int jSrf, int jEdge,
1036                                             bool coplanar, bool neighflag)
1037 {
1038 
1039     if(neighflag)
1040     {
1041 
1042         if(nNeighs_(iSrf) == NUM_NEIGH_MAX)
1043         {
1044             nTooManyNeighs_++;
1045             if(TrackingMesh<NUM_NODES>::verbose())
1046                 fprintf(this->screen,"Mesh %s: element id %d (line %d) has %d neighs, but only %d expected\n",
1047                         this->mesh_id_,TrackingMesh<NUM_NODES>::id(iSrf),TrackingMesh<NUM_NODES>::lineNo(iSrf),nNeighs_(iSrf)+1,NUM_NEIGH_MAX);
1048             if(MultiNodeMesh<NUM_NODES>::elementExclusionList())
1049             {
1050 
1051                 fprintf(MultiNodeMesh<NUM_NODES>::elementExclusionList(),"%d\n",TrackingMesh<NUM_NODES>::lineNo(iSrf));
1052             }
1053 
1054         }
1055         if(nNeighs_(jSrf) == NUM_NEIGH_MAX)
1056         {
1057             nTooManyNeighs_++;
1058             if(TrackingMesh<NUM_NODES>::verbose())
1059                 fprintf(this->screen,"Mesh %s: element id %d (line %d) has %d neighs, but only %d expected\n",
1060                         this->mesh_id_,TrackingMesh<NUM_NODES>::id(jSrf),TrackingMesh<NUM_NODES>::lineNo(jSrf),nNeighs_(jSrf)+1,NUM_NEIGH_MAX);
1061             if(MultiNodeMesh<NUM_NODES>::elementExclusionList())
1062             {
1063 
1064                 fprintf(MultiNodeMesh<NUM_NODES>::elementExclusionList(),"%d\n",TrackingMesh<NUM_NODES>::lineNo(jSrf));
1065             }
1066 
1067         }
1068 
1069         // set neighbor topology
1070         if(nNeighs_(iSrf) < NUM_NEIGH_MAX)
1071             neighFaces_(iSrf)[nNeighs_(iSrf)] = TrackingMesh<NUM_NODES>::id(jSrf);
1072         if(nNeighs_(jSrf) < NUM_NEIGH_MAX)
1073             neighFaces_(jSrf)[nNeighs_(jSrf)] = TrackingMesh<NUM_NODES>::id(iSrf);
1074         nNeighs_(iSrf)++;
1075         nNeighs_(jSrf)++;
1076 
1077     }
1078 
1079     // deactivate one egde
1080     // other as well if coplanar
1081 
1082     if(!coplanar || coplanarNeighsOverlap(iSrf,iEdge,jSrf,jEdge))
1083     {
1084         if(TrackingMesh<NUM_NODES>::id(iSrf) < TrackingMesh<NUM_NODES>::id(jSrf))
1085         {
1086 
1087             edgeActive(iSrf)[iEdge] = false;
1088             edgeActive(jSrf)[jEdge] = true;
1089         }
1090         else
1091         {
1092 
1093             edgeActive(iSrf)[iEdge] = true;
1094             edgeActive(jSrf)[jEdge] = false;
1095         }
1096     }
1097     else // coplanar
1098     {
1099         if(!coplanar) this->error->one(FLERR,"internal error");
1100 
1101         edgeActive(iSrf)[iEdge] = false;
1102         edgeActive(jSrf)[jEdge] = false;
1103     }
1104 }
1105 
1106 /* ---------------------------------------------------------------------- */
1107 
1108 template<int NUM_NODES, int NUM_NEIGH_MAX>
handleCorner(int iSrf,int iNode,int * idListVisited,int * idListHasNode,double ** edgeList,double ** edgeEndPoint)1109 int SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::handleCorner(int iSrf, int iNode,
1110         int *idListVisited,int *idListHasNode,double **edgeList,double **edgeEndPoint)
1111 {
1112     double nodeToCheck[3];
1113     bool hasTwoColinearEdges, anyActiveEdge;
1114     int nIdListVisited = 0, nIdListHasNode = 0, maxId = -1, nEdgeList;
1115 
1116     this->node(iSrf,iNode,nodeToCheck);
1117     anyActiveEdge = false;
1118     checkNodeRecursive(iSrf,nodeToCheck,nIdListVisited,idListVisited,
1119         nIdListHasNode,idListHasNode,edgeList,edgeEndPoint,anyActiveEdge);
1120 
1121     if (!this->domain->is_in_subdomain(nodeToCheck))
1122         return nIdListHasNode;
1123 
1124     // each element that shares the node contributes two edges
1125     nEdgeList = 2*nIdListHasNode;
1126 
1127     // get max ID
1128     for(int i = 0; i < nIdListHasNode; i++)
1129         maxId = std::max(maxId,idListHasNode[i]);
1130 
1131     // check if any 2 edges coplanar
1132 
1133     hasTwoColinearEdges = false;
1134     for(int i = 0; i < nEdgeList; i++)
1135     {
1136         for(int j = i+1; j < nEdgeList; j++)
1137         {
1138 
1139             if(edgeVecsColinear(edgeList[i],edgeList[j]) && !this->nodesAreEqual(edgeEndPoint[i],edgeEndPoint[j]))
1140                 hasTwoColinearEdges = true;
1141         }
1142     }
1143 
1144     if(hasTwoColinearEdges || !anyActiveEdge)
1145         cornerActive(iSrf)[iNode] = false;
1146 
1147     else if(TrackingMesh<NUM_NODES>::id(iSrf) == maxId)
1148         cornerActive(iSrf)[iNode] = true;
1149     else
1150         cornerActive(iSrf)[iNode] = false;
1151 
1152     return nIdListHasNode;
1153 }
1154 
1155 /* ---------------------------------------------------------------------- */
1156 
1157 template<int NUM_NODES, int NUM_NEIGH_MAX>
checkNodeRecursive(int iSrf,double * nodeToCheck,int & nIdListVisited,int * idListVisited,int & nIdListHasNode,int * idListHasNode,double ** edgeList,double ** edgeEndPoint,bool & anyActiveEdge)1158 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::checkNodeRecursive(int iSrf,double *nodeToCheck,
1159         int &nIdListVisited,int *idListVisited,int &nIdListHasNode,int *idListHasNode,
1160         double **edgeList,double **edgeEndPoint,bool &anyActiveEdge)
1161 {
1162     int idNeigh, iNeigh, nEdgeList = 2*nIdListHasNode, nEdgeEndPoint = 2*nIdListHasNode;
1163 
1164     // check if I have been here already
1165     for(int i = 0; i < nIdListVisited; i++)
1166         if(idListVisited[i] == TrackingMesh<NUM_NODES>::id(iSrf)) return;
1167 
1168     // add to visited list
1169     idListVisited[nIdListVisited++] = TrackingMesh<NUM_NODES>::id(iSrf);
1170 
1171     // if contains node, add to list and call neighbors
1172     int iNode = this->containsNode(iSrf, nodeToCheck);
1173     if(iNode >= 0)
1174     {
1175 
1176         idListHasNode[nIdListHasNode++] = TrackingMesh<NUM_NODES>::id(iSrf);
1177         // node iNode is associated with edge iNode and iNode-1
1178         vectorCopy3D(edgeVec(iSrf)[iNode],edgeList[nEdgeList++]);
1179         vectorCopy3D(edgeVec(iSrf)[(iNode-1+NUM_NODES)%NUM_NODES],edgeList[nEdgeList++]);
1180         vectorCopy3D(this->node_(iSrf)[(iNode+1)%NUM_NODES],edgeEndPoint[nEdgeEndPoint++]);
1181         vectorCopy3D(this->node_(iSrf)[(iNode-1+NUM_NODES)%NUM_NODES],edgeEndPoint[nEdgeEndPoint++]);
1182         if(edgeActive(iSrf)[iNode]) anyActiveEdge = true;
1183         else if(edgeActive(iSrf)[(iNode-1+NUM_NODES)%NUM_NODES]) anyActiveEdge = true;
1184 
1185         // only call recursive if have neighbor and if I have data of neigh element (own or ghost)
1186 
1187         for(int iN = 0; iN < nNeighs_(iSrf); iN++)
1188         {
1189             idNeigh = neighFaces_(iSrf)[iN];
1190             if(idNeigh < 0) return;
1191             const int nTri_j = this->map_size(idNeigh);
1192             for (int j = 0; j < nTri_j; j++)
1193             {
1194                 iNeigh = this->map(idNeigh, j);
1195                 if(iNeigh >= 0)
1196                     checkNodeRecursive(iNeigh,nodeToCheck,nIdListVisited,idListVisited,nIdListHasNode,
1197                                        idListHasNode,edgeList,edgeEndPoint,anyActiveEdge);
1198             }
1199         }
1200     }
1201 
1202 }
1203 
1204 /* ----------------------------------------------------------------------
1205    move mesh
1206 ------------------------------------------------------------------------- */
1207 
1208 template<int NUM_NODES, int NUM_NEIGH_MAX>
move(const double * const vecTotal,const double * const vecIncremental)1209 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::move(const double * const vecTotal, const double * const vecIncremental)
1210 {
1211     TrackingMesh<NUM_NODES>::move(vecTotal,vecIncremental);
1212 }
1213 
1214 template<int NUM_NODES, int NUM_NEIGH_MAX>
move(const double * const vecIncremental)1215 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::move(const double * const vecIncremental)
1216 {
1217     TrackingMesh<NUM_NODES>::move(vecIncremental);
1218 }
1219 
1220 /* ----------------------------------------------------------------------
1221    scale mesh
1222 ------------------------------------------------------------------------- */
1223 
1224 template<int NUM_NODES, int NUM_NEIGH_MAX>
scale(double factor)1225 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::scale(double factor)
1226 {
1227     TrackingMesh<NUM_NODES>::scale(factor);
1228 
1229 }
1230 
1231 /* ----------------------------------------------------------------------
1232    rotate mesh
1233 ------------------------------------------------------------------------- */
1234 
1235 template<int NUM_NODES, int NUM_NEIGH_MAX>
rotate(const double * const totalQ,const double * const dQ,const double * const origin)1236 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::rotate(const double * const totalQ, const double * const dQ, const double * const origin)
1237 {
1238     TrackingMesh<NUM_NODES>::rotate(totalQ,dQ,origin);
1239 
1240     // find out if rotating every property is cheaper than
1241     // re-calculating them from the new nodes
1242 
1243 }
1244 
1245 template<int NUM_NODES, int NUM_NEIGH_MAX>
rotate(const double * const dQ,const double * const origin)1246 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::rotate(const double * const dQ, const double * const origin)
1247 {
1248     TrackingMesh<NUM_NODES>::rotate(dQ,origin);
1249 
1250     // find out if rotating every property is cheaper than
1251     // re-calculating them from the new nodes
1252 
1253 }
1254 
1255 /* ----------------------------------------------------------------------
1256    check if faces is planar
1257    used to check if a face can be used for particle insertion
1258 ------------------------------------------------------------------------- */
1259 
1260 template<int NUM_NODES, int NUM_NEIGH_MAX>
isPlanar()1261 bool SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::isPlanar()
1262 {
1263     int id_j;
1264     int flag = 0;
1265 
1266     int nlocal = this->sizeLocal();
1267 
1268     for(int i = 0; i < nlocal; i++)
1269     {
1270         if(flag) break;
1271 
1272         for(int ineigh = 0; ineigh < nNeighs_(i); ineigh++)
1273         {
1274             id_j = neighFaces_(i)[ineigh];
1275             if(!areCoplanarNeighs(TrackingMesh<NUM_NODES>::id(i),id_j))
1276                 flag = 1;
1277         }
1278     }
1279 
1280     MPI_Max_Scalar(flag,this->world);
1281 
1282     if(flag) return false;
1283     return true;
1284 }
1285 
1286 /* ----------------------------------------------------------------------
1287    check if point on surface - only valid if pos is in my subbox
1288 ------------------------------------------------------------------------- */
1289 
1290 template<int NUM_NODES, int NUM_NEIGH_MAX>
isOnSurface(double * pos)1291 bool SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::isOnSurface(double *pos)
1292 {
1293     bool on_surf = false;
1294 
1295     int nall = this->sizeLocal()+this->sizeGhost();
1296 
1297     // brute force
1298     // loop over ghosts as well as they might overlap my subbox
1299     for(int i = 0; i < nall; i++)
1300     {
1301         on_surf = on_surf || isInElement(pos,i);
1302         if(on_surf) break;
1303     }
1304 
1305     return on_surf;
1306 }
1307 
1308 /* ----------------------------------------------------------------------
1309    return number of active edges and corners for debugging
1310 ------------------------------------------------------------------------- */
1311 
1312 template<int NUM_NODES, int NUM_NEIGH_MAX>
n_active_edges(int i)1313 int SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::n_active_edges(int i)
1314 {
1315     int n = 0;
1316     if(i > this->size()) return n;
1317 
1318     if(edgeActive(i)[0]) n++;
1319     if(edgeActive(i)[1]) n++;
1320     if(edgeActive(i)[2]) n++;
1321     return n;
1322 }
1323 
1324 template<int NUM_NODES, int NUM_NEIGH_MAX>
n_active_corners(int i)1325 int SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::n_active_corners(int i)
1326 {
1327     int n = 0;
1328     if(i > this->size()) return n;
1329 
1330     if(cornerActive(i)[0]) n++;
1331     if(cornerActive(i)[1]) n++;
1332     if(cornerActive(i)[2]) n++;
1333     return n;
1334 }
1335 
1336 /* ----------------------------------------------------------------------
1337    edge-edge, edge-node, edge-point distance
1338 ------------------------------------------------------------------------- */
1339 
1340 template<int NUM_NODES, int NUM_NEIGH_MAX>
edgeEdgeDist(int iSrf,int iEdge,int jSrf,int jEdge)1341 double SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::edgeEdgeDist(int iSrf, int iEdge, int jSrf, int jEdge)
1342 {
1343     double d1,d2,d3,d4;
1344     d1 = edgeNodeDist(iSrf,iEdge,jSrf,jEdge);
1345     d2 = edgeNodeDist(iSrf,iEdge,jSrf,(jEdge+1)%NUM_NODES);
1346     d3 = edgeNodeDist(jSrf,jEdge,iSrf,(iEdge+1)%NUM_NODES);
1347     d4 = edgeNodeDist(jSrf,jEdge,iSrf,(iEdge+1)%NUM_NODES);
1348     return MathExtraLiggghts::min(d1,d2,d3,d4);
1349 }
1350 
1351 template<int NUM_NODES, int NUM_NEIGH_MAX>
edgeNodeDist(int iSrf,int iEdge,int jSrf,int jNode)1352 double SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::edgeNodeDist(int iSrf, int iEdge, int jSrf, int jNode)
1353 {
1354     return edgePointDist(iSrf, iEdge, MultiNodeMesh<NUM_NODES>::node_(jSrf)[jNode]);
1355 }
1356 
1357 template<int NUM_NODES, int NUM_NEIGH_MAX>
edgePointDist(int iSrf,int iEdge,double * point)1358 double SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::edgePointDist(int iSrf, int iEdge, double *point)
1359 {
1360         double nodeToP[3], dot;
1361 
1362         vectorSubtract3D(point,MultiNodeMesh<NUM_NODES>::node_(iSrf)[iEdge],nodeToP);
1363         dot = vectorDot3D(edgeVec(iSrf)[iEdge],nodeToP);
1364 
1365         if(dot < 0)
1366             return vectorMag3D(nodeToP);
1367 
1368         else if(dot > edgeLen(iSrf)[iEdge])
1369         {
1370             vectorSubtract3D(point,MultiNodeMesh<NUM_NODES>::node_(iSrf)[(iEdge+1)%NUM_NODES],nodeToP);
1371             return vectorMag3D(nodeToP);
1372         }
1373 
1374         else
1375             return MathExtraLiggghts::abs(vectorDot3D(edgeNorm(iSrf)[iEdge],nodeToP));
1376 }
1377 
1378 /* ----------------------------------------------------------------------
1379     Extrude a planar mesh in direction of the normal by length
1380 ------------------------------------------------------------------------- */
1381 
1382 template<int NUM_NODES, int NUM_NEIGH_MAX>
extrudePlanarMesh(const double length,double * & extrusion_tri_nodes,int & extrusion_tri_count)1383 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::extrudePlanarMesh(const double length, double * &extrusion_tri_nodes, int &extrusion_tri_count)
1384 {
1385     if(!this->isPlanar())
1386         this->error->all(FLERR, "Cannot extrude non-planar mesh");
1387 
1388     const int nlocal = this->sizeLocal();
1389 
1390     if (nlocal == 0 && this->comm->nprocs == 1)
1391         return;
1392 
1393     double extrudeVec[3];
1394     vectorCopy3D(surfaceNorm(0), extrudeVec);
1395     vectorScalarMult3D(extrudeVec, -length);
1396 
1397     extrusion_tri_count = 0;
1398     for(int i = 0; i < nlocal; i++)
1399         extrusion_tri_count += n_active_edges(i);
1400     extrusion_tri_count *= 2;
1401 
1402     // number of triangles local
1403     int count_local = extrusion_tri_count;
1404     // offset in nodes array for local proc
1405     int offset = 0;
1406     // number of tris for each processor
1407     int *n_tris_proc = new int[this->comm->nprocs];
1408     int *offsets_proc= new int[this->comm->nprocs];
1409     if (this->comm->nprocs > 1)
1410     {
1411         // get number of triangles of each processor
1412         MPI_Allgather(&count_local, 1, MPI_INT, &n_tris_proc[0], 1, MPI_INT, this->world);
1413         // compute total number of triangles and offset
1414         extrusion_tri_count = 0;
1415         for (int i = 0; i < this->comm->nprocs; i++)
1416         {
1417             extrusion_tri_count += n_tris_proc[i];
1418             if (i < this->comm->me)
1419                 offset += n_tris_proc[i];
1420             // now this becomes the number of vector_components * nodes * triangles
1421             n_tris_proc[i] *= 3*3;
1422             if (i > 0)
1423                 offsets_proc[i] = offsets_proc[i-1] + n_tris_proc[i-1];
1424             else
1425                 offsets_proc[i] = 0;
1426         }
1427     }
1428 
1429     if (extrusion_tri_count == 0)
1430         return;
1431 
1432     offset *= 3;
1433 
1434     // number of new triangles = number of active edges * 2 (for rectangle)
1435     extrusion_tri_nodes = new double[extrusion_tri_count*3*3];
1436     double * outbuf = new double[count_local*3*3];
1437 
1438     // loop over all triangles
1439     int count = 0;
1440     for(int i = 0; i < nlocal; i++)
1441     {
1442         // check if which edges of it are active (i.e. are on the boundary)
1443         // those will be extended along the normal by two triangles
1444         for (int j = 0; j < NUM_NODES; j++)
1445         {
1446             if(edgeActive(i)[j])
1447                 extrudeEdge(i, j, extrudeVec, count, outbuf);
1448         }
1449     }
1450 
1451     if (this->comm->nprocs > 1)
1452         MPI_Allgatherv(outbuf, count_local*3*3, MPI_DOUBLE,
1453                        &extrusion_tri_nodes[0], n_tris_proc, offsets_proc,
1454                        MPI_DOUBLE, this->world);
1455     else
1456         vectorCopyN(outbuf, extrusion_tri_nodes, count_local*3*3);
1457     delete[] outbuf;
1458     delete[] n_tris_proc;
1459     delete[] offsets_proc;
1460 }
1461 
1462 template<int NUM_NODES, int NUM_NEIGH_MAX>
extrudeEdge(const int nElem,const int edge,const double * const extrudeVec,int & count,double * extrusion_tri_nodes)1463 void SurfaceMesh<NUM_NODES,NUM_NEIGH_MAX>::extrudeEdge(const int nElem, const int edge, const double * const extrudeVec, int &count, double * extrusion_tri_nodes)
1464 {
1465     vectorCopy3D(MultiNodeMesh<NUM_NODES>::node_(nElem)[(edge+1)%NUM_NODES],
1466                  &extrusion_tri_nodes[count*3]);
1467     count++;
1468     vectorCopy3D(MultiNodeMesh<NUM_NODES>::node_(nElem)[edge],
1469                  &extrusion_tri_nodes[count*3]);
1470     count++;
1471     vectorAdd3D(&extrusion_tri_nodes[(count-2)*3], extrudeVec,
1472                 &extrusion_tri_nodes[count*3]);
1473     count++;
1474     vectorCopy3D(&extrusion_tri_nodes[(count-2)*3],
1475                  &extrusion_tri_nodes[count*3]);
1476     count++;
1477     vectorAdd3D(&extrusion_tri_nodes[(count-1)*3], extrudeVec,
1478                 &extrusion_tri_nodes[count*3]);
1479     count++;
1480     vectorCopy3D(&extrusion_tri_nodes[(count-3)*3],
1481                  &extrusion_tri_nodes[count*3]);
1482     count++;
1483 }
1484 
1485 #endif
1486