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