35 const Foam::scalar Foam::conformalVoronoiMesh::searchConeAngle
38 const Foam::scalar Foam::conformalVoronoiMesh::searchAngleOppositeSurface
44 void Foam::conformalVoronoiMesh::conformToSurface()
46 this->resetCellCount();
50 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
51 cit != finite_cells_end();
58 if (!reconformToSurface())
61 reinsertSurfaceConformation();
65 sync(decomposition().procBounds());
73 buildSurfaceConformation();
75 if (distributeBackground(*
this))
79 sync(decomposition().procBounds());
85 storeSurfaceConformation();
92 bool Foam::conformalVoronoiMesh::reconformToSurface()
const 97 % foamyHexMeshControls().surfaceConformationRebuildFrequency() == 0
108 Foam::label Foam::conformalVoronoiMesh::findVerticesNearBoundaries()
110 label countNearBoundaryVertices = 0;
114 Delaunay::Finite_facets_iterator fit = finite_facets_begin();
115 fit != finite_facets_end();
119 Cell_handle c1 = fit->first;
120 Cell_handle c2 = fit->first->neighbor(fit->second);
122 if (is_infinite(c1) || is_infinite(c2))
130 if (!geometryToConformTo_.findSurfaceAnyIntersection(dE0, dE1))
135 for (
label celli = 0; celli < 4; ++celli)
137 Vertex_handle v = c1->vertex(celli);
142 && v->internalPoint()
143 && fit->second != celli
150 for (
label celli = 0; celli < 4; ++celli)
152 Vertex_handle v = c2->vertex(celli);
157 && v->internalPoint()
158 && fit->second != celli
161 v->setNearBoundary();
168 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
169 vit != finite_vertices_end();
173 if (vit->nearBoundary())
175 countNearBoundaryVertices++;
216 return countNearBoundaryVertices;
220 void Foam::conformalVoronoiMesh::buildSurfaceConformation()
222 timeCheck(
"Start buildSurfaceConformation");
225 <<
"Rebuilding surface conformation for more iterations" 228 existingEdgeLocations_.clearStorage();
229 existingSurfacePtLocations_.clearStorage();
231 buildEdgeLocationTree(existingEdgeLocations_);
232 buildSurfacePtLocationTree(existingSurfacePtLocations_);
234 label initialTotalHits = 0;
265 label countNearBoundaryVertices = findVerticesNearBoundaries();
267 Info<<
" Vertices marked as being near a boundary: " 268 <<
returnReduce(countNearBoundaryVertices, sumOp<label>())
269 <<
" (estimated)" <<
endl;
271 timeCheck(
"After set near boundary");
273 const scalar edgeSearchDistCoeffSqr =
274 foamyHexMeshControls().edgeSearchDistCoeffSqr();
276 const scalar surfacePtReplaceDistCoeffSqr =
277 foamyHexMeshControls().surfacePtReplaceDistCoeffSqr();
283 pointIndexHitAndFeatureDynList featureEdgeHits(AtoV/4);
284 pointIndexHitAndFeatureDynList surfaceHits(AtoV);
285 DynamicList<label> edgeToTreeShape(AtoV/4);
286 DynamicList<label> surfaceToTreeShape(AtoV);
288 Map<scalar> surfacePtToEdgePtDist(AtoV/4);
292 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
293 vit != finite_vertices_end();
297 if (vit->nearBoundary())
299 pointIndexHitAndFeatureDynList surfaceIntersections(AtoV);
303 dualCellSurfaceAllIntersections
314 addSurfaceAndEdgeHits
317 surfaceIntersections,
318 surfacePtReplaceDistCoeffSqr,
319 edgeSearchDistCoeffSqr,
324 surfacePtToEdgePtDist,
331 countNearBoundaryVertices--;
336 Info<<
" Vertices marked as being near a boundary: " 337 <<
returnReduce(countNearBoundaryVertices, sumOp<label>())
338 <<
" (after dual surface intersection)" <<
endl;
340 label nVerts = number_of_vertices();
341 label nSurfHits = surfaceHits.size();
342 label nFeatEdHits = featureEdgeHits.size();
346 reduce(nVerts, sumOp<label>());
347 reduce(nSurfHits, sumOp<label>());
348 reduce(nFeatEdHits, sumOp<label>());
351 Info<<
nl <<
"Initial conformation" <<
nl 352 <<
" Number of vertices " << nVerts <<
nl 353 <<
" Number of surface hits " << nSurfHits <<
nl 354 <<
" Number of edge hits " << nFeatEdHits
360 synchroniseSurfaceTrees(surfaceToTreeShape, surfaceHits);
363 DynamicList<Vb> pts(2*surfaceHits.size() + 3*featureEdgeHits.size());
365 insertSurfacePointPairs
368 "surfaceConformationLocations_initial.obj",
375 synchroniseEdgeTrees(edgeToTreeShape, featureEdgeHits);
378 insertEdgePointGroups
381 "edgeConformationLocations_initial.obj",
387 Map<label> oldToNewIndices = insertPointPairs(pts,
true,
true);
390 ptPairs_.reIndex(oldToNewIndices);
396 timeCheck(
"After initial conformation");
398 initialTotalHits = nSurfHits + nFeatEdHits;
407 autoPtr<labelPairHashSet> receivedVertices;
411 forAll(referralVertices, proci)
431 decomposition_().procBounds(),
437 label iterationNo = 0;
439 label maxIterations = foamyHexMeshControls().maxConformationIterations();
441 scalar iterationToInitialHitRatioLimit =
442 foamyHexMeshControls().iterationToInitialHitRatioLimit();
444 label hitLimit =
label(iterationToInitialHitRatioLimit*initialTotalHits);
446 Info<<
nl <<
"Stopping iterations when: " <<
nl 447 <<
" total number of hits drops below " 448 << iterationToInitialHitRatioLimit
449 <<
" of initial hits (" << hitLimit <<
")" <<
nl 451 <<
" maximum number of iterations (" << maxIterations
457 label totalHits = initialTotalHits;
462 && totalHits >= hitLimit
463 && iterationNo < maxIterations
466 pointIndexHitAndFeatureDynList surfaceHits(0.5*AtoV);
467 pointIndexHitAndFeatureDynList featureEdgeHits(0.25*AtoV);
468 DynamicList<label> surfaceToTreeShape(AtoV/2);
469 DynamicList<label> edgeToTreeShape(AtoV/4);
471 Map<scalar> surfacePtToEdgePtDist;
475 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
476 vit != finite_vertices_end();
487 || vit->internalBoundaryPoint()
488 || (vit->internalOrBoundaryPoint() && vit->referred())
491 pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
498 dualCellLargestSurfaceProtrusion(vit, surfHit, hitSurface);
502 surfaceIntersections.append
504 pointIndexHitAndFeature(surfHit, hitSurface)
507 addSurfaceAndEdgeHits
510 surfaceIntersections,
511 surfacePtReplaceDistCoeffSqr,
512 edgeSearchDistCoeffSqr,
517 surfacePtToEdgePtDist,
525 if (vit->nearBoundary())
533 vit->externalBoundaryPoint()
534 || (vit->externalBoundaryPoint() && vit->referred())
537 pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
544 dualCellLargestSurfaceIncursion(vit, surfHit, hitSurface);
548 surfaceIntersections.append
550 pointIndexHitAndFeature(surfHit, hitSurface)
553 addSurfaceAndEdgeHits
556 surfaceIntersections,
557 surfacePtReplaceDistCoeffSqr,
558 edgeSearchDistCoeffSqr,
563 surfacePtToEdgePtDist,
570 label nVerts = number_of_vertices();
571 label nSurfHits = surfaceHits.size();
572 label nFeatEdHits = featureEdgeHits.size();
576 reduce(nVerts, sumOp<label>());
577 reduce(nSurfHits, sumOp<label>());
578 reduce(nFeatEdHits, sumOp<label>());
581 Info<<
nl <<
"Conformation iteration " << iterationNo <<
nl 582 <<
" Number of vertices " << nVerts <<
nl 583 <<
" Number of surface hits " << nSurfHits <<
nl 584 <<
" Number of edge hits " << nFeatEdHits
587 totalHits = nSurfHits + nFeatEdHits;
589 label nNotInserted = 0;
597 synchroniseSurfaceTrees(surfaceToTreeShape, surfaceHits);
602 2*surfaceHits.size() + 3*featureEdgeHits.size()
605 insertSurfacePointPairs
608 "surfaceConformationLocations_" +
name(iterationNo) +
".obj",
616 synchroniseEdgeTrees(edgeToTreeShape, featureEdgeHits);
619 insertEdgePointGroups
622 "edgeConformationLocations_" +
name(iterationNo) +
".obj",
628 Map<label> oldToNewIndices = insertPointPairs(pts,
true,
true);
631 ptPairs_.reIndex(oldToNewIndices);
639 decomposition_().procBounds(),
646 timeCheck(
"Conformation iteration " +
name(iterationNo));
650 if (iterationNo == maxIterations)
653 <<
"Maximum surface conformation iterations (" 654 << maxIterations <<
") reached." <<
endl;
657 if (totalHits <= nNotInserted)
659 Info<<
nl <<
"Total hits (" << totalHits
660 <<
") less than number of failed insertions (" << nNotInserted
661 <<
"), stopping iterations" <<
endl;
665 if (totalHits < hitLimit)
667 Info<<
nl <<
"Total hits (" << totalHits
668 <<
") less than limit (" << hitLimit
669 <<
"), stopping iterations" <<
endl;
673 edgeLocationTreePtr_.clear();
674 surfacePtLocationTreePtr_.clear();
678 Foam::label Foam::conformalVoronoiMesh::synchroniseSurfaceTrees
680 const DynamicList<label>& surfaceToTreeShape,
681 pointIndexHitAndFeatureList& surfaceHits
684 Info<<
" Surface tree synchronisation" <<
endl;
686 pointIndexHitAndFeatureDynList synchronisedSurfLocations
691 List<pointIndexHitAndFeatureDynList> procSurfLocations(
Pstream::nProcs());
700 label nStoppedInsertion = 0;
711 const pointIndexHitAndFeatureList& otherSurfEdges =
712 procSurfLocations[proci];
714 forAll(otherSurfEdges, peI)
716 const Foam::point& pt = otherSurfEdges[peI].first().hitPoint();
719 pointIsNearSurfaceLocation(pt, nearest);
722 pointIsNearFeatureEdgeLocation(pt, nearestEdge);
724 if (nearest.hit() || nearestEdge.hit())
728 if (!hits[proci].
found(peI))
730 hits[proci].insert(peI);
743 synchronisedSurfLocations.append(surfaceHits[eI]);
747 surfacePtLocationTreePtr_().remove(surfaceToTreeShape[eI]);
761 Info<<
" Not inserting total of " << nNotInserted <<
" locations" 764 surfaceHits = synchronisedSurfLocations;
770 Foam::label Foam::conformalVoronoiMesh::synchroniseEdgeTrees
772 const DynamicList<label>& edgeToTreeShape,
773 pointIndexHitAndFeatureList& featureEdgeHits
776 Info<<
" Edge tree synchronisation" <<
endl;
778 pointIndexHitAndFeatureDynList synchronisedEdgeLocations
780 featureEdgeHits.size()
783 List<pointIndexHitAndFeatureDynList> procEdgeLocations(
Pstream::nProcs());
792 label nStoppedInsertion = 0;
803 pointIndexHitAndFeatureList& otherProcEdges = procEdgeLocations[proci];
805 forAll(otherProcEdges, peI)
807 const Foam::point& pt = otherProcEdges[peI].first().hitPoint();
810 pointIsNearFeatureEdgeLocation(pt, nearest);
823 if (!hits[proci].
found(peI))
825 hits[proci].insert(peI);
834 forAll(featureEdgeHits, eI)
838 synchronisedEdgeLocations.append(featureEdgeHits[eI]);
842 edgeLocationTreePtr_().remove(edgeToTreeShape[eI]);
856 Info<<
" Not inserting total of " << nNotInserted <<
" locations" 859 featureEdgeHits = synchronisedEdgeLocations;
865 bool Foam::conformalVoronoiMesh::surfaceLocationConformsToInside
867 const pointIndexHitAndFeature& info
870 if (info.first().hit())
874 geometryToConformTo_.getNormal
877 List<pointIndexHit>(1, info.first()),
881 const vector& n = norm[0];
883 const scalar ppDist = pointPairDistance(info.first().hitPoint());
885 const Foam::point innerPoint = info.first().hitPoint() - ppDist*
n;
887 if (!geometryToConformTo_.inside(innerPoint))
899 bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
901 const Delaunay::Finite_vertices_iterator& vit
904 std::list<Facet> facets;
905 incident_facets(vit, std::back_inserter(facets));
909 std::list<Facet>::iterator fit=facets.begin();
916 is_infinite(fit->first)
917 || is_infinite(fit->first->neighbor(fit->second))
918 || !fit->first->hasInternalPoint()
919 || !fit->first->neighbor(fit->second)->hasInternalPoint()
926 Foam::point dE1 = fit->first->neighbor(fit->second)->dual();
933 bool inProc = clipLineToProc(
topoint(vit->point()), a, b);
939 && geometryToConformTo_.findSurfaceAnyIntersection(a, b)
947 if (geometryToConformTo_.findSurfaceAnyIntersection(dE0, dE1))
958 bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
960 const Delaunay::Finite_vertices_iterator& vit,
961 pointIndexHitAndFeatureDynList& infoList
964 bool flagIntersection =
false;
966 std::list<Facet> facets;
967 incident_facets(vit, std::back_inserter(facets));
971 std::list<Facet>::iterator fit = facets.begin();
978 is_infinite(fit->first)
979 || is_infinite(fit->first->neighbor(fit->second))
980 || !fit->first->hasInternalPoint()
981 || !fit->first->neighbor(fit->second)->hasInternalPoint()
990 Foam::point dE1 = fit->first->neighbor(fit->second)->dual();
993 label hitSurfaceIntersection = -1;
997 bool inProc = clipLineToProc(
topoint(vit->point()), dE0, dE1);
1005 geometryToConformTo_.findSurfaceNearestIntersection
1010 hitSurfaceIntersection
1013 if (infoIntersection.hit())
1017 geometryToConformTo_.getNormal
1019 hitSurfaceIntersection,
1020 List<pointIndexHit>(1, infoIntersection),
1024 const vector& n = norm[0];
1028 const plane
p(infoIntersection.hitPoint(),
n);
1030 const plane::ray r(vertex, n);
1032 const scalar d = p.normalIntersect(r);
1036 pointIndexHitAndFeature info;
1037 geometryToConformTo_.findSurfaceNearest
1040 4.0*
magSqr(newPoint - vertex),
1045 bool rejectPoint =
false;
1047 if (!surfaceLocationConformsToInside(info))
1052 if (!rejectPoint && info.first().hit())
1054 if (!infoList.empty())
1061 infoList[hitI].first().index()
1062 == info.first().index()
1070 = infoList[hitI].first().hitPoint();
1072 const scalar separationDistance =
1073 mag(p - info.first().hitPoint());
1075 const scalar minSepDist =
1078 foamyHexMeshControls().removalDistCoeff()
1085 if (separationDistance < minSepDist)
1096 if (!rejectPoint && info.first().hit())
1098 flagIntersection =
true;
1099 infoList.append(info);
1104 return flagIntersection;
1108 bool Foam::conformalVoronoiMesh::clipLineToProc
1115 bool inProc =
false;
1117 pointIndexHit findAnyIntersection = decomposition_().findLine(a, b);
1119 if (!findAnyIntersection.hit())
1139 b = findAnyIntersection.hitPoint();
1144 a = findAnyIntersection.hitPoint();
1152 void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
1154 const Delaunay::Finite_vertices_iterator& vit,
1156 label& hitSurfaceLargest
1161 hitSurfaceLargest = -1;
1163 std::list<Facet> facets;
1164 finite_incident_facets(vit, std::back_inserter(facets));
1168 scalar maxProtrusionDistance = maxSurfaceProtrusion(vert);
1172 std::list<Facet>::iterator fit = facets.begin();
1173 fit != facets.end();
1177 Cell_handle c1 = fit->first;
1178 Cell_handle c2 = fit->first->neighbor(fit->second);
1182 is_infinite(c1) || is_infinite(c2)
1184 !c1->internalOrBoundaryDualVertex()
1185 || !c2->internalOrBoundaryDualVertex()
1187 || !c1->real() || !c2->real()
1196 if (
magSqr(vert - c1->dual()) <
magSqr(vert - c2->dual()))
1204 >
magSqr(geometryToConformTo().globalBounds().
mag())
1213 geometryToConformTo_.findSurfaceNearestIntersection
1225 allGeometry_[hitSurface].getNormal
1227 List<pointIndexHit>(1, surfHit),
1231 const vector& n = norm[0];
1233 const scalar normalProtrusionDistance
1235 (endPt - surfHit.hitPoint()) & n
1238 if (normalProtrusionDistance > maxProtrusionDistance)
1240 const plane
p(surfHit.hitPoint(),
n);
1242 const plane::ray r(endPt, -n);
1244 const scalar d = p.normalIntersect(r);
1248 pointIndexHitAndFeature info;
1249 geometryToConformTo_.findSurfaceNearest
1252 4.0*
magSqr(newPoint - endPt),
1257 if (info.first().hit())
1261 surfaceLocationConformsToInside
1263 pointIndexHitAndFeature(info.first(), info.second())
1267 surfHitLargest = info.first();
1268 hitSurfaceLargest = info.second();
1270 maxProtrusionDistance = normalProtrusionDistance;
1281 surfHitLargest.hit()
1284 && !decomposition().positionOnThisProcessor(surfHitLargest.hitPoint())
1292 hitSurfaceLargest = -1;
1297 void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
1299 const Delaunay::Finite_vertices_iterator& vit,
1301 label& hitSurfaceLargest
1306 hitSurfaceLargest = -1;
1308 std::list<Facet> facets;
1309 finite_incident_facets(vit, std::back_inserter(facets));
1313 scalar minIncursionDistance = -maxSurfaceProtrusion(vert);
1317 std::list<Facet>::iterator fit = facets.begin();
1318 fit != facets.end();
1322 Cell_handle c1 = fit->first;
1323 Cell_handle c2 = fit->first->neighbor(fit->second);
1327 is_infinite(c1) || is_infinite(c2)
1329 !c1->internalOrBoundaryDualVertex()
1330 || !c2->internalOrBoundaryDualVertex()
1332 || !c1->real() || !c2->real()
1341 if (
magSqr(vert - c1->dual()) <
magSqr(vert - c2->dual()))
1349 >
magSqr(geometryToConformTo().globalBounds().
mag())
1358 geometryToConformTo_.findSurfaceNearestIntersection
1370 allGeometry_[hitSurface].getNormal
1372 List<pointIndexHit>(1, surfHit),
1376 const vector& n = norm[0];
1378 scalar normalIncursionDistance
1380 (endPt - surfHit.hitPoint()) & n
1383 if (normalIncursionDistance < minIncursionDistance)
1385 const plane
p(surfHit.hitPoint(),
n);
1387 const plane::ray r(endPt, n);
1389 const scalar d = p.normalIntersect(r);
1393 pointIndexHitAndFeature info;
1394 geometryToConformTo_.findSurfaceNearest
1397 4.0*
magSqr(newPoint - endPt),
1402 if (info.first().hit())
1406 surfaceLocationConformsToInside
1408 pointIndexHitAndFeature(info.first(), info.second())
1412 surfHitLargest = info.first();
1413 hitSurfaceLargest = info.second();
1415 minIncursionDistance = normalIncursionDistance;
1426 surfHitLargest.hit()
1429 && !decomposition().positionOnThisProcessor(surfHitLargest.hitPoint())
1437 hitSurfaceLargest = -1;
1442 void Foam::conformalVoronoiMesh::reportProcessorOccupancy()
1446 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1447 vit != finite_vertices_end();
1456 && !decomposition().positionOnThisProcessor(
topoint(vit->point()))
1459 Pout<<
topoint(vit->point()) <<
" is not on this processor " 1572 void Foam::conformalVoronoiMesh::limitDisplacement
1574 const Delaunay::Finite_vertices_iterator& vit,
1584 displacement =
Zero;
1596 if (!geometryToConformTo_.globalBounds().contains(dispPt))
1601 else if (geometryToConformTo_.findSurfaceAnyIntersection(pt, dispPt))
1612 scalar searchDistanceSqr =
sqr 1614 2*vit->targetCellSize()
1615 *foamyHexMeshControls().pointPairDistanceCoeff()
1618 geometryToConformTo_.findSurfaceNearest
1630 if (
magSqr(pt - surfHit.hitPoint()) <= searchDistanceSqr)
1633 displacement =
Zero;
1644 displacement *= 0.5;
1646 limitDisplacement(vit, displacement, callCount);
1651 Foam::scalar Foam::conformalVoronoiMesh::angleBetweenSurfacePoints
1658 label pAsurfaceHit = -1;
1660 const scalar searchDist = 5.0*targetCellSize(pA);
1662 geometryToConformTo_.findSurfaceNearest
1677 allGeometry_[pAsurfaceHit].getNormal
1679 List<pointIndexHit>(1, pAhit),
1683 const vector nA = norm[0];
1686 label pBsurfaceHit = -1;
1688 geometryToConformTo_.findSurfaceNearest
1701 allGeometry_[pBsurfaceHit].getNormal
1703 List<pointIndexHit>(1, pBhit),
1707 const vector nB = norm[0];
1713 bool Foam::conformalVoronoiMesh::nearSurfacePoint
1715 pointIndexHitAndFeature& pHit
1721 const bool closeToSurfacePt = pointIsNearSurfaceLocation(pt, closePoint);
1727 magSqr(pt - closePoint.hitPoint())
1728 >
sqr(pointPairDistance(pt))
1732 const scalar cosAngle =
1733 angleBetweenSurfacePoints(pt, closePoint.hitPoint());
1736 if (cosAngle < searchAngleOppositeSurface)
1739 label pCloseSurfaceHit = -1;
1741 const scalar searchDist = targetCellSize(closePoint.hitPoint());
1743 geometryToConformTo_.findSurfaceNearest
1745 closePoint.hitPoint(),
1753 allGeometry_[pCloseSurfaceHit].getNormal
1755 List<pointIndexHit>(1, pCloseHit),
1759 const vector& nA = norm[0];
1762 label oppositeSurfaceHit = -1;
1764 geometryToConformTo_.findSurfaceNearestIntersection
1766 closePoint.hitPoint() + 0.5*pointPairDistance(pt)*nA,
1767 closePoint.hitPoint() + 5*targetCellSize(pt)*nA,
1772 if (oppositeHit.hit())
1775 pHit.first() = oppositeHit;
1776 pHit.second() = oppositeSurfaceHit;
1778 return !closeToSurfacePt;
1783 return closeToSurfacePt;
1787 bool Foam::conformalVoronoiMesh::appendToSurfacePtTree
1792 label startIndex = existingSurfacePtLocations_.size();
1794 existingSurfacePtLocations_.append(pt);
1796 label endIndex = existingSurfacePtLocations_.size();
1798 return surfacePtLocationTreePtr_().insert(startIndex, endIndex);
1802 bool Foam::conformalVoronoiMesh::appendToEdgeLocationTree
1807 label startIndex = existingEdgeLocations_.size();
1809 existingEdgeLocations_.append(pt);
1811 label endIndex = existingEdgeLocations_.size();
1813 return edgeLocationTreePtr_().insert(startIndex, endIndex);
1818 Foam::conformalVoronoiMesh::nearestFeatureEdgeLocations
1823 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1826 = edgeLocationTreePtr_().findSphere(pt, exclusionRangeSqr);
1828 DynamicList<pointIndexHit> dynPointHit;
1832 label index = elems[elemI];
1835 = edgeLocationTreePtr_().shapes().shapePoints()[index];
1839 dynPointHit.append(nearHit);
1846 bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation
1851 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1854 = edgeLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
1860 bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation
1866 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1868 info = edgeLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
1874 bool Foam::conformalVoronoiMesh::pointIsNearSurfaceLocation
1881 pointIsNearSurfaceLocation(pt, info);
1887 bool Foam::conformalVoronoiMesh::pointIsNearSurfaceLocation
1893 const scalar exclusionRangeSqr = surfacePtExclusionDistanceSqr(pt);
1895 info = surfacePtLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
1901 bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
1909 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1911 bool closeToFeatureEdge =
1912 pointIsNearFeatureEdgeLocation(pt, nearestEdgeHit);
1914 if (closeToFeatureEdge)
1916 List<pointIndexHit> nearHits = nearestFeatureEdgeLocations(pt);
1925 label featureHit = -1;
1927 geometryToConformTo_.findEdgeNearest
1935 const extendedFeatureEdgeMesh& eMesh
1936 = geometryToConformTo_.features()[featureHit];
1938 const vector& edgeDir = eMesh.edgeDirections()[edgeHit.index()];
1940 const vector lineBetweenPoints = pt - info.hitPoint();
1942 const scalar cosAngle
1954 mag(cosAngle) < searchConeAngle
1955 && (
mag(lineBetweenPoints) > pointPairDistance(pt))
1960 closeToFeatureEdge =
false;
1964 closeToFeatureEdge =
true;
1970 return closeToFeatureEdge;
1974 void Foam::conformalVoronoiMesh::buildEdgeLocationTree
1976 const DynamicList<Foam::point>& existingEdgeLocations
1979 treeBoundBox overallBb
1981 geometryToConformTo_.globalBounds().extend(1
e-4)
1984 edgeLocationTreePtr_.reset
1986 new dynamicIndexedOctree<dynamicTreeDataPoint>
1988 dynamicTreeDataPoint(existingEdgeLocations),
1998 void Foam::conformalVoronoiMesh::buildSurfacePtLocationTree
2000 const DynamicList<Foam::point>& existingSurfacePtLocations
2003 treeBoundBox overallBb
2005 geometryToConformTo_.globalBounds().extend(1
e-4)
2008 surfacePtLocationTreePtr_.reset
2010 new dynamicIndexedOctree<dynamicTreeDataPoint>
2012 dynamicTreeDataPoint(existingSurfacePtLocations),
2022 void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
2025 const pointIndexHitAndFeatureDynList& surfaceIntersections,
2026 scalar surfacePtReplaceDistCoeffSqr,
2027 scalar edgeSearchDistCoeffSqr,
2028 pointIndexHitAndFeatureDynList& surfaceHits,
2029 pointIndexHitAndFeatureDynList& featureEdgeHits,
2030 DynamicList<label>& surfaceToTreeShape,
2031 DynamicList<label>& edgeToTreeShape,
2032 Map<scalar>& surfacePtToEdgePtDist,
2036 const scalar cellSize = targetCellSize(vit);
2037 const scalar cellSizeSqr =
sqr(cellSize);
2039 forAll(surfaceIntersections, sI)
2041 pointIndexHitAndFeature surfHitI = surfaceIntersections[sI];
2043 bool keepSurfacePoint =
true;
2045 if (!surfHitI.first().hit())
2050 const Foam::point& surfPt = surfHitI.first().hitPoint();
2052 bool isNearFeaturePt = nearFeaturePt(surfPt);
2054 bool isNearFeatureEdge = surfacePtNearFeatureEdge(surfPt);
2056 bool isNearSurfacePt = nearSurfacePoint(surfHitI);
2058 if (isNearFeaturePt || isNearSurfacePt || isNearFeatureEdge)
2060 keepSurfacePoint =
false;
2063 List<List<pointIndexHit>> edHitsByFeature;
2067 const scalar searchRadiusSqr = edgeSearchDistCoeffSqr*cellSizeSqr;
2069 geometryToConformTo_.findAllNearestEdges
2077 forAll(edHitsByFeature, i)
2079 const label featureHit = featuresHit[i];
2081 List<pointIndexHit>& edHits = edHitsByFeature[i];
2094 && !decomposition().positionOnThisProcessor(edPt)
2101 if (!nearFeaturePt(edPt))
2106 < surfacePtReplaceDistCoeffSqr*cellSizeSqr
2116 keepSurfacePoint =
false;
2128 !nearFeatureEdgeLocation(edHit, nearestEdgeHit)
2131 appendToEdgeLocationTree(edPt);
2133 edgeToTreeShape.append
2135 existingEdgeLocations_.size() - 1
2140 featureEdgeHits.append
2142 pointIndexHitAndFeature(edHit, featureHit)
2148 surfacePtToEdgePtDist.insert
2150 existingEdgeLocations_.size() - 1,
2156 label hitIndex = nearestEdgeHit.index();
2163 < surfacePtToEdgePtDist[hitIndex]
2166 featureEdgeHits[hitIndex] =
2167 pointIndexHitAndFeature(edHit, featureHit);
2169 existingEdgeLocations_[hitIndex] =
2171 surfacePtToEdgePtDist[hitIndex] =
2177 edgeLocationTreePtr_().remove(hitIndex);
2178 edgeLocationTreePtr_().insert
2190 if (keepSurfacePoint)
2192 surfaceHits.append(surfHitI);
2193 appendToSurfacePtTree(surfPt);
2194 surfaceToTreeShape.append(existingSurfacePtLocations_.size() - 1);
2206 void Foam::conformalVoronoiMesh::storeSurfaceConformation()
2208 Info<<
nl <<
"Storing surface conformation" <<
endl;
2210 surfaceConformationVertices_.clear();
2213 DynamicList<Vb> tempSurfaceVertices(number_of_vertices()/10);
2217 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
2218 vit != finite_vertices_end();
2227 && vit->boundaryPoint()
2228 && !vit->featurePoint()
2229 && !vit->constrained()
2232 tempSurfaceVertices.append
2245 tempSurfaceVertices.shrink();
2247 surfaceConformationVertices_.transfer(tempSurfaceVertices);
2252 label(surfaceConformationVertices_.size()),
2255 <<
" vertices" <<
nl << endl;
2259 void Foam::conformalVoronoiMesh::reinsertSurfaceConformation()
2261 Info<<
nl <<
"Reinserting stored surface conformation" <<
endl;
2263 Map<label> oldToNewIndices =
2264 insertPointPairs(surfaceConformationVertices_,
true,
true);
2266 ptPairs_.reIndex(oldToNewIndices);
2268 PackedBoolList selectedElems(surfaceConformationVertices_.size(),
true);
2270 forAll(surfaceConformationVertices_, vI)
2272 Vb& v = surfaceConformationVertices_[vI];
2277 if (iter != oldToNewIndices.end())
2279 const label newIndex = iter();
2287 selectedElems[vI] =
false;
2292 inplaceSubset<PackedBoolList, List<Vb>>
2295 surfaceConformationVertices_
void setNearBoundary()
Set the point to be near the boundary.
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
pointFromPoint topoint(const Point &P)
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Ostream & endl(Ostream &os)
Add newline and flush stream.
HashTable< label, label, Hash< label > >::const_iterator const_iterator
PointIndexHit< point > pointIndexHit
Vector< scalar > vector
A scalar version of the templated Vector.
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
dimensionedScalar cos(const dimensionedScalar &ds)
List< label > labelList
A List of labels.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
Istream and Ostream manipulators taking arguments.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
static bool & parRun()
Is this a parallel run?
static label nProcs(const label communicator=0)
Number of processes in parallel run.
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
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.
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.