conformalVoronoiMesh.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2012-2018 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "conformalVoronoiMesh.H"
27 #include "initialPointsMethod.H"
28 #include "relaxationModel.H"
29 #include "faceAreaWeightModel.H"
30 #include "meshSearch.H"
31 #include "vectorTools.H"
32 #include "IOmanip.H"
33 #include "indexedCellChecks.H"
34 #include "controlMeshRefinement.H"
35 #include "smoothAlignmentSolver.H"
36 #include "OBJstream.H"
37 #include "indexedVertexOps.H"
38 #include "DelaunayMeshTools.H"
39 
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 
42 namespace Foam
43 {
44  defineTypeNameAndDebug(conformalVoronoiMesh, 0);
45 
46  template<>
47  const char* NamedEnum
48  <
50  5
51  >::names[] =
52  {
53  "internal",
54  "surface",
55  "featureEdge",
56  "featurePoint",
57  "constrained"
58  };
59 }
60 
63 
64 
65 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
66 
67 void Foam::conformalVoronoiMesh::cellSizeMeshOverlapsBackground() const
68 {
69  const cellShapeControlMesh& cellSizeMesh =
70  cellShapeControl_.shapeControlMesh();
71 
72  DynamicList<Foam::point> pts(number_of_vertices());
73 
74  for
75  (
76  Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
77  vit != finite_vertices_end();
78  ++vit
79  )
80  {
81  if (vit->internalOrBoundaryPoint() && !vit->referred())
82  {
83  pts.append(topoint(vit->point()));
84  }
85  }
86 
87  boundBox bb(pts);
88 
89  boundBox cellSizeMeshBb = cellSizeMesh.bounds();
90 
91  bool fullyContained = true;
92 
93  if (!cellSizeMeshBb.contains(bb))
94  {
95  Pout<< "Triangulation not fully contained in cell size mesh."
96  << endl;
97 
98  Pout<< "Cell Size Mesh Bounds = " << cellSizeMesh.bounds() << endl;
99  Pout<< "foamyHexMesh Bounds = " << bb << endl;
100 
101  fullyContained = false;
102  }
103 
104  reduce(fullyContained, andOp<unsigned int>());
105 
106  Info<< "Triangulation is "
107  << (fullyContained ? "fully" : "not fully")
108  << " contained in the cell size mesh"
109  << endl;
110 }
111 
112 
113 void Foam::conformalVoronoiMesh::insertInternalPoints
114 (
115  List<Point>& points,
116  bool distribute
117 )
118 {
119  label nPoints = points.size();
120 
121  if (Pstream::parRun())
122  {
123  reduce(nPoints, sumOp<label>());
124  }
125 
126  Info<< " " << nPoints << " points to insert..." << endl;
127 
128  if (Pstream::parRun() && distribute)
129  {
130  List<Foam::point> transferPoints(points.size());
131 
132  forAll(points, pI)
133  {
134  transferPoints[pI] = topoint(points[pI]);
135  }
136 
137  // Send the points that are not on this processor to the appropriate
138  // place
140  (
141  decomposition_().distributePoints(transferPoints)
142  );
143 
144  transferPoints.clear();
145 
146  map().distribute(points);
147  }
148 
149  label nVert = number_of_vertices();
150 
151  insert(points.begin(), points.end());
152 
153  label nInserted(number_of_vertices() - nVert);
154 
155  if (Pstream::parRun())
156  {
157  reduce(nInserted, sumOp<label>());
158  }
159 
160  Info<< " " << nInserted << " points inserted"
161  << ", failed to insert " << nPoints - nInserted
162  << " ("
163  << 100.0*(nPoints - nInserted)/(nInserted + small)
164  << " %)"<< endl;
165 
166  for
167  (
168  Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
169  vit != finite_vertices_end();
170  ++vit
171  )
172  {
174  {
175  vit->index() = getNewVertexIndex();
176  vit->type() = Vb::vtInternal;
177  }
178  }
179 }
180 
181 
182 Foam::Map<Foam::label> Foam::conformalVoronoiMesh::insertPointPairs
183 (
184  List<Vb>& vertices,
185  bool distribute,
186  bool reIndex
187 )
188 {
189  if (Pstream::parRun() && distribute)
190  {
191  autoPtr<mapDistribute> mapDist =
192  decomposition_().distributePoints(vertices);
193 
194  // Re-index the point pairs if one or both have been distributed.
195  // If both, remove
196 
197  // If added a point, then need to know its point pair
198  // If one moved, then update procIndex locally
199 
200  forAll(vertices, vI)
201  {
202  vertices[vI].procIndex() = Pstream::myProcNo();
203  }
204  }
205 
206  label preReinsertionSize(number_of_vertices());
207 
208  Map<label> oldToNewIndices =
209  this->DelaunayMesh<Delaunay>::insertPoints(vertices, reIndex);
210 
211  const label nReinserted = returnReduce
212  (
213  label(number_of_vertices()) - preReinsertionSize,
214  sumOp<label>()
215  );
216 
217  Info<< " Reinserted " << nReinserted << " vertices out of "
218  << returnReduce(vertices.size(), sumOp<label>())
219  << endl;
220 
221  return oldToNewIndices;
222 }
223 
224 
225 void Foam::conformalVoronoiMesh::insertSurfacePointPairs
226 (
227  const pointIndexHitAndFeatureList& surfaceHits,
228  const fileName fName,
229  DynamicList<Vb>& pts
230 )
231 {
232  forAll(surfaceHits, i)
233  {
234  vectorField norm(1);
235 
236  const pointIndexHit surfaceHit = surfaceHits[i].first();
237  const label featureIndex = surfaceHits[i].second();
238 
239  allGeometry_[featureIndex].getNormal
240  (
241  List<pointIndexHit>(1, surfaceHit),
242  norm
243  );
244 
245  const vector& normal = norm[0];
246 
247  const Foam::point& surfacePt(surfaceHit.hitPoint());
248 
250  geometryToConformTo_.meshableSide(featureIndex, surfaceHit);
251 
252  if (meshableSide == extendedFeatureEdgeMesh::BOTH)
253  {
254  createBafflePointPair
255  (
256  pointPairDistance(surfacePt),
257  surfacePt,
258  normal,
259  true,
260  pts
261  );
262  }
263  else if (meshableSide == extendedFeatureEdgeMesh::INSIDE)
264  {
265  createPointPair
266  (
267  pointPairDistance(surfacePt),
268  surfacePt,
269  normal,
270  true,
271  pts
272  );
273  }
274  else if (meshableSide == extendedFeatureEdgeMesh::OUTSIDE)
275  {
276  createPointPair
277  (
278  pointPairDistance(surfacePt),
279  surfacePt,
280  -normal,
281  true,
282  pts
283  );
284  }
285  else
286  {
288  << meshableSide << ", bad"
289  << endl;
290  }
291  }
292 
293  if (foamyHexMeshControls().objOutput() && fName != fileName::null)
294  {
295  DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
296  }
297 }
298 
299 
300 void Foam::conformalVoronoiMesh::insertEdgePointGroups
301 (
302  const pointIndexHitAndFeatureList& edgeHits,
303  const fileName fName,
304  DynamicList<Vb>& pts
305 )
306 {
307  forAll(edgeHits, i)
308  {
309  if (edgeHits[i].first().hit())
310  {
311  const extendedFeatureEdgeMesh& feMesh
312  (
313  geometryToConformTo_.features()[edgeHits[i].second()]
314  );
315 
316 // const bool isBaffle =
317 // geometryToConformTo_.isBaffleFeature(edgeHits[i].second());
318 
319  createEdgePointGroup
320  (
321  feMesh,
322  edgeHits[i].first(),
323  pts
324  );
325  }
326  }
327 
328  if (foamyHexMeshControls().objOutput() && fName != fileName::null)
329  {
330  DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
331  }
332 }
333 
334 
335 bool Foam::conformalVoronoiMesh::nearFeaturePt(const Foam::point& pt) const
336 {
337  scalar exclusionRangeSqr = featurePointExclusionDistanceSqr(pt);
338 
339  pointIndexHit info;
340  label featureHit;
341 
342  geometryToConformTo_.findFeaturePointNearest
343  (
344  pt,
345  exclusionRangeSqr,
346  info,
347  featureHit
348  );
349 
350  return info.hit();
351 }
352 
353 
354 bool Foam::conformalVoronoiMesh::surfacePtNearFeatureEdge
355 (
356  const Foam::point& pt
357 ) const
358 {
359  scalar exclusionRangeSqr = surfacePtExclusionDistanceSqr(pt);
360 
361  pointIndexHit info;
362  label featureHit;
363 
364  geometryToConformTo_.findEdgeNearest
365  (
366  pt,
367  exclusionRangeSqr,
368  info,
369  featureHit
370  );
371 
372  return info.hit();
373 }
374 
375 
376 void Foam::conformalVoronoiMesh::insertInitialPoints()
377 {
378  Info<< nl << "Inserting initial points" << endl;
379 
380  timeCheck("Before initial points call");
381 
382  List<Point> initPts = initialPointsMethod_->initialPoints();
383 
384  timeCheck("After initial points call");
385 
386  // Assume that the initial points method made the correct decision for
387  // which processor each point should be on, so give distribute = false
388  insertInternalPoints(initPts);
389 
390  if (initialPointsMethod_->fixInitialPoints())
391  {
392  for
393  (
394  Finite_vertices_iterator vit = finite_vertices_begin();
395  vit != finite_vertices_end();
396  ++vit
397  )
398  {
399  vit->fixed() = true;
400  }
401  }
402 
403  if (foamyHexMeshControls().objOutput())
404  {
406  (
407  time().path()/"initialPoints.obj",
408  *this,
410  );
411  }
412 }
413 
414 
415 void Foam::conformalVoronoiMesh::distribute()
416 {
417  if (!Pstream::parRun())
418  {
419  return;
420  }
421 
422  DynamicList<Foam::point> points(number_of_vertices());
423  DynamicList<Foam::indexedVertexEnum::vertexType> types
424  (
425  number_of_vertices()
426  );
427  DynamicList<scalar> sizes(number_of_vertices());
428  DynamicList<tensor> alignments(number_of_vertices());
429 
430  for
431  (
432  Finite_vertices_iterator vit = finite_vertices_begin();
433  vit != finite_vertices_end();
434  ++vit
435  )
436  {
437  if (vit->real())
438  {
439  points.append(topoint(vit->point()));
440  types.append(vit->type());
441  sizes.append(vit->targetCellSize());
442  alignments.append(vit->alignment());
443  }
444  }
445 
446  autoPtr<mapDistribute> mapDist =
447  DistributedDelaunayMesh<Delaunay>::distribute(decomposition_(), points);
448 
449  mapDist().distribute(types);
450  mapDist().distribute(sizes);
451  mapDist().distribute(alignments);
452 
453  // Reset the entire tessellation
455 
456  Info<< nl << " Inserting distributed tessellation" << endl;
457 
458  // Internal points have to be inserted first
459 
460  DynamicList<Vb> verticesToInsert(points.size());
461 
462  forAll(points, pI)
463  {
464  verticesToInsert.append
465  (
466  Vb
467  (
468  toPoint(points[pI]),
469  -1,
470  types[pI],
472  )
473  );
474 
475  verticesToInsert.last().targetCellSize() = sizes[pI];
476  verticesToInsert.last().alignment() = alignments[pI];
477  }
478 
479  this->rangeInsertWithInfo
480  (
481  verticesToInsert.begin(),
482  verticesToInsert.end(),
483  true
484  );
485 
486  Info<< " Total number of vertices after redistribution "
487  << returnReduce
488  (
489  label(number_of_vertices()), sumOp<label>()
490  )
491  << endl;
492 }
493 
494 
495 void Foam::conformalVoronoiMesh::buildCellSizeAndAlignmentMesh()
496 {
497  controlMeshRefinement meshRefinement
498  (
499  cellShapeControl_
500  );
501 
502  smoothAlignmentSolver meshAlignmentSmoother
503  (
504  cellShapeControl_.shapeControlMesh()
505  );
506 
507  meshRefinement.initialMeshPopulation(decomposition_);
508 
509  cellShapeControlMesh& cellSizeMesh = cellShapeControl_.shapeControlMesh();
510 
511  if (Pstream::parRun())
512  {
513  if (!distributeBackground(cellSizeMesh))
514  {
515  // Synchronise the cell size mesh if it has not been distributed
516  cellSizeMesh.distribute(decomposition_);
517  }
518  }
519 
520  const dictionary& motionControlDict
521  = foamyHexMeshControls().foamyHexMeshDict().subDict("motionControl");
522 
523  label nMaxIter = readLabel
524  (
525  motionControlDict.lookup("maxRefinementIterations")
526  );
527 
528  Info<< "Maximum number of refinement iterations : " << nMaxIter << endl;
529 
530  for (label i = 0; i < nMaxIter; ++i)
531  {
532  label nAdded = meshRefinement.refineMesh(decomposition_);
533  // cellShapeControl_.refineMesh(decomposition_);
534  reduce(nAdded, sumOp<label>());
535 
536  if (Pstream::parRun())
537  {
538  cellSizeMesh.distribute(decomposition_);
539  }
540 
541  Info<< " Iteration " << i
542  << " Added = " << nAdded << " points"
543  << endl;
544 
545  if (nAdded == 0)
546  {
547  break;
548  }
549  }
550 
551  if (Pstream::parRun())
552  {
553  // Need to distribute the cell size mesh to cover the background mesh
554  if (!distributeBackground(cellSizeMesh))
555  {
556  cellSizeMesh.distribute(decomposition_);
557  }
558  }
559 
560  label maxSmoothingIterations = readLabel
561  (
562  motionControlDict.lookup("maxSmoothingIterations")
563  );
564  meshAlignmentSmoother.smoothAlignments(maxSmoothingIterations);
565 
566  Info<< "Background cell size and alignment mesh:" << endl;
567  cellSizeMesh.printInfo(Info);
568 
569  Info<< "Triangulation is "
570  << (cellSizeMesh.is_valid() ? "valid" : "not valid!" )
571  << endl;
572 
573  if (foamyHexMeshControls().writeCellShapeControlMesh())
574  {
575  // cellSizeMesh.writeTriangulation();
576  cellSizeMesh.write();
577  }
578 
579  if (foamyHexMeshControls().printVertexInfo())
580  {
581  cellSizeMesh.printVertexInfo(Info);
582  }
583 
584 // Info<< "Estimated number of cells in final mesh = "
585 // << returnReduce
586 // (
587 // cellSizeMesh.estimateCellCount(decomposition_),
588 // sumOp<label>()
589 // )
590 // << endl;
591 }
592 
593 
594 void Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()
595 {
596  Info<< nl << "Calculating target cell alignment and size" << endl;
597 
598  for
599  (
600  Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
601  vit != finite_vertices_end();
602  vit++
603  )
604  {
605  if (vit->internalOrBoundaryPoint())
606  {
607  pointFromPoint pt = topoint(vit->point());
608 
609  cellShapeControls().cellSizeAndAlignment
610  (
611  pt,
612  vit->targetCellSize(),
613  vit->alignment()
614  );
615  }
616  }
617 }
618 
619 
620 Foam::face Foam::conformalVoronoiMesh::buildDualFace
621 (
622  const Delaunay::Finite_edges_iterator& eit
623 ) const
624 {
625  Cell_circulator ccStart = incident_cells(*eit);
626  Cell_circulator cc1 = ccStart;
627  Cell_circulator cc2 = cc1;
628 
629  // Advance the second circulator so that it always stays on the next
630  // cell around the edge;
631  cc2++;
632 
633  DynamicList<label> verticesOnFace;
634 
635  label nUniqueVertices = 0;
636 
637  do
638  {
639  if
640  (
641  cc1->hasFarPoint() || cc2->hasFarPoint()
642  || is_infinite(cc1) || is_infinite(cc2)
643  )
644  {
645  Cell_handle c = eit->first;
646  Vertex_handle vA = c->vertex(eit->second);
647  Vertex_handle vB = c->vertex(eit->third);
648 
649 // DelaunayMeshTools::drawDelaunayCell(Pout, cc1);
650 // DelaunayMeshTools::drawDelaunayCell(Pout, cc2);
651 
653  << "Dual face uses circumcenter defined by a "
654  << "Delaunay tetrahedron with no internal "
655  << "or boundary points. Defining Delaunay edge ends: "
656  << vA->info() << " "
657  << vB->info() << nl
658  <<endl;//<< exit(FatalError);
659  }
660  else
661  {
662  label cc1I = cc1->cellIndex();
663  label cc2I = cc2->cellIndex();
664 
665  if (cc1I != cc2I)
666  {
667  if (findIndex(verticesOnFace, cc1I) == -1)
668  {
669  nUniqueVertices++;
670  }
671 
672  verticesOnFace.append(cc1I);
673  }
674  }
675 
676  cc1++;
677  cc2++;
678 
679  } while (cc1 != ccStart);
680 
681  verticesOnFace.shrink();
682 
683  if (verticesOnFace.size() >= 3 && nUniqueVertices < 3)
684  {
685  // There are not enough unique vertices on this face to
686  // justify its size, it may have a form like:
687 
688  // Vertices:
689  // A B
690  // A B
691 
692  // Face:
693  // ABAB
694 
695  // Setting the size to be below 3, so that it will not be
696  // created
697 
698  verticesOnFace.setSize(nUniqueVertices);
699  }
700 
701  return face(verticesOnFace);
702 }
703 
704 
705 Foam::label Foam::conformalVoronoiMesh::maxFilterCount
706 (
707  const Delaunay::Finite_edges_iterator& eit
708 ) const
709 {
710  Cell_circulator ccStart = incident_cells(*eit);
711  Cell_circulator cc = ccStart;
712 
713  label maxFC = 0;
714 
715  do
716  {
717  if (cc->hasFarPoint())
718  {
719  Cell_handle c = eit->first;
720  Vertex_handle vA = c->vertex(eit->second);
721  Vertex_handle vB = c->vertex(eit->third);
722 
724  << "Dual face uses circumcenter defined by a "
725  << "Delaunay tetrahedron with no internal "
726  << "or boundary points. Defining Delaunay edge ends: "
727  << topoint(vA->point()) << " "
728  << topoint(vB->point()) << nl
729  << exit(FatalError);
730  }
731 
732  if (cc->filterCount() > maxFC)
733  {
734  maxFC = cc->filterCount();
735  }
736 
737  cc++;
738 
739  } while (cc != ccStart);
740 
741  return maxFC;
742 }
743 
744 
745 bool Foam::conformalVoronoiMesh::ownerAndNeighbour
746 (
747  Vertex_handle vA,
748  Vertex_handle vB,
749  label& owner,
750  label& neighbour
751 ) const
752 {
753  bool reverse = false;
754 
755  owner = -1;
756  neighbour = -1;
757 
758  label dualCellIndexA = vA->index();
759 
760  if (!vA->internalOrBoundaryPoint() || vA->referred())
761  {
762  if (!vA->constrained())
763  {
764  dualCellIndexA = -1;
765  }
766  }
767 
768  label dualCellIndexB = vB->index();
769 
770  if (!vB->internalOrBoundaryPoint() || vB->referred())
771  {
772  if (!vB->constrained())
773  {
774  dualCellIndexB = -1;
775  }
776  }
777 
778  if (dualCellIndexA == -1 && dualCellIndexB == -1)
779  {
781  << "Attempting to create a face joining "
782  << "two unindexed dual cells "
783  << exit(FatalError);
784  }
785  else if (dualCellIndexA == -1 || dualCellIndexB == -1)
786  {
787  // boundary face, find which is the owner
788 
789  if (dualCellIndexA == -1)
790  {
791  owner = dualCellIndexB;
792 
793  reverse = true;
794  }
795  else
796  {
797  owner = dualCellIndexA;
798  }
799  }
800  else
801  {
802  // internal face, find the lower cell to be the owner
803 
804  if (dualCellIndexB > dualCellIndexA)
805  {
806  owner = dualCellIndexA;
807  neighbour = dualCellIndexB;
808  }
809  else
810  {
811  owner = dualCellIndexB;
812  neighbour = dualCellIndexA;
813 
814  // reverse face order to correctly orientate normal
815  reverse = true;
816  }
817  }
818 
819  return reverse;
820 }
821 
822 
823 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
824 
825 Foam::conformalVoronoiMesh::conformalVoronoiMesh
826 (
827  const Time& runTime,
828  const dictionary& foamyHexMeshDict
829 )
830 :
831  DistributedDelaunayMesh<Delaunay>(runTime),
832  runTime_(runTime),
833  rndGen_(64293*Pstream::myProcNo()),
834  foamyHexMeshControls_(foamyHexMeshDict),
835  allGeometry_
836  (
837  IOobject
838  (
839  "cvSearchableSurfaces",
840  runTime_.constant(),
841  "triSurface",
842  runTime_,
843  IOobject::MUST_READ,
844  IOobject::NO_WRITE
845  ),
846  foamyHexMeshDict.subDict("geometry"),
847  foamyHexMeshDict.lookupOrDefault("singleRegionName", true)
848  ),
849  geometryToConformTo_
850  (
851  runTime_,
852  rndGen_,
853  allGeometry_,
854  foamyHexMeshDict.subDict("surfaceConformation")
855  ),
856  decomposition_
857  (
858  Pstream::parRun()
859  ? new backgroundMeshDecomposition
860  (
861  runTime_,
862  rndGen_,
863  geometryToConformTo_,
864  foamyHexMeshControls().foamyHexMeshDict().subDict
865  (
866  "backgroundMeshDecomposition"
867  )
868  )
869  : nullptr
870  ),
871  cellShapeControl_
872  (
873  runTime_,
874  foamyHexMeshControls_,
875  allGeometry_,
876  geometryToConformTo_
877  ),
878  limitBounds_(),
879  ptPairs_(*this),
880  ftPtConformer_(*this),
881  edgeLocationTreePtr_(),
882  surfacePtLocationTreePtr_(),
883  surfaceConformationVertices_(),
884  initialPointsMethod_
885  (
886  initialPointsMethod::New
887  (
888  foamyHexMeshDict.subDict("initialPoints"),
889  runTime_,
890  rndGen_,
891  geometryToConformTo_,
892  cellShapeControl_,
893  decomposition_
894  )
895  ),
896  relaxationModel_
897  (
898  relaxationModel::New
899  (
900  foamyHexMeshDict.subDict("motionControl"),
901  runTime_
902  )
903  ),
904  faceAreaWeightModel_
905  (
906  faceAreaWeightModel::New
907  (
908  foamyHexMeshDict.subDict("motionControl")
909  )
910  )
911 {}
912 
913 
914 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
915 
917 {}
918 
919 
920 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
921 
923 {
924  if (foamyHexMeshControls().objOutput())
925  {
926  geometryToConformTo_.writeFeatureObj("foamyHexMesh");
927  }
928 
929  buildCellSizeAndAlignmentMesh();
930 
931  insertInitialPoints();
932 
933  insertFeaturePoints(true);
934 
935  setVertexSizeAndAlignment();
936 
937  cellSizeMeshOverlapsBackground();
938 
939  // Improve the guess that the backgroundMeshDecomposition makes with the
940  // initial positions. Use before building the surface conformation to
941  // better balance the surface conformation load.
942  distributeBackground(*this);
943 
944  buildSurfaceConformation();
945 
946  // The introduction of the surface conformation may have distorted the
947  // balance of vertices, distribute if necessary.
948  distributeBackground(*this);
949 
950  if (Pstream::parRun())
951  {
952  sync(decomposition_().procBounds());
953  }
954 
955  // Do not store the surface conformation until after it has been
956  // (potentially) redistributed.
957  storeSurfaceConformation();
958 
959  // Report any Delaunay vertices that do not think that they are in the
960  // domain the processor they are on.
961  // reportProcessorOccupancy();
962 
963  cellSizeMeshOverlapsBackground();
964 
966  {
968  }
969 
970  if (foamyHexMeshControls().objOutput())
971  {
973  (
974  time().path()/"internalPoints_" + time().timeName() + ".obj",
975  *this,
978  );
979  }
980 }
981 
982 
984 {
985  if (Pstream::parRun())
986  {
987  decomposition_.reset
988  (
989  new backgroundMeshDecomposition
990  (
991  runTime_,
992  rndGen_,
993  geometryToConformTo_,
994  foamyHexMeshControls().foamyHexMeshDict().subDict
995  (
996  "backgroundMeshDecomposition"
997  )
998  )
999  );
1000  }
1001 
1002  insertInitialPoints();
1003 
1004  insertFeaturePoints();
1005 
1006  // Improve the guess that the backgroundMeshDecomposition makes with the
1007  // initial positions. Use before building the surface conformation to
1008  // better balance the surface conformation load.
1009  distributeBackground(*this);
1010 
1011  buildSurfaceConformation();
1012 
1013  // The introduction of the surface conformation may have distorted the
1014  // balance of vertices, distribute if necessary.
1015  distributeBackground(*this);
1016 
1017  if (Pstream::parRun())
1018  {
1019  sync(decomposition_().procBounds());
1020  }
1021 
1022  cellSizeMeshOverlapsBackground();
1023 
1025  {
1027  }
1028 }
1029 
1030 
1032 {
1033  timeCheck("Start of move");
1034 
1035  scalar relaxation = relaxationModel_->relaxation();
1036 
1037  Info<< nl << "Relaxation = " << relaxation << endl;
1038 
1039  pointField dualVertices(number_of_finite_cells());
1040 
1041  this->resetCellCount();
1042 
1043  // Find the dual point of each tetrahedron and assign it an index.
1044  for
1045  (
1046  Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1047  cit != finite_cells_end();
1048  ++cit
1049  )
1050  {
1051  cit->cellIndex() = Cb::ctUnassigned;
1052 
1053  if (cit->anyInternalOrBoundaryDualVertex())
1054  {
1055  cit->cellIndex() = getNewCellIndex();
1056 
1057  dualVertices[cit->cellIndex()] = cit->dual();
1058  }
1059 
1060  if (cit->hasFarPoint())
1061  {
1062  cit->cellIndex() = Cb::ctFar;
1063  }
1064  }
1065 
1066  dualVertices.setSize(cellCount());
1067 
1068  setVertexSizeAndAlignment();
1069 
1070  timeCheck("Determined sizes and alignments");
1071 
1072  Info<< nl << "Determining vertex displacements" << endl;
1073 
1074  vectorField cartesianDirections(3);
1075 
1076  cartesianDirections[0] = vector(1, 0, 0);
1077  cartesianDirections[1] = vector(0, 1, 0);
1078  cartesianDirections[2] = vector(0, 0, 1);
1079 
1080  vectorField displacementAccumulator
1081  (
1082  number_of_vertices(),
1083  Zero
1084  );
1085 
1086  PackedBoolList pointToBeRetained
1087  (
1088  number_of_vertices(),
1089  true
1090  );
1091 
1092  DynamicList<Point> pointsToInsert(number_of_vertices());
1093 
1094  for
1095  (
1096  Delaunay::Finite_edges_iterator eit = finite_edges_begin();
1097  eit != finite_edges_end();
1098  ++eit
1099  )
1100  {
1101  Cell_handle c = eit->first;
1102  Vertex_handle vA = c->vertex(eit->second);
1103  Vertex_handle vB = c->vertex(eit->third);
1104 
1105  if
1106  (
1107  (
1108  vA->internalPoint() && !vA->referred()
1109  && vB->internalOrBoundaryPoint()
1110  )
1111  || (
1112  vB->internalPoint() && !vB->referred()
1113  && vA->internalOrBoundaryPoint()
1114  )
1115  )
1116  {
1117  pointFromPoint dVA = topoint(vA->point());
1118  pointFromPoint dVB = topoint(vB->point());
1119 
1120  Field<vector> alignmentDirsA
1121  (
1122  vA->alignment().T() & cartesianDirections
1123  );
1124  Field<vector> alignmentDirsB
1125  (
1126  vB->alignment().T() & cartesianDirections
1127  );
1128 
1129  Field<vector> alignmentDirs(alignmentDirsA);
1130 
1131  forAll(alignmentDirsA, aA)
1132  {
1133  const vector& a = alignmentDirsA[aA];
1134 
1135  scalar maxDotProduct = 0.0;
1136 
1137  forAll(alignmentDirsB, aB)
1138  {
1139  const vector& b = alignmentDirsB[aB];
1140 
1141  const scalar dotProduct = a & b;
1142 
1143  if (mag(dotProduct) > maxDotProduct)
1144  {
1145  maxDotProduct = mag(dotProduct);
1146 
1147  alignmentDirs[aA] = a + sign(dotProduct)*b;
1148 
1149  alignmentDirs[aA] /= mag(alignmentDirs[aA]);
1150  }
1151  }
1152  }
1153 
1154  vector rAB = dVA - dVB;
1155 
1156  scalar rABMag = mag(rAB);
1157 
1158  if (rABMag < small)
1159  {
1160  // Removal of close points
1161 
1162  if
1163  (
1164  vA->internalPoint() && !vA->referred() && !vA->fixed()
1165  && vB->internalPoint() && !vB->referred() && !vB->fixed()
1166  )
1167  {
1168  // Only insert a point at the midpoint of
1169  // the short edge if neither attached
1170  // point has already been identified to be
1171  // removed.
1172 
1173  if
1174  (
1175  pointToBeRetained[vA->index()] == true
1176  && pointToBeRetained[vB->index()] == true
1177  )
1178  {
1179  const Foam::point pt(0.5*(dVA + dVB));
1180 
1181  if (internalPointIsInside(pt))
1182  {
1183  pointsToInsert.append(toPoint(pt));
1184  }
1185  }
1186  }
1187 
1188  if (vA->internalPoint() && !vA->referred() && !vA->fixed())
1189  {
1190  pointToBeRetained[vA->index()] = false;
1191  }
1192 
1193  if (vB->internalPoint() && !vB->referred() && !vB->fixed())
1194  {
1195  pointToBeRetained[vB->index()] = false;
1196  }
1197 
1198  // Do not consider this Delaunay edge any further
1199 
1200  continue;
1201  }
1202 
1203  forAll(alignmentDirs, aD)
1204  {
1205  vector& alignmentDir = alignmentDirs[aD];
1206 
1207  scalar dotProd = rAB & alignmentDir;
1208 
1209  if (dotProd < 0)
1210  {
1211  // swap the direction of the alignment so that has the
1212  // same sense as rAB
1213  alignmentDir *= -1;
1214  dotProd *= -1;
1215  }
1216 
1217  const scalar alignmentDotProd = dotProd/rABMag;
1218 
1219  if
1220  (
1221  alignmentDotProd
1222  > foamyHexMeshControls().cosAlignmentAcceptanceAngle()
1223  )
1224  {
1225  scalar targetCellSize =
1227 
1228  scalar targetFaceArea = sqr(targetCellSize);
1229 
1230  const vector originalAlignmentDir = alignmentDir;
1231 
1232  // Update cell size and face area
1234  (
1235  alignmentDir,
1236  targetFaceArea,
1237  targetCellSize
1238  );
1239 
1240  // Vector to move end points around middle of vector
1241  // to align edge (i.e. dual face normal) with alignment
1242  // directions.
1243  vector delta = alignmentDir - 0.5*rAB;
1244 
1245  face dualFace = buildDualFace(eit);
1246 
1247 // Pout<< dualFace << endl;
1248 // Pout<< " " << vA->info() << endl;
1249 // Pout<< " " << vB->info() << endl;
1250 
1251  const scalar faceArea = dualFace.mag(dualVertices);
1252 
1253  // Update delta vector
1255  (
1256  originalAlignmentDir,
1257  targetCellSize,
1258  rABMag,
1259  delta
1260  );
1261 
1262 // if (targetFaceArea == 0)
1263 // {
1264 // Pout<< vA->info() << vB->info();
1265 //
1266 // Cell_handle ch = locate(vA->point());
1267 // if (is_infinite(ch))
1268 // {
1269 // Pout<< "vA " << vA->targetCellSize() << endl;
1270 // }
1271 //
1272 // ch = locate(vB->point());
1273 // if (is_infinite(ch))
1274 // {
1275 // Pout<< "vB " << vB->targetCellSize() << endl;
1276 // }
1277 // }
1278 
1279  delta *= faceAreaWeightModel_->faceAreaWeight
1280  (
1281  faceArea/targetFaceArea
1282  );
1283 
1284  if
1285  (
1286  (
1287  (vA->internalPoint() && vB->internalPoint())
1288  && (!vA->referred() || !vB->referred())
1289 // ||
1290 // (
1291 // vA->referredInternalPoint()
1292 // && vB->referredInternalPoint()
1293 // )
1294  )
1295  && rABMag
1297  *targetCellSize
1298  && faceArea
1300  *targetFaceArea
1301  && alignmentDotProd
1303  )
1304  {
1305  // Point insertion
1306  if
1307  (
1308  !geometryToConformTo_.findSurfaceAnyIntersection
1309  (
1310  dVA,
1311  dVB
1312  )
1313  )
1314  {
1315  const Foam::point newPt(0.5*(dVA + dVB));
1316 
1317  // Prevent insertions spanning surfaces
1318  if (internalPointIsInside(newPt))
1319  {
1320  if (Pstream::parRun())
1321  {
1322  if
1323  (
1324  decomposition().positionOnThisProcessor
1325  (
1326  newPt
1327  )
1328  )
1329  {
1330  pointsToInsert.append(toPoint(newPt));
1331  }
1332  }
1333  else
1334  {
1335  pointsToInsert.append(toPoint(newPt));
1336  }
1337  }
1338  }
1339  }
1340  else if
1341  (
1342  (
1343  (vA->internalPoint() && !vA->referred())
1344  || (vB->internalPoint() && !vB->referred())
1345  )
1346  && rABMag
1348  *targetCellSize
1349  )
1350  {
1351  // Point removal
1352  if
1353  (
1354  (
1355  vA->internalPoint()
1356  && !vA->referred()
1357  && !vA->fixed()
1358  )
1359  &&
1360  (
1361  vB->internalPoint()
1362  && !vB->referred()
1363  && !vB->fixed()
1364  )
1365  )
1366  {
1367  // Only insert a point at the midpoint of
1368  // the short edge if neither attached
1369  // point has already been identified to be
1370  // removed.
1371  if
1372  (
1373  pointToBeRetained[vA->index()] == true
1374  && pointToBeRetained[vB->index()] == true
1375  )
1376  {
1377  const Foam::point pt(0.5*(dVA + dVB));
1378 
1379  if (internalPointIsInside(pt))
1380  {
1381  pointsToInsert.append(toPoint(pt));
1382  }
1383  }
1384  }
1385 
1386  if
1387  (
1388  vA->internalPoint()
1389  && !vA->referred()
1390  && !vA->fixed()
1391  )
1392  {
1393  pointToBeRetained[vA->index()] = false;
1394  }
1395 
1396  if
1397  (
1398  vB->internalPoint()
1399  && !vB->referred()
1400  && !vB->fixed()
1401  )
1402  {
1403  pointToBeRetained[vB->index()] = false;
1404  }
1405  }
1406  else
1407  {
1408  if
1409  (
1410  vA->internalPoint()
1411  && !vA->referred()
1412  && !vA->fixed()
1413  )
1414  {
1415  if (vB->fixed())
1416  {
1417  displacementAccumulator[vA->index()] += 2*delta;
1418  }
1419  else
1420  {
1421  displacementAccumulator[vA->index()] += delta;
1422  }
1423  }
1424 
1425  if
1426  (
1427  vB->internalPoint()
1428  && !vB->referred()
1429  && !vB->fixed()
1430  )
1431  {
1432  if (vA->fixed())
1433  {
1434  displacementAccumulator[vB->index()] -= 2*delta;
1435  }
1436  else
1437  {
1438  displacementAccumulator[vB->index()] -= delta;
1439  }
1440  }
1441  }
1442  }
1443  }
1444  }
1445  }
1446 
1447  Info<< "Limit displacements" << endl;
1448 
1449  // Limit displacements that pierce, or get too close to the surface
1450  for
1451  (
1452  Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1453  vit != finite_vertices_end();
1454  ++vit
1455  )
1456  {
1457  if (vit->internalPoint() && !vit->referred() && !vit->fixed())
1458  {
1459  if (pointToBeRetained[vit->index()] == true)
1460  {
1461  limitDisplacement
1462  (
1463  vit,
1464  displacementAccumulator[vit->index()]
1465  );
1466  }
1467  }
1468  }
1469 
1470  vector totalDisp = gSum(displacementAccumulator);
1471  scalar totalDist = gSum(mag(displacementAccumulator));
1472 
1473  displacementAccumulator *= relaxation;
1474 
1475  Info<< "Sum displacements" << endl;
1476 
1477  label nPointsToRetain = 0;
1478  label nPointsToRemove = 0;
1479 
1480  for
1481  (
1482  Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1483  vit != finite_vertices_end();
1484  ++vit
1485  )
1486  {
1487  if (vit->internalPoint() && !vit->referred() && !vit->fixed())
1488  {
1489  if (pointToBeRetained[vit->index()] == true)
1490  {
1491  // Convert vit->point() to FOAM vector (double) to do addition,
1492  // avoids memory increase because a record of the constructions
1493  // would be kept otherwise.
1494  // See cgal-discuss@lists-sop.inria.fr:
1495  // "Memory issue with openSUSE 11.3, exact kernel, adding
1496  // points/vectors"
1497  // 14/1/2011.
1498  // Only necessary if using an exact constructions kernel
1499  // (extended precision)
1500  Foam::point pt
1501  (
1502  topoint(vit->point())
1503  + displacementAccumulator[vit->index()]
1504  );
1505 
1506  if (internalPointIsInside(pt))
1507  {
1508  pointsToInsert.append(toPoint(pt));
1509  nPointsToRemove++;
1510  }
1511 
1512  nPointsToRetain++;
1513  }
1514  }
1515  }
1516 
1517  pointsToInsert.shrink();
1518 
1519  Info<< returnReduce
1520  (
1521  nPointsToRetain - nPointsToRemove,
1522  sumOp<label>()
1523  )
1524  << " internal points are outside the domain. "
1525  << "They will not be inserted." << endl;
1526 
1527  // Save displacements to file.
1528  if (foamyHexMeshControls().objOutput() && time().writeTime())
1529  {
1530  Info<< "Writing point displacement vectors to file." << endl;
1531  OFstream str
1532  (
1533  time().path()/"displacements_" + runTime_.timeName() + ".obj"
1534  );
1535 
1536  for
1537  (
1538  Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1539  vit != finite_vertices_end();
1540  ++vit
1541  )
1542  {
1543  if (vit->internalPoint() && !vit->referred())
1544  {
1545  if (pointToBeRetained[vit->index()] == true)
1546  {
1547  meshTools::writeOBJ(str, topoint(vit->point()));
1548 
1549  str << "vn "
1550  << displacementAccumulator[vit->index()][0] << " "
1551  << displacementAccumulator[vit->index()][1] << " "
1552  << displacementAccumulator[vit->index()][2] << " "
1553  << endl;
1554  }
1555  }
1556  }
1557  }
1558 
1559  // Remove the entire tessellation
1561 
1562  timeCheck("Displacement calculated");
1563 
1564  Info<< nl<< "Inserting displaced tessellation" << endl;
1565 
1566  insertFeaturePoints(true);
1567 
1568  insertInternalPoints(pointsToInsert, true);
1569 
1570  // Fix points that have not been significantly displaced
1571 // for
1572 // (
1573 // Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1574 // vit != finite_vertices_end();
1575 // ++vit
1576 // )
1577 // {
1578 // if (vit->internalPoint())
1579 // {
1580 // if
1581 // (
1582 // mag(displacementAccumulator[vit->index()])
1583 // < 0.1*targetCellSize(topoint(vit->point()))
1584 // )
1585 // {
1586 // vit->setVertexFixed();
1587 // }
1588 // }
1589 // }
1590 
1591  timeCheck("Internal points inserted");
1592 
1593  {
1594  // Check that no index is shared between any of the local points
1595  labelHashSet usedIndices;
1596  for
1597  (
1598  Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1599  vit != finite_vertices_end();
1600  ++vit
1601  )
1602  {
1603  if (!vit->referred() && !usedIndices.insert(vit->index()))
1604  {
1606  << "Index already used! Could not insert: " << nl
1607  << vit->info()
1608  << abort(FatalError);
1609  }
1610  }
1611  }
1612 
1613  conformToSurface();
1614 
1615  if (foamyHexMeshControls().objOutput())
1616  {
1618  (
1619  time().path()/"internalPoints_" + time().timeName() + ".obj",
1620  *this,
1622  );
1623 
1624  if (reconformToSurface())
1625  {
1627  (
1628  time().path()/"boundaryPoints_" + time().timeName() + ".obj",
1629  *this
1630  );
1631 
1633  (
1634  time().path()/"internalBoundaryPoints_" + time().timeName()
1635  + ".obj",
1636  *this,
1639  );
1640 
1642  (
1643  time().path()/"externalBoundaryPoints_" + time().timeName()
1644  + ".obj",
1645  *this,
1648  );
1649 
1650  OBJstream multipleIntersections
1651  (
1652  "multipleIntersections_"
1653  + time().timeName()
1654  + ".obj"
1655  );
1656 
1657  for
1658  (
1659  Delaunay::Finite_edges_iterator eit = finite_edges_begin();
1660  eit != finite_edges_end();
1661  ++eit
1662  )
1663  {
1664  Cell_handle c = eit->first;
1665  Vertex_handle vA = c->vertex(eit->second);
1666  Vertex_handle vB = c->vertex(eit->third);
1667 
1668  Foam::point ptA(topoint(vA->point()));
1669  Foam::point ptB(topoint(vB->point()));
1670 
1671  List<pointIndexHit> surfHits;
1672  labelList hitSurfaces;
1673 
1675  (
1676  ptA,
1677  ptB,
1678  surfHits,
1679  hitSurfaces
1680  );
1681 
1682  if
1683  (
1684  surfHits.size() >= 2
1685  || (
1686  surfHits.size() == 0
1687  && (
1688  (vA->externalBoundaryPoint()
1689  && vB->internalBoundaryPoint())
1690  || (vB->externalBoundaryPoint()
1691  && vA->internalBoundaryPoint())
1692  )
1693  )
1694  )
1695  {
1696  multipleIntersections.write(linePointRef(ptA, ptB));
1697  }
1698  }
1699  }
1700  }
1701 
1702  timeCheck("After conformToSurface");
1703 
1705  {
1707  }
1708 
1709  if (time().writeTime())
1710  {
1711  writeMesh(time().timeName());
1712  }
1713 
1714  Info<< nl
1715  << "Total displacement = " << totalDisp << nl
1716  << "Total distance = " << totalDist << nl
1717  << endl;
1718 }
1719 
1720 
1721 void Foam::conformalVoronoiMesh::checkCoPlanarCells() const
1722 {
1723  typedef CGAL::Exact_predicates_exact_constructions_kernel Kexact;
1724  typedef CGAL::Point_3<Kexact> PointExact;
1725 
1726  if (!is_valid())
1727  {
1728  Pout<< "Triangulation is invalid!" << endl;
1729  }
1730 
1731  OFstream str("badCells.obj");
1732 
1733  label badCells = 0;
1734 
1735  for
1736  (
1737  Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1738  cit != finite_cells_end();
1739  ++cit
1740  )
1741  {
1742  const scalar quality = foamyHexMeshChecks::coplanarTet(cit, 1e-16);
1743 
1744  if (quality == 0)
1745  {
1746  Pout<< "COPLANAR: " << cit->info() << nl
1747  << " quality = " << quality << nl
1748  << " dual = " << topoint(cit->dual()) << endl;
1749 
1750  DelaunayMeshTools::drawDelaunayCell(str, cit, badCells++);
1751 
1752  FixedList<PointExact, 4> cellVerticesExact(PointExact(0,0,0));
1753  forAll(cellVerticesExact, vI)
1754  {
1755  cellVerticesExact[vI] = PointExact
1756  (
1757  cit->vertex(vI)->point().x(),
1758  cit->vertex(vI)->point().y(),
1759  cit->vertex(vI)->point().z()
1760  );
1761  }
1762 
1763  PointExact synchronisedDual = CGAL::circumcenter<Kexact>
1764  (
1765  cellVerticesExact[0],
1766  cellVerticesExact[1],
1767  cellVerticesExact[2],
1768  cellVerticesExact[3]
1769  );
1770 
1771  Foam::point exactPt
1772  (
1773  CGAL::to_double(synchronisedDual.x()),
1774  CGAL::to_double(synchronisedDual.y()),
1775  CGAL::to_double(synchronisedDual.z())
1776  );
1777 
1778  Info<< "inexact = " << cit->dual() << nl
1779  << "exact = " << exactPt << endl;
1780  }
1781  }
1782 
1783  Pout<< "There are " << badCells << " bad cells out of "
1784  << number_of_finite_cells() << endl;
1785 
1786 
1787  label nNonGabriel = 0;
1788  for
1789  (
1790  Delaunay::Finite_facets_iterator fit = finite_facets_begin();
1791  fit != finite_facets_end();
1792  ++fit
1793  )
1794  {
1795  if (!is_Gabriel(*fit))
1796  {
1797  nNonGabriel++;//Pout<< "Non-gabriel face" << endl;
1798  }
1799  }
1800 
1801  Pout<< "There are " << nNonGabriel << " non-Gabriel faces out of "
1802  << number_of_finite_facets() << endl;
1803 }
1804 
1805 
1806 // ************************************************************************* //
scalar delta
dimensionedScalar sign(const dimensionedScalar &ds)
Delaunay::Cell_handle Cell_handle
CGAL::Delaunay_triangulation_3< K, Tds, CompactLocator > Delaunay
label cellCount() const
Return the cell count (the next unique cell index)
Definition: DelaunayMeshI.H:93
~conformalVoronoiMesh()
Destructor.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Delaunay::Vertex_handle Vertex_handle
pointFromPoint topoint(const Point &P)
void writeOBJ(const fileName &fName, const List< Foam::point > &points)
Write list of points to file.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const cellShapeControl & cellShapeControls() const
Return the cellShapeControl object.
static const fileName null
An empty fileName.
Definition: fileName.H:97
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:427
engineTime & runTime
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
bool distribute(const boundBox &bb)
scalar faceAreaRatioCoeff() const
Return the faceAreaRatioCoeff.
Definition: cvControlsI.H:185
static void timeCheck(const Time &runTime, const string &description=string::null, const bool check=true)
Write the elapsedCpuTime and memory usage, with an optional.
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:626
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const conformationSurfaces & geometryToConformTo() const
Return the conformationSurfaces object.
void drawDelaunayCell(Ostream &os, const CellHandle &c, label offset=0)
Draws a tet cell to an output stream. The offset is supplied as the tet.
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:51
scalar removalDistCoeff() const
Return removalDistCoeff.
Definition: cvControlsI.H:197
void writeMesh(const fileName &instance)
Prepare data and call writeMesh for polyMesh and.
void printVertexInfo(Ostream &os) const
Write vertex statistics in the form of a table to stream.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Definition: indexedVertex.H:51
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
void move()
Move the vertices according to the controller, re-conforming to the.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
bool uninitialised(const VertexType &v)
Type gSum(const FieldField< Field, Type > &f)
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
const pointField & points
void updateCellSizeAndFaceArea(vector &alignmentDir, scalar &targetFaceArea, scalar &targetCellSize) const
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
bool findSurfaceAnyIntersection(const point &start, const point &end) const
timeIndices insert(timeIndex, timeDirs[timeI].value())
const cellAspectRatioControl & aspectRatio() const
sideVolumeType
Normals point to the outside.
void reset()
Clear the entire triangulation.
PointFrompoint toPoint(const Foam::point &p)
word timeName
Definition: getTimeIndex.H:3
List< label > labelList
A List of labels.
Definition: labelList.H:56
void writeFeatureObj(const fileName &prefix) const
Write all components of all the extendedFeatureEdgeMeshes as.
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
label readLabel(Istream &is)
Definition: label.H:64
Istream and Ostream manipulators taking arguments.
const backgroundMeshDecomposition & decomposition() const
Return the backgroundMeshDecomposition.
void reverse(UList< T > &, const label n)
Definition: UListI.H:322
static const char nl
Definition: Ostream.H:265
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
defineTypeNameAndDebug(combustionModel, 0)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
scalar cosInsertionAcceptanceAngle() const
Return the cosInsertionAcceptanceAngle.
Definition: cvControlsI.H:191
void findSurfaceAllIntersections(const point &start, const point &end, List< pointIndexHit > &surfHit, labelList &hitSurface) const
void writeBoundaryPoints(const fileName &fName, const Triangulation &t)
Write the boundary Delaunay points to an OBJ file.
void resetCellCount()
Set the cell count to zero.
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:397
Foam::scalar averageCellSize(const VertexType &vA, const VertexType &vB)
Return the target cell size from that stored on a pair of Delaunay vertices,.
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
void updateDeltaVector(const vector &alignmentDir, const scalar targetCellSize, const scalar rABMag, vector &delta) const
label getNewCellIndex() const
Create a new unique cell index and return.
Definition: DelaunayMeshI.H:63
messageStream Info
static const NamedEnum< dualMeshPointType, 5 > dualMeshPointTypeNames_
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:98
const cvControls & foamyHexMeshControls() const
Return the foamyHexMeshControls object.
const Time & time() const
Return the Time object.
scalar coplanarTet(Cell &c, const scalar tol=1e-12)
Map< label > insertPoints(const List< Vb > &vertices, const bool reIndex)
Insert the list of vertices (calls rangeInsertWithInfo)
scalar insertionDistCoeff() const
Return the insertionDistCoeff.
Definition: cvControlsI.H:179
void sync(const boundBox &bb)
Refer vertices so that the processor interfaces are consistent.
InfoProxy< IOstream > info() const
Return info proxy.
Definition: IOstream.H:531
Namespace for OpenFOAM.
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49