38 const Foam::scalar Foam::conformalVoronoiMesh::searchConeAngle
41 const Foam::scalar Foam::conformalVoronoiMesh::searchAngleOppositeSurface
47 void Foam::conformalVoronoiMesh::conformToSurface()
49 this->resetCellCount();
53 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
54 cit != finite_cells_end();
58 cit->cellIndex() = Cb::ctUnassigned;
61 if (!reconformToSurface())
64 reinsertSurfaceConformation();
68 sync(decomposition().procBounds());
76 buildSurfaceConformation();
78 if (distributeBackground(*
this))
82 sync(decomposition().procBounds());
88 storeSurfaceConformation();
95 bool Foam::conformalVoronoiMesh::reconformToSurface()
const
100 % foamyHexMeshControls().surfaceConformationRebuildFrequency() == 0
111 Foam::label Foam::conformalVoronoiMesh::findVerticesNearBoundaries()
113 label countNearBoundaryVertices = 0;
117 Delaunay::Finite_facets_iterator fit = finite_facets_begin();
118 fit != finite_facets_end();
122 Cell_handle
c1 = fit->first;
123 Cell_handle
c2 = fit->first->neighbor(fit->second);
125 if (is_infinite(
c1) || is_infinite(
c2))
133 if (!geometryToConformTo_.findSurfaceAnyIntersection(dE0, dE1))
138 for (label celli = 0; celli < 4; ++celli)
140 Vertex_handle v =
c1->vertex(celli);
145 && v->internalPoint()
146 && fit->second != celli
149 v->setNearBoundary();
153 for (label celli = 0; celli < 4; ++celli)
155 Vertex_handle v =
c2->vertex(celli);
160 && v->internalPoint()
161 && fit->second != celli
164 v->setNearBoundary();
171 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
172 vit != finite_vertices_end();
176 if (vit->nearBoundary())
178 countNearBoundaryVertices++;
219 return countNearBoundaryVertices;
223 void Foam::conformalVoronoiMesh::buildSurfaceConformation()
225 timeCheck(
"Start buildSurfaceConformation");
228 <<
"Rebuilding surface conformation for more iterations"
231 existingEdgeLocations_.clearStorage();
232 existingSurfacePtLocations_.clearStorage();
234 buildEdgeLocationTree(existingEdgeLocations_);
235 buildSurfacePtLocationTree(existingSurfacePtLocations_);
237 label initialTotalHits = 0;
268 label countNearBoundaryVertices = findVerticesNearBoundaries();
270 Info<<
" Vertices marked as being near a boundary: "
271 <<
returnReduce(countNearBoundaryVertices, sumOp<label>())
272 <<
" (estimated)" <<
endl;
274 timeCheck(
"After set near boundary");
276 const scalar edgeSearchDistCoeffSqr =
277 foamyHexMeshControls().edgeSearchDistCoeffSqr();
279 const scalar surfacePtReplaceDistCoeffSqr =
280 foamyHexMeshControls().surfacePtReplaceDistCoeffSqr();
282 const label AtoV = label(6/
Foam::pow(scalar(number_of_vertices()), 3));
286 pointIndexHitAndFeatureDynList featureEdgeHits(AtoV/4);
287 pointIndexHitAndFeatureDynList surfaceHits(AtoV);
288 DynamicList<label> edgeToTreeShape(AtoV/4);
289 DynamicList<label> surfaceToTreeShape(AtoV);
291 Map<scalar> surfacePtToEdgePtDist(AtoV/4);
295 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
296 vit != finite_vertices_end();
300 if (vit->nearBoundary())
302 pointIndexHitAndFeatureDynList surfaceIntersections(AtoV);
306 dualCellSurfaceAllIntersections
317 addSurfaceAndEdgeHits
320 surfaceIntersections,
321 surfacePtReplaceDistCoeffSqr,
322 edgeSearchDistCoeffSqr,
327 surfacePtToEdgePtDist,
334 countNearBoundaryVertices--;
339 Info<<
" Vertices marked as being near a boundary: "
340 <<
returnReduce(countNearBoundaryVertices, sumOp<label>())
341 <<
" (after dual surface intersection)" <<
endl;
343 label nVerts = number_of_vertices();
344 label nSurfHits = surfaceHits.size();
345 label nFeatEdHits = featureEdgeHits.size();
349 reduce(nVerts, sumOp<label>());
350 reduce(nSurfHits, sumOp<label>());
351 reduce(nFeatEdHits, sumOp<label>());
354 Info<<
nl <<
"Initial conformation" <<
nl
355 <<
" Number of vertices " << nVerts <<
nl
356 <<
" Number of surface hits " << nSurfHits <<
nl
357 <<
" Number of edge hits " << nFeatEdHits
363 synchroniseSurfaceTrees(surfaceToTreeShape, surfaceHits);
366 DynamicList<Vb> pts(2*surfaceHits.size() + 3*featureEdgeHits.size());
368 insertSurfacePointPairs
371 "surfaceConformationLocations_initial.obj",
378 synchroniseEdgeTrees(edgeToTreeShape, featureEdgeHits);
381 insertEdgePointGroups
384 "edgeConformationLocations_initial.obj",
390 Map<label> oldToNewIndices = insertPointPairs(pts,
true,
true);
393 ptPairs_.reIndex(oldToNewIndices);
399 timeCheck(
"After initial conformation");
401 initialTotalHits = nSurfHits + nFeatEdHits;
410 autoPtr<labelPairHashSet> receivedVertices;
414 forAll(referralVertices, proci)
426 receivedVertices.reset
434 decomposition_().procBounds(),
440 label iterationNo = 0;
442 label maxIterations = foamyHexMeshControls().maxConformationIterations();
444 scalar iterationToInitialHitRatioLimit =
445 foamyHexMeshControls().iterationToInitialHitRatioLimit();
447 label hitLimit = label(iterationToInitialHitRatioLimit*initialTotalHits);
449 Info<<
nl <<
"Stopping iterations when: " <<
nl
450 <<
" total number of hits drops below "
451 << iterationToInitialHitRatioLimit
452 <<
" of initial hits (" << hitLimit <<
")" <<
nl
454 <<
" maximum number of iterations (" << maxIterations
460 label totalHits = initialTotalHits;
465 && totalHits >= hitLimit
466 && iterationNo < maxIterations
469 pointIndexHitAndFeatureDynList surfaceHits(0.5*AtoV);
470 pointIndexHitAndFeatureDynList featureEdgeHits(0.25*AtoV);
471 DynamicList<label> surfaceToTreeShape(AtoV/2);
472 DynamicList<label> edgeToTreeShape(AtoV/4);
474 Map<scalar> surfacePtToEdgePtDist;
478 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
479 vit != finite_vertices_end();
490 || vit->internalBoundaryPoint()
491 || (vit->internalOrBoundaryPoint() && vit->referred())
494 pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
501 dualCellLargestSurfaceProtrusion(vit, surfHit, hitSurface);
505 surfaceIntersections.append
507 pointIndexHitAndFeature(surfHit, hitSurface)
510 addSurfaceAndEdgeHits
513 surfaceIntersections,
514 surfacePtReplaceDistCoeffSqr,
515 edgeSearchDistCoeffSqr,
520 surfacePtToEdgePtDist,
528 if (vit->nearBoundary())
536 vit->externalBoundaryPoint()
537 || (vit->externalBoundaryPoint() && vit->referred())
540 pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
547 dualCellLargestSurfaceIncursion(vit, surfHit, hitSurface);
551 surfaceIntersections.append
553 pointIndexHitAndFeature(surfHit, hitSurface)
556 addSurfaceAndEdgeHits
559 surfaceIntersections,
560 surfacePtReplaceDistCoeffSqr,
561 edgeSearchDistCoeffSqr,
566 surfacePtToEdgePtDist,
573 label nVerts = number_of_vertices();
574 label nSurfHits = surfaceHits.size();
575 label nFeatEdHits = featureEdgeHits.size();
579 reduce(nVerts, sumOp<label>());
580 reduce(nSurfHits, sumOp<label>());
581 reduce(nFeatEdHits, sumOp<label>());
584 Info<<
nl <<
"Conformation iteration " << iterationNo <<
nl
585 <<
" Number of vertices " << nVerts <<
nl
586 <<
" Number of surface hits " << nSurfHits <<
nl
587 <<
" Number of edge hits " << nFeatEdHits
590 totalHits = nSurfHits + nFeatEdHits;
592 label nNotInserted = 0;
600 synchroniseSurfaceTrees(surfaceToTreeShape, surfaceHits);
605 2*surfaceHits.size() + 3*featureEdgeHits.size()
608 insertSurfacePointPairs
611 "surfaceConformationLocations_" +
name(iterationNo) +
".obj",
619 synchroniseEdgeTrees(edgeToTreeShape, featureEdgeHits);
622 insertEdgePointGroups
625 "edgeConformationLocations_" +
name(iterationNo) +
".obj",
631 Map<label> oldToNewIndices = insertPointPairs(pts,
true,
true);
634 ptPairs_.reIndex(oldToNewIndices);
642 decomposition_().procBounds(),
649 timeCheck(
"Conformation iteration " +
name(iterationNo));
653 if (iterationNo == maxIterations)
656 <<
"Maximum surface conformation iterations ("
657 << maxIterations <<
") reached." <<
endl;
660 if (totalHits <= nNotInserted)
662 Info<<
nl <<
"Total hits (" << totalHits
663 <<
") less than number of failed insertions (" << nNotInserted
664 <<
"), stopping iterations" <<
endl;
668 if (totalHits < hitLimit)
670 Info<<
nl <<
"Total hits (" << totalHits
671 <<
") less than limit (" << hitLimit
672 <<
"), stopping iterations" <<
endl;
676 edgeLocationTreePtr_.clear();
677 surfacePtLocationTreePtr_.clear();
681 Foam::label Foam::conformalVoronoiMesh::synchroniseSurfaceTrees
683 const DynamicList<label>& surfaceToTreeShape,
684 pointIndexHitAndFeatureList& surfaceHits
687 Info<<
" Surface tree synchronisation" <<
endl;
689 pointIndexHitAndFeatureDynList synchronisedSurfLocations
694 List<pointIndexHitAndFeatureDynList> procSurfLocations(
Pstream::nProcs());
703 label nStoppedInsertion = 0;
714 const pointIndexHitAndFeatureList& otherSurfEdges =
715 procSurfLocations[proci];
717 forAll(otherSurfEdges, peI)
719 const Foam::point& pt = otherSurfEdges[peI].first().hitPoint();
722 pointIsNearSurfaceLocation(pt, nearest);
725 pointIsNearFeatureEdgeLocation(pt, nearestEdge);
727 if (nearest.hit() || nearestEdge.hit())
730 hits[proci].insert(peI);
742 synchronisedSurfLocations.append(surfaceHits[eI]);
746 surfacePtLocationTreePtr_().remove(surfaceToTreeShape[eI]);
758 const label nNotInserted =
returnReduce(nStoppedInsertion, sumOp<label>());
760 Info<<
" Not inserting total of " << nNotInserted <<
" locations"
763 surfaceHits = synchronisedSurfLocations;
769 Foam::label Foam::conformalVoronoiMesh::synchroniseEdgeTrees
771 const DynamicList<label>& edgeToTreeShape,
772 pointIndexHitAndFeatureList& featureEdgeHits
775 Info<<
" Edge tree synchronisation" <<
endl;
777 pointIndexHitAndFeatureDynList synchronisedEdgeLocations
779 featureEdgeHits.size()
782 List<pointIndexHitAndFeatureDynList> procEdgeLocations(
Pstream::nProcs());
791 label nStoppedInsertion = 0;
802 pointIndexHitAndFeatureList& otherProcEdges = procEdgeLocations[proci];
804 forAll(otherProcEdges, peI)
806 const Foam::point& pt = otherProcEdges[peI].first().hitPoint();
809 pointIsNearFeatureEdgeLocation(pt, nearest);
821 hits[proci].insert(peI);
829 forAll(featureEdgeHits, eI)
833 synchronisedEdgeLocations.append(featureEdgeHits[eI]);
837 edgeLocationTreePtr_().remove(edgeToTreeShape[eI]);
849 const label nNotInserted =
returnReduce(nStoppedInsertion, sumOp<label>());
851 Info<<
" Not inserting total of " << nNotInserted <<
" locations"
854 featureEdgeHits = synchronisedEdgeLocations;
860 bool Foam::conformalVoronoiMesh::surfaceLocationConformsToInside
862 const pointIndexHitAndFeature& info
865 if (info.first().hit())
869 geometryToConformTo_.getNormal
872 List<pointIndexHit>(1, info.first()),
878 const scalar ppDist = pointPairDistance(info.first().hitPoint());
880 const Foam::point innerPoint = info.first().hitPoint() - ppDist*
n;
882 if (!geometryToConformTo_.inside(innerPoint))
894 bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
896 const Delaunay::Finite_vertices_iterator& vit
899 std::list<Facet> facets;
900 incident_facets(vit, std::back_inserter(facets));
904 std::list<Facet>::iterator fit=facets.begin();
911 is_infinite(fit->first)
912 || is_infinite(fit->first->neighbor(fit->second))
913 || !fit->first->hasInternalPoint()
914 || !fit->first->neighbor(fit->second)->hasInternalPoint()
921 Foam::point dE1 = fit->first->neighbor(fit->second)->dual();
928 bool inProc = clipLineToProc(
topoint(vit->point()), a,
b);
934 && geometryToConformTo_.findSurfaceAnyIntersection(a,
b)
942 if (geometryToConformTo_.findSurfaceAnyIntersection(dE0, dE1))
953 bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
955 const Delaunay::Finite_vertices_iterator& vit,
956 pointIndexHitAndFeatureDynList& infoList
959 bool flagIntersection =
false;
961 std::list<Facet> facets;
962 incident_facets(vit, std::back_inserter(facets));
966 std::list<Facet>::iterator fit = facets.begin();
973 is_infinite(fit->first)
974 || is_infinite(fit->first->neighbor(fit->second))
975 || !fit->first->hasInternalPoint()
976 || !fit->first->neighbor(fit->second)->hasInternalPoint()
985 Foam::point dE1 = fit->first->neighbor(fit->second)->dual();
988 label hitSurfaceIntersection = -1;
992 bool inProc = clipLineToProc(
topoint(vit->point()), dE0, dE1);
1000 geometryToConformTo_.findSurfaceNearestIntersection
1005 hitSurfaceIntersection
1008 if (infoIntersection.hit())
1012 geometryToConformTo_.getNormal
1014 hitSurfaceIntersection,
1015 List<pointIndexHit>(1, infoIntersection),
1023 const plane
p(infoIntersection.hitPoint(),
n);
1025 const plane::ray r(vertex,
n);
1027 const scalar d =
p.normalIntersect(r);
1031 pointIndexHitAndFeature info;
1032 geometryToConformTo_.findSurfaceNearest
1035 4.0*
magSqr(newPoint - vertex),
1040 bool rejectPoint =
false;
1042 if (!surfaceLocationConformsToInside(info))
1047 if (!rejectPoint && info.first().hit())
1049 if (!infoList.empty())
1056 infoList[hitI].first().index()
1057 == info.first().index()
1065 = infoList[hitI].first().hitPoint();
1067 const scalar separationDistance =
1068 mag(
p - info.first().hitPoint());
1070 const scalar minSepDist =
1073 foamyHexMeshControls().removalDistCoeff()
1080 if (separationDistance < minSepDist)
1091 if (!rejectPoint && info.first().hit())
1093 flagIntersection =
true;
1094 infoList.append(info);
1099 return flagIntersection;
1103 bool Foam::conformalVoronoiMesh::clipLineToProc
1110 bool inProc =
false;
1112 pointIndexHit findAnyIntersection = decomposition_().findLine(a,
b);
1114 if (!findAnyIntersection.hit())
1134 b = findAnyIntersection.hitPoint();
1139 a = findAnyIntersection.hitPoint();
1147 void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
1149 const Delaunay::Finite_vertices_iterator& vit,
1151 label& hitSurfaceLargest
1156 hitSurfaceLargest = -1;
1158 std::list<Facet> facets;
1159 finite_incident_facets(vit, std::back_inserter(facets));
1163 scalar maxProtrusionDistance = maxSurfaceProtrusion(vert);
1167 std::list<Facet>::iterator fit = facets.begin();
1168 fit != facets.end();
1172 Cell_handle
c1 = fit->first;
1173 Cell_handle
c2 = fit->first->neighbor(fit->second);
1177 is_infinite(
c1) || is_infinite(
c2)
1179 !
c1->internalOrBoundaryDualVertex()
1180 || !
c2->internalOrBoundaryDualVertex()
1182 || !
c1->real() || !
c2->real()
1199 >
magSqr(geometryToConformTo().globalBounds().
mag())
1208 geometryToConformTo_.findSurfaceNearestIntersection
1220 allGeometry_[hitSurface].getNormal
1222 List<pointIndexHit>(1, surfHit),
1228 const scalar normalProtrusionDistance
1230 (endPt - surfHit.hitPoint()) &
n
1233 if (normalProtrusionDistance > maxProtrusionDistance)
1235 const plane
p(surfHit.hitPoint(),
n);
1237 const plane::ray r(endPt, -
n);
1239 const scalar d =
p.normalIntersect(r);
1243 pointIndexHitAndFeature info;
1244 geometryToConformTo_.findSurfaceNearest
1247 4.0*
magSqr(newPoint - endPt),
1252 if (info.first().hit())
1256 surfaceLocationConformsToInside
1258 pointIndexHitAndFeature(info.first(), info.second())
1262 surfHitLargest = info.first();
1263 hitSurfaceLargest = info.second();
1265 maxProtrusionDistance = normalProtrusionDistance;
1276 surfHitLargest.hit()
1279 && !decomposition().positionOnThisProcessor(surfHitLargest.hitPoint())
1287 hitSurfaceLargest = -1;
1292 void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
1294 const Delaunay::Finite_vertices_iterator& vit,
1296 label& hitSurfaceLargest
1301 hitSurfaceLargest = -1;
1303 std::list<Facet> facets;
1304 finite_incident_facets(vit, std::back_inserter(facets));
1308 scalar minIncursionDistance = -maxSurfaceProtrusion(vert);
1312 std::list<Facet>::iterator fit = facets.begin();
1313 fit != facets.end();
1317 Cell_handle
c1 = fit->first;
1318 Cell_handle
c2 = fit->first->neighbor(fit->second);
1322 is_infinite(
c1) || is_infinite(
c2)
1324 !
c1->internalOrBoundaryDualVertex()
1325 || !
c2->internalOrBoundaryDualVertex()
1327 || !
c1->real() || !
c2->real()
1344 >
magSqr(geometryToConformTo().globalBounds().
mag())
1353 geometryToConformTo_.findSurfaceNearestIntersection
1365 allGeometry_[hitSurface].getNormal
1367 List<pointIndexHit>(1, surfHit),
1373 scalar normalIncursionDistance
1375 (endPt - surfHit.hitPoint()) &
n
1378 if (normalIncursionDistance < minIncursionDistance)
1380 const plane
p(surfHit.hitPoint(),
n);
1382 const plane::ray r(endPt,
n);
1384 const scalar d =
p.normalIntersect(r);
1388 pointIndexHitAndFeature info;
1389 geometryToConformTo_.findSurfaceNearest
1392 4.0*
magSqr(newPoint - endPt),
1397 if (info.first().hit())
1401 surfaceLocationConformsToInside
1403 pointIndexHitAndFeature(info.first(), info.second())
1407 surfHitLargest = info.first();
1408 hitSurfaceLargest = info.second();
1410 minIncursionDistance = normalIncursionDistance;
1421 surfHitLargest.hit()
1424 && !decomposition().positionOnThisProcessor(surfHitLargest.hitPoint())
1432 hitSurfaceLargest = -1;
1437 void Foam::conformalVoronoiMesh::reportProcessorOccupancy()
1441 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1442 vit != finite_vertices_end();
1451 && !decomposition().positionOnThisProcessor(
topoint(vit->point()))
1454 Pout<<
topoint(vit->point()) <<
" is not on this processor "
1567 void Foam::conformalVoronoiMesh::limitDisplacement
1569 const Delaunay::Finite_vertices_iterator& vit,
1579 displacement =
Zero;
1591 if (!geometryToConformTo_.globalBounds().contains(dispPt))
1596 else if (geometryToConformTo_.findSurfaceAnyIntersection(pt, dispPt))
1607 scalar searchDistanceSqr =
sqr
1609 2*vit->targetCellSize()
1610 *foamyHexMeshControls().pointPairDistanceCoeff()
1613 geometryToConformTo_.findSurfaceNearest
1625 if (
magSqr(pt - surfHit.hitPoint()) <= searchDistanceSqr)
1628 displacement =
Zero;
1639 displacement *= 0.5;
1641 limitDisplacement(vit, displacement, callCount);
1646 Foam::scalar Foam::conformalVoronoiMesh::angleBetweenSurfacePoints
1653 label pAsurfaceHit = -1;
1655 const scalar searchDist = 5.0*targetCellSize(pA);
1657 geometryToConformTo_.findSurfaceNearest
1672 allGeometry_[pAsurfaceHit].getNormal
1674 List<pointIndexHit>(1, pAhit),
1678 const vector nA = norm[0];
1681 label pBsurfaceHit = -1;
1683 geometryToConformTo_.findSurfaceNearest
1696 allGeometry_[pBsurfaceHit].getNormal
1698 List<pointIndexHit>(1, pBhit),
1702 const vector nB = norm[0];
1708 bool Foam::conformalVoronoiMesh::nearSurfacePoint
1710 pointIndexHitAndFeature& pHit
1716 const bool closeToSurfacePt = pointIsNearSurfaceLocation(pt, closePoint);
1722 magSqr(pt - closePoint.hitPoint())
1723 >
sqr(pointPairDistance(pt))
1727 const scalar cosAngle =
1728 angleBetweenSurfacePoints(pt, closePoint.hitPoint());
1731 if (cosAngle < searchAngleOppositeSurface)
1734 label pCloseSurfaceHit = -1;
1736 const scalar searchDist = targetCellSize(closePoint.hitPoint());
1738 geometryToConformTo_.findSurfaceNearest
1740 closePoint.hitPoint(),
1748 allGeometry_[pCloseSurfaceHit].getNormal
1750 List<pointIndexHit>(1, pCloseHit),
1754 const vector& nA = norm[0];
1757 label oppositeSurfaceHit = -1;
1759 geometryToConformTo_.findSurfaceNearestIntersection
1761 closePoint.hitPoint() + 0.5*pointPairDistance(pt)*nA,
1762 closePoint.hitPoint() + 5*targetCellSize(pt)*nA,
1767 if (oppositeHit.hit())
1770 pHit.first() = oppositeHit;
1771 pHit.second() = oppositeSurfaceHit;
1773 return !closeToSurfacePt;
1778 return closeToSurfacePt;
1782 bool Foam::conformalVoronoiMesh::appendToSurfacePtTree
1787 label startIndex = existingSurfacePtLocations_.size();
1789 existingSurfacePtLocations_.append(pt);
1791 label endIndex = existingSurfacePtLocations_.size();
1793 return surfacePtLocationTreePtr_().insert(startIndex, endIndex);
1797 bool Foam::conformalVoronoiMesh::appendToEdgeLocationTree
1802 label startIndex = existingEdgeLocations_.size();
1804 existingEdgeLocations_.append(pt);
1806 label endIndex = existingEdgeLocations_.size();
1808 return edgeLocationTreePtr_().insert(startIndex, endIndex);
1813 Foam::conformalVoronoiMesh::nearestFeatureEdgeLocations
1818 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1821 = edgeLocationTreePtr_().findSphere(pt, exclusionRangeSqr);
1823 DynamicList<pointIndexHit> dynPointHit;
1827 label index = elems[elemI];
1830 = edgeLocationTreePtr_().shapes().shapePoints()[index];
1834 dynPointHit.append(nearHit);
1841 bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation
1846 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1849 = edgeLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
1855 bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation
1861 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1863 info = edgeLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
1869 bool Foam::conformalVoronoiMesh::pointIsNearSurfaceLocation
1876 pointIsNearSurfaceLocation(pt, info);
1882 bool Foam::conformalVoronoiMesh::pointIsNearSurfaceLocation
1888 const scalar exclusionRangeSqr = surfacePtExclusionDistanceSqr(pt);
1890 info = surfacePtLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
1896 bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
1904 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1906 bool closeToFeatureEdge =
1907 pointIsNearFeatureEdgeLocation(pt, nearestEdgeHit);
1909 if (closeToFeatureEdge)
1911 List<pointIndexHit> nearHits = nearestFeatureEdgeLocations(pt);
1920 label featureHit = -1;
1922 geometryToConformTo_.findEdgeNearest
1930 const extendedFeatureEdgeMesh& eMesh
1931 = geometryToConformTo_.features()[featureHit];
1933 const vector& edgeDir = eMesh.edgeDirections()[edgeHit.index()];
1935 const vector lineBetweenPoints = pt - info.hitPoint();
1937 const scalar cosAngle
1949 mag(cosAngle) < searchConeAngle
1950 && (
mag(lineBetweenPoints) > pointPairDistance(pt))
1955 closeToFeatureEdge =
false;
1959 closeToFeatureEdge =
true;
1965 return closeToFeatureEdge;
1969 void Foam::conformalVoronoiMesh::buildEdgeLocationTree
1971 const DynamicList<Foam::point>& existingEdgeLocations
1974 treeBoundBox overallBb
1976 geometryToConformTo_.globalBounds().extend(rndGen_, 1
e-4)
1982 edgeLocationTreePtr_.reset
1984 new dynamicIndexedOctree<dynamicTreeDataPoint>
1986 dynamicTreeDataPoint(existingEdgeLocations),
1996 void Foam::conformalVoronoiMesh::buildSurfacePtLocationTree
1998 const DynamicList<Foam::point>& existingSurfacePtLocations
2001 treeBoundBox overallBb
2003 geometryToConformTo_.globalBounds().extend(rndGen_, 1
e-4)
2009 surfacePtLocationTreePtr_.reset
2011 new dynamicIndexedOctree<dynamicTreeDataPoint>
2013 dynamicTreeDataPoint(existingSurfacePtLocations),
2023 void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
2026 const pointIndexHitAndFeatureDynList& surfaceIntersections,
2027 scalar surfacePtReplaceDistCoeffSqr,
2028 scalar edgeSearchDistCoeffSqr,
2029 pointIndexHitAndFeatureDynList& surfaceHits,
2030 pointIndexHitAndFeatureDynList& featureEdgeHits,
2031 DynamicList<label>& surfaceToTreeShape,
2032 DynamicList<label>& edgeToTreeShape,
2033 Map<scalar>& surfacePtToEdgePtDist,
2037 const scalar cellSize = targetCellSize(vit);
2038 const scalar cellSizeSqr =
sqr(cellSize);
2040 forAll(surfaceIntersections, sI)
2042 pointIndexHitAndFeature surfHitI = surfaceIntersections[sI];
2044 bool keepSurfacePoint =
true;
2046 if (!surfHitI.first().hit())
2051 const Foam::point& surfPt = surfHitI.first().hitPoint();
2053 bool isNearFeaturePt = nearFeaturePt(surfPt);
2055 bool isNearFeatureEdge = surfacePtNearFeatureEdge(surfPt);
2057 bool isNearSurfacePt = nearSurfacePoint(surfHitI);
2059 if (isNearFeaturePt || isNearSurfacePt || isNearFeatureEdge)
2061 keepSurfacePoint =
false;
2064 List<List<pointIndexHit>> edHitsByFeature;
2068 const scalar searchRadiusSqr = edgeSearchDistCoeffSqr*cellSizeSqr;
2070 geometryToConformTo_.findAllNearestEdges
2078 forAll(edHitsByFeature, i)
2080 const label featureHit = featuresHit[i];
2082 List<pointIndexHit>& edHits = edHitsByFeature[i];
2095 && !decomposition().positionOnThisProcessor(edPt)
2102 if (!nearFeaturePt(edPt))
2107 < surfacePtReplaceDistCoeffSqr*cellSizeSqr
2117 keepSurfacePoint =
false;
2129 !nearFeatureEdgeLocation(edHit, nearestEdgeHit)
2132 appendToEdgeLocationTree(edPt);
2134 edgeToTreeShape.append
2136 existingEdgeLocations_.size() - 1
2141 featureEdgeHits.append
2143 pointIndexHitAndFeature(edHit, featureHit)
2149 surfacePtToEdgePtDist.insert
2151 existingEdgeLocations_.size() - 1,
2157 label hitIndex = nearestEdgeHit.index();
2164 < surfacePtToEdgePtDist[hitIndex]
2167 featureEdgeHits[hitIndex] =
2168 pointIndexHitAndFeature(edHit, featureHit);
2170 existingEdgeLocations_[hitIndex] =
2172 surfacePtToEdgePtDist[hitIndex] =
2178 edgeLocationTreePtr_().remove(hitIndex);
2179 edgeLocationTreePtr_().insert
2191 if (keepSurfacePoint)
2193 surfaceHits.append(surfHitI);
2194 appendToSurfacePtTree(surfPt);
2195 surfaceToTreeShape.append(existingSurfacePtLocations_.size() - 1);
2207 void Foam::conformalVoronoiMesh::storeSurfaceConformation()
2209 Info<<
nl <<
"Storing surface conformation" <<
endl;
2211 surfaceConformationVertices_.clear();
2214 DynamicList<Vb> tempSurfaceVertices(number_of_vertices()/10);
2218 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
2219 vit != finite_vertices_end();
2228 && vit->boundaryPoint()
2229 && !vit->featurePoint()
2230 && !vit->constrained()
2233 tempSurfaceVertices.append
2246 tempSurfaceVertices.shrink();
2248 surfaceConformationVertices_.transfer(tempSurfaceVertices);
2253 label(surfaceConformationVertices_.size()),
2256 <<
" vertices" <<
nl <<
endl;
2260 void Foam::conformalVoronoiMesh::reinsertSurfaceConformation()
2262 Info<<
nl <<
"Reinserting stored surface conformation" <<
endl;
2264 Map<label> oldToNewIndices =
2265 insertPointPairs(surfaceConformationVertices_,
true,
true);
2267 ptPairs_.reIndex(oldToNewIndices);
2269 bitSet selectedElems(surfaceConformationVertices_.size(),
true);
2271 forAll(surfaceConformationVertices_, vI)
2273 Vb& v = surfaceConformationVertices_[vI];
2274 label& vIndex = v.
index();
2276 const auto iter = oldToNewIndices.cfind(vIndex);
2280 const label newIndex = *iter;
2288 selectedElems.unset(vI);
2293 inplaceSubset<bitSet, List<Vb>>
2296 surfaceConformationVertices_