57 { dualMeshPointType::surface,
"surface" },
58 { dualMeshPointType::featureEdge,
"featureEdge" },
59 { dualMeshPointType::featurePoint,
"featurePoint" },
60 { dualMeshPointType::constrained,
"constrained" },
66 void Foam::conformalVoronoiMesh::cellSizeMeshOverlapsBackground()
const
68 const cellShapeControlMesh& cellSizeMesh =
71 DynamicList<Foam::point> pts(number_of_vertices());
75 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
76 vit != finite_vertices_end();
80 if (vit->internalOrBoundaryPoint() && !vit->referred())
82 pts.append(
topoint(vit->point()));
88 boundBox cellSizeMeshBb = cellSizeMesh.bounds();
90 bool fullyContained =
true;
92 if (!cellSizeMeshBb.contains(bb))
94 Pout<<
"Triangulation not fully contained in cell size mesh."
97 Pout<<
"Cell Size Mesh Bounds = " << cellSizeMesh.bounds() <<
endl;
98 Pout<<
"foamyHexMesh Bounds = " << bb <<
endl;
100 fullyContained =
false;
103 reduce(fullyContained, andOp<unsigned int>());
105 Info<<
"Triangulation is "
106 << (fullyContained ?
"fully" :
"not fully")
107 <<
" contained in the cell size mesh"
112 void Foam::conformalVoronoiMesh::insertInternalPoints
129 List<Foam::point> transferPoints(
points.size());
140 decomposition_().distributePoints(transferPoints)
143 transferPoints.clear();
148 label nVert = number_of_vertices();
152 label nInserted(number_of_vertices() - nVert);
156 reduce(nInserted, sumOp<label>());
159 Info<<
" " << nInserted <<
" points inserted"
160 <<
", failed to insert " <<
nPoints - nInserted
162 << 100.0*(
nPoints - nInserted)/(nInserted + SMALL)
167 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
168 vit != finite_vertices_end();
174 vit->index() = getNewVertexIndex();
175 vit->type() = Vb::vtInternal;
190 autoPtr<mapDistribute> mapDist =
191 decomposition_().distributePoints(
vertices);
205 label preReinsertionSize(number_of_vertices());
207 Map<label> oldToNewIndices =
212 label(number_of_vertices()) - preReinsertionSize,
216 Info<<
" Reinserted " << nReinserted <<
" vertices out of "
220 return oldToNewIndices;
224 void Foam::conformalVoronoiMesh::insertSurfacePointPairs
226 const pointIndexHitAndFeatureList& surfaceHits,
227 const fileName fName,
236 const label featureIndex = surfaceHits[i].second();
238 allGeometry_[featureIndex].getNormal
240 List<pointIndexHit>(1, surfaceHit),
244 const vector& normal = norm[0];
246 const Foam::point& surfacePt(surfaceHit.hitPoint());
249 geometryToConformTo_.meshableSide(featureIndex, surfaceHit);
253 createBafflePointPair
255 pointPairDistance(surfacePt),
266 pointPairDistance(surfacePt),
277 pointPairDistance(surfacePt),
287 << meshableSide <<
", bad"
292 if (foamyHexMeshControls().objOutput() && !fName.empty())
299 void Foam::conformalVoronoiMesh::insertEdgePointGroups
301 const pointIndexHitAndFeatureList& edgeHits,
302 const fileName fName,
308 if (edgeHits[i].first().hit())
310 const extendedFeatureEdgeMesh& feMesh
312 geometryToConformTo_.features()[edgeHits[i].second()]
327 if (foamyHexMeshControls().objOutput() && !fName.empty())
334 bool Foam::conformalVoronoiMesh::nearFeaturePt(
const Foam::point& pt)
const
336 scalar exclusionRangeSqr = featurePointExclusionDistanceSqr(pt);
341 geometryToConformTo_.findFeaturePointNearest
353 bool Foam::conformalVoronoiMesh::surfacePtNearFeatureEdge
358 scalar exclusionRangeSqr = surfacePtExclusionDistanceSqr(pt);
363 geometryToConformTo_.findEdgeNearest
375 void Foam::conformalVoronoiMesh::insertInitialPoints()
377 Info<<
nl <<
"Inserting initial points" <<
endl;
379 timeCheck(
"Before initial points call");
381 List<Point> initPts = initialPointsMethod_->initialPoints();
383 timeCheck(
"After initial points call");
387 insertInternalPoints(initPts);
389 if (initialPointsMethod_->fixInitialPoints())
393 Finite_vertices_iterator vit = finite_vertices_begin();
394 vit != finite_vertices_end();
402 if (foamyHexMeshControls().objOutput())
406 time().
path()/
"initialPoints.obj",
414 void Foam::conformalVoronoiMesh::distribute()
421 DynamicList<Foam::point>
points(number_of_vertices());
422 DynamicList<Foam::indexedVertexEnum::vertexType> types
426 DynamicList<scalar> sizes(number_of_vertices());
427 DynamicList<tensor> alignments(number_of_vertices());
431 Finite_vertices_iterator vit = finite_vertices_begin();
432 vit != finite_vertices_end();
439 types.append(vit->type());
440 sizes.append(vit->targetCellSize());
441 alignments.append(vit->alignment());
445 autoPtr<mapDistribute> mapDist =
448 mapDist().distribute(types);
449 mapDist().distribute(sizes);
450 mapDist().distribute(alignments);
455 Info<<
nl <<
" Inserting distributed tessellation" <<
endl;
459 DynamicList<Vb> verticesToInsert(
points.size());
463 verticesToInsert.append
474 verticesToInsert.last().targetCellSize() = sizes[pI];
475 verticesToInsert.last().alignment() = alignments[pI];
478 this->rangeInsertWithInfo
480 verticesToInsert.begin(),
481 verticesToInsert.end(),
485 Info<<
" Total number of vertices after redistribution "
488 label(number_of_vertices()), sumOp<label>()
494 void Foam::conformalVoronoiMesh::buildCellSizeAndAlignmentMesh()
496 controlMeshRefinement meshRefinement
501 smoothAlignmentSolver meshAlignmentSmoother
503 cellShapeControl_.shapeControlMesh()
506 meshRefinement.initialMeshPopulation(decomposition_);
508 cellShapeControlMesh& cellSizeMesh = cellShapeControl_.shapeControlMesh();
512 if (!distributeBackground(cellSizeMesh))
515 cellSizeMesh.distribute(decomposition_());
519 const dictionary& motionControlDict
520 = foamyHexMeshControls().foamyHexMeshDict().subDict(
"motionControl");
522 const label nMaxIter =
523 motionControlDict.get<label>(
"maxRefinementIterations");
525 Info<<
"Maximum number of refinement iterations : " << nMaxIter <<
endl;
527 for (label i = 0; i < nMaxIter; ++i)
529 label nAdded = meshRefinement.refineMesh(decomposition_);
531 reduce(nAdded, sumOp<label>());
535 cellSizeMesh.distribute(decomposition_());
538 Info<<
" Iteration " << i
539 <<
" Added = " << nAdded <<
" points"
551 if (!distributeBackground(cellSizeMesh))
553 cellSizeMesh.distribute(decomposition_());
557 meshAlignmentSmoother.smoothAlignments
559 motionControlDict.get<label>(
"maxSmoothingIterations")
562 Info<<
"Background cell size and alignment mesh:" <<
endl;
563 cellSizeMesh.printInfo(
Info);
565 Info<<
"Triangulation is "
566 << (cellSizeMesh.is_valid() ?
"valid" :
"not valid!" )
569 if (foamyHexMeshControls().writeCellShapeControlMesh())
572 cellSizeMesh.
write();
575 if (foamyHexMeshControls().printVertexInfo())
577 cellSizeMesh.printVertexInfo(
Info);
590 void Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()
592 Info<<
nl <<
"Calculating target cell alignment and size" <<
endl;
596 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
597 vit != finite_vertices_end();
601 if (vit->internalOrBoundaryPoint())
605 cellShapeControls().cellSizeAndAlignment
608 vit->targetCellSize(),
616 Foam::face Foam::conformalVoronoiMesh::buildDualFace
618 const Delaunay::Finite_edges_iterator& eit
621 Cell_circulator ccStart = incident_cells(*eit);
622 Cell_circulator cc1 = ccStart;
623 Cell_circulator cc2 = cc1;
629 DynamicList<label> verticesOnFace;
631 label nUniqueVertices = 0;
637 cc1->hasFarPoint() || cc2->hasFarPoint()
638 || is_infinite(cc1) || is_infinite(cc2)
641 Cell_handle
c = eit->first;
642 Vertex_handle vA =
c->vertex(eit->second);
643 Vertex_handle vB =
c->vertex(eit->third);
649 <<
"Dual face uses circumcenter defined by a "
650 <<
"Delaunay tetrahedron with no internal "
651 <<
"or boundary points. Defining Delaunay edge ends: "
658 label cc1I = cc1->cellIndex();
659 label cc2I = cc2->cellIndex();
663 if (!verticesOnFace.found(cc1I))
668 verticesOnFace.append(cc1I);
675 }
while (cc1 != ccStart);
677 verticesOnFace.shrink();
679 if (verticesOnFace.size() >= 3 && nUniqueVertices < 3)
694 verticesOnFace.setSize(nUniqueVertices);
697 return face(verticesOnFace);
701 Foam::label Foam::conformalVoronoiMesh::maxFilterCount
703 const Delaunay::Finite_edges_iterator& eit
706 Cell_circulator ccStart = incident_cells(*eit);
707 Cell_circulator cc = ccStart;
713 if (cc->hasFarPoint())
715 Cell_handle
c = eit->first;
716 Vertex_handle vA =
c->vertex(eit->second);
717 Vertex_handle vB =
c->vertex(eit->third);
720 <<
"Dual face uses circumcenter defined by a "
721 <<
"Delaunay tetrahedron with no internal "
722 <<
"or boundary points. Defining Delaunay edge ends: "
728 if (cc->filterCount() > maxFC)
730 maxFC = cc->filterCount();
735 }
while (cc != ccStart);
741 bool Foam::conformalVoronoiMesh::ownerAndNeighbour
754 label dualCellIndexA = vA->index();
756 if (!vA->internalOrBoundaryPoint() || vA->referred())
758 if (!vA->constrained())
764 label dualCellIndexB = vB->index();
766 if (!vB->internalOrBoundaryPoint() || vB->referred())
768 if (!vB->constrained())
774 if (dualCellIndexA == -1 && dualCellIndexB == -1)
777 <<
"Attempting to create a face joining "
778 <<
"two unindexed dual cells "
781 else if (dualCellIndexA == -1 || dualCellIndexB == -1)
785 if (dualCellIndexA == -1)
787 owner = dualCellIndexB;
793 owner = dualCellIndexA;
800 if (dualCellIndexB > dualCellIndexA)
802 owner = dualCellIndexA;
803 neighbour = dualCellIndexB;
807 owner = dualCellIndexB;
808 neighbour = dualCellIndexA;
821 Foam::conformalVoronoiMesh::conformalVoronoiMesh
824 const dictionary& foamyHexMeshDict,
825 const fileName& decompDictFile
830 rndGen_(64293*Pstream::myProcNo()),
831 foamyHexMeshControls_(foamyHexMeshDict),
836 "cvSearchableSurfaces",
843 foamyHexMeshDict.subDict(
"geometry"),
844 foamyHexMeshDict.getOrDefault(
"singleRegionName", true)
851 foamyHexMeshDict.subDict(
"surfaceConformation")
856 ? new backgroundMeshDecomposition
860 geometryToConformTo_,
861 foamyHexMeshControls().foamyHexMeshDict().subDict
863 "backgroundMeshDecomposition"
872 foamyHexMeshControls_,
878 ftPtConformer_(*this),
879 edgeLocationTreePtr_(),
880 surfacePtLocationTreePtr_(),
881 surfaceConformationVertices_(),
884 initialPointsMethod::
New
886 foamyHexMeshDict.subDict(
"initialPoints"),
889 geometryToConformTo_,
898 foamyHexMeshDict.subDict(
"motionControl"),
904 faceAreaWeightModel::
New
906 foamyHexMeshDict.subDict(
"motionControl")
922 if (foamyHexMeshControls().objOutput())
924 geometryToConformTo_.writeFeatureObj(
"foamyHexMesh");
927 buildCellSizeAndAlignmentMesh();
929 insertInitialPoints();
931 insertFeaturePoints(
true);
933 setVertexSizeAndAlignment();
935 cellSizeMeshOverlapsBackground();
940 distributeBackground(*
this);
942 buildSurfaceConformation();
946 distributeBackground(*
this);
950 sync(decomposition_().procBounds());
955 storeSurfaceConformation();
961 cellSizeMeshOverlapsBackground();
963 if (foamyHexMeshControls().printVertexInfo())
965 printVertexInfo(
Info);
968 if (foamyHexMeshControls().objOutput())
972 time().
path()/
"internalPoints_" + time().
timeName() +
".obj",
987 new backgroundMeshDecomposition
991 geometryToConformTo_,
992 foamyHexMeshControls().foamyHexMeshDict().subDict
994 "backgroundMeshDecomposition"
1000 insertInitialPoints();
1002 insertFeaturePoints();
1007 distributeBackground(*
this);
1009 buildSurfaceConformation();
1013 distributeBackground(*
this);
1017 sync(decomposition_().procBounds());
1020 cellSizeMeshOverlapsBackground();
1022 if (foamyHexMeshControls().printVertexInfo())
1024 printVertexInfo(
Info);
1031 timeCheck(
"Start of move");
1033 scalar relaxation = relaxationModel_->relaxation();
1035 Info<<
nl <<
"Relaxation = " << relaxation <<
endl;
1037 pointField dualVertices(number_of_finite_cells());
1039 this->resetCellCount();
1044 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1045 cit != finite_cells_end();
1049 cit->cellIndex() = Cb::ctUnassigned;
1051 if (cit->anyInternalOrBoundaryDualVertex())
1053 cit->cellIndex() = getNewCellIndex();
1055 dualVertices[cit->cellIndex()] = cit->dual();
1058 if (cit->hasFarPoint())
1060 cit->cellIndex() = Cb::ctFar;
1064 dualVertices.setSize(cellCount());
1066 setVertexSizeAndAlignment();
1068 timeCheck(
"Determined sizes and alignments");
1070 Info<<
nl <<
"Determining vertex displacements" <<
endl;
1074 cartesianDirections[0] =
vector(1, 0, 0);
1075 cartesianDirections[1] =
vector(0, 1, 0);
1076 cartesianDirections[2] =
vector(0, 0, 1);
1080 number_of_vertices(),
1084 bitSet pointToBeRetained(number_of_vertices(),
true);
1086 DynamicList<Point> pointsToInsert(number_of_vertices());
1090 Delaunay::Finite_edges_iterator eit = finite_edges_begin();
1091 eit != finite_edges_end();
1095 Cell_handle
c = eit->first;
1096 Vertex_handle vA =
c->vertex(eit->second);
1097 Vertex_handle vB =
c->vertex(eit->third);
1102 vA->internalPoint() && !vA->referred()
1103 && vB->internalOrBoundaryPoint()
1106 vB->internalPoint() && !vB->referred()
1107 && vA->internalOrBoundaryPoint()
1114 Field<vector> alignmentDirsA
1116 vA->alignment().T() & cartesianDirections
1118 Field<vector> alignmentDirsB
1120 vB->alignment().T() & cartesianDirections
1123 Field<vector> alignmentDirs(alignmentDirsA);
1125 forAll(alignmentDirsA, aA)
1127 const vector& a = alignmentDirsA[aA];
1129 scalar maxDotProduct = 0.0;
1131 forAll(alignmentDirsB, aB)
1133 const vector&
b = alignmentDirsB[aB];
1135 const scalar dotProduct = a &
b;
1137 if (
mag(dotProduct) > maxDotProduct)
1139 maxDotProduct =
mag(dotProduct);
1141 alignmentDirs[aA] = a +
sign(dotProduct)*
b;
1143 alignmentDirs[aA].normalise();
1150 scalar rABMag =
mag(rAB);
1158 vA->internalPoint() && !vA->referred() && !vA->fixed()
1159 && vB->internalPoint() && !vB->referred() && !vB->fixed()
1169 pointToBeRetained.test(vA->index())
1170 && pointToBeRetained.test(vB->index())
1175 if (internalPointIsInside(pt))
1177 pointsToInsert.append(
toPoint(pt));
1182 if (vA->internalPoint() && !vA->referred() && !vA->fixed())
1184 pointToBeRetained.unset(vA->index());
1187 if (vB->internalPoint() && !vB->referred() && !vB->fixed())
1189 pointToBeRetained.unset(vB->index());
1197 forAll(alignmentDirs, aD)
1199 vector& alignmentDir = alignmentDirs[aD];
1201 scalar dotProd = rAB & alignmentDir;
1211 const scalar alignmentDotProd = dotProd/rABMag;
1216 > foamyHexMeshControls().cosAlignmentAcceptanceAngle()
1219 scalar targetCellSize =
1222 scalar targetFaceArea =
sqr(targetCellSize);
1224 const vector originalAlignmentDir = alignmentDir;
1227 cellShapeControls().aspectRatio().updateCellSizeAndFaceArea
1239 face dualFace = buildDualFace(eit);
1245 const scalar faceArea = dualFace.mag(dualVertices);
1248 cellShapeControls().aspectRatio().updateDeltaVector
1250 originalAlignmentDir,
1273 delta *= faceAreaWeightModel_->faceAreaWeight
1275 faceArea/targetFaceArea
1281 (vA->internalPoint() && vB->internalPoint())
1282 && (!vA->referred() || !vB->referred())
1290 > foamyHexMeshControls().insertionDistCoeff()
1293 > foamyHexMeshControls().faceAreaRatioCoeff()
1296 > foamyHexMeshControls().cosInsertionAcceptanceAngle()
1302 !geometryToConformTo_.findSurfaceAnyIntersection
1312 if (internalPointIsInside(newPt))
1318 decomposition().positionOnThisProcessor
1324 pointsToInsert.append(
toPoint(newPt));
1329 pointsToInsert.append(
toPoint(newPt));
1337 (vA->internalPoint() && !vA->referred())
1338 || (vB->internalPoint() && !vB->referred())
1341 < foamyHexMeshControls().removalDistCoeff()
1367 pointToBeRetained.test(vA->index())
1368 && pointToBeRetained.test(vB->index())
1373 if (internalPointIsInside(pt))
1375 pointsToInsert.append(
toPoint(pt));
1387 pointToBeRetained.unset(vA->index());
1397 pointToBeRetained.unset(vB->index());
1411 displacementAccumulator[vA->index()] += 2*
delta;
1415 displacementAccumulator[vA->index()] +=
delta;
1428 displacementAccumulator[vB->index()] -= 2*
delta;
1432 displacementAccumulator[vB->index()] -=
delta;
1441 Info<<
"Limit displacements" <<
endl;
1446 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1447 vit != finite_vertices_end();
1451 if (vit->internalPoint() && !vit->referred() && !vit->fixed())
1453 if (pointToBeRetained.test(vit->index()))
1458 displacementAccumulator[vit->index()]
1464 vector totalDisp =
gSum(displacementAccumulator);
1465 scalar totalDist =
gSum(
mag(displacementAccumulator));
1467 displacementAccumulator *= relaxation;
1469 Info<<
"Sum displacements" <<
endl;
1471 label nPointsToRetain = 0;
1472 label nPointsToRemove = 0;
1476 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1477 vit != finite_vertices_end();
1481 if (vit->internalPoint() && !vit->referred() && !vit->fixed())
1483 if (pointToBeRetained.test(vit->index()))
1497 + displacementAccumulator[vit->index()]
1500 if (internalPointIsInside(pt))
1502 pointsToInsert.append(
toPoint(pt));
1511 pointsToInsert.shrink();
1515 nPointsToRetain - nPointsToRemove,
1518 <<
" internal points are outside the domain. "
1519 <<
"They will not be inserted." <<
endl;
1522 if (foamyHexMeshControls().objOutput() && time().writeTime())
1524 Info<<
"Writing point displacement vectors to file." <<
endl;
1527 time().
path()/
"displacements_" + runTime_.timeName() +
".obj"
1532 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1533 vit != finite_vertices_end();
1537 if (vit->internalPoint() && !vit->referred())
1539 if (pointToBeRetained.test(vit->index()))
1544 << displacementAccumulator[vit->index()][0] <<
" "
1545 << displacementAccumulator[vit->index()][1] <<
" "
1546 << displacementAccumulator[vit->index()][2] <<
" "
1556 timeCheck(
"Displacement calculated");
1558 Info<<
nl<<
"Inserting displaced tessellation" <<
endl;
1560 insertFeaturePoints(
true);
1562 insertInternalPoints(pointsToInsert,
true);
1585 timeCheck(
"Internal points inserted");
1592 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1593 vit != finite_vertices_end();
1597 if (!vit->referred() && !usedIndices.insert(vit->index()))
1600 <<
"Index already used! Could not insert: " <<
nl
1609 if (foamyHexMeshControls().objOutput())
1613 time().
path()/
"internalPoints_" + time().
timeName() +
".obj",
1618 if (reconformToSurface())
1622 time().
path()/
"boundaryPoints_" + time().
timeName() +
".obj",
1628 time().
path()/
"internalBoundaryPoints_" + time().
timeName()
1637 time().
path()/
"externalBoundaryPoints_" + time().
timeName()
1644 OBJstream multipleIntersections
1646 "multipleIntersections_"
1653 Delaunay::Finite_edges_iterator eit = finite_edges_begin();
1654 eit != finite_edges_end();
1658 Cell_handle
c = eit->first;
1659 Vertex_handle vA =
c->vertex(eit->second);
1660 Vertex_handle vB =
c->vertex(eit->third);
1665 List<pointIndexHit> surfHits;
1668 geometryToConformTo().findSurfaceAllIntersections
1678 surfHits.size() >= 2
1680 surfHits.size() == 0
1682 (vA->externalBoundaryPoint()
1683 && vB->internalBoundaryPoint())
1684 || (vB->externalBoundaryPoint()
1685 && vA->internalBoundaryPoint())
1696 timeCheck(
"After conformToSurface");
1698 if (foamyHexMeshControls().printVertexInfo())
1700 printVertexInfo(
Info);
1703 if (time().writeTime())
1709 <<
"Total displacement = " << totalDisp <<
nl
1710 <<
"Total distance = " << totalDist <<
nl
1715 void Foam::conformalVoronoiMesh::checkCoPlanarCells()
const
1717 typedef CGAL::Exact_predicates_exact_constructions_kernel Kexact;
1718 typedef CGAL::Point_3<Kexact> PointExact;
1722 Pout<<
"Triangulation is invalid!" <<
endl;
1725 OFstream str(
"badCells.obj");
1731 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1732 cit != finite_cells_end();
1741 <<
" quality = " << quality <<
nl
1746 FixedList<PointExact, 4> cellVerticesExact(PointExact(0,0,0));
1747 forAll(cellVerticesExact, vI)
1749 cellVerticesExact[vI] = PointExact
1751 cit->vertex(vI)->point().x(),
1752 cit->vertex(vI)->point().y(),
1753 cit->vertex(vI)->point().z()
1757 PointExact synchronisedDual = CGAL::circumcenter<Kexact>
1759 cellVerticesExact[0],
1760 cellVerticesExact[1],
1761 cellVerticesExact[2],
1762 cellVerticesExact[3]
1767 CGAL::to_double(synchronisedDual.x()),
1768 CGAL::to_double(synchronisedDual.y()),
1769 CGAL::to_double(synchronisedDual.z())
1772 Info<<
"inexact = " << cit->dual() <<
nl
1773 <<
"exact = " << exactPt <<
endl;
1777 Pout<<
"There are " << badCells <<
" bad cells out of "
1778 << number_of_finite_cells() <<
endl;
1781 label nNonGabriel = 0;
1784 Delaunay::Finite_facets_iterator fit = finite_facets_begin();
1785 fit != finite_facets_end();
1789 if (!is_Gabriel(*fit))
1795 Pout<<
"There are " << nNonGabriel <<
" non-Gabriel faces out of "
1796 << number_of_finite_facets() <<
endl;