46 this->resetCellCount();
50 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
51 cit != finite_cells_end();
55 cit->cellIndex() = Cb::ctUnassigned;
58 if (!reconformToSurface())
61 reinsertSurfaceConformation();
73 buildSurfaceConformation();
75 if (distributeBackground(*
this))
85 storeSurfaceConformation();
97 % foamyHexMeshControls().surfaceConformationRebuildFrequency() == 0
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
146 v->setNearBoundary();
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;
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)
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);
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();
680 const DynamicList<label>& surfaceToTreeShape,
681 pointIndexHitAndFeatureList& surfaceHits
684 Info<<
" Surface tree synchronisation" <<
endl;
686 pointIndexHitAndFeatureDynList synchronisedSurfLocations
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;
772 const DynamicList<label>& edgeToTreeShape,
773 pointIndexHitAndFeatureList& featureEdgeHits
776 Info<<
" Edge tree synchronisation" <<
endl;
778 pointIndexHitAndFeatureDynList synchronisedEdgeLocations
780 featureEdgeHits.size()
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;
867 const pointIndexHitAndFeature& info
870 if (info.first().hit())
874 geometryToConformTo_.getNormal
883 const scalar ppDist = pointPairDistance(info.first().hitPoint());
885 const Foam::point innerPoint = info.first().hitPoint() - ppDist*
n;
887 if (!geometryToConformTo_.inside(innerPoint))
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))
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,
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;
1115 bool inProc =
false;
1117 pointIndexHit findAnyIntersection = decomposition_().findLine(a,
b);
1119 if (!findAnyIntersection.hit())
1139 b = findAnyIntersection.hitPoint();
1144 a = findAnyIntersection.hitPoint();
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()
1204 >
magSqr(geometryToConformTo().globalBounds().
mag())
1213 geometryToConformTo_.findSurfaceNearestIntersection
1225 allGeometry_[hitSurface].getNormal
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;
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()
1349 >
magSqr(geometryToConformTo().globalBounds().
mag())
1358 geometryToConformTo_.findSurfaceNearestIntersection
1370 allGeometry_[hitSurface].getNormal
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;
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 "
1574 const Delaunay::Finite_vertices_iterator& vit,
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)
1644 displacement *= 0.5;
1646 limitDisplacement(vit, displacement, callCount);
1658 label pAsurfaceHit = -1;
1660 const scalar searchDist = 5.0*targetCellSize(pA);
1662 geometryToConformTo_.findSurfaceNearest
1677 allGeometry_[pAsurfaceHit].getNormal
1683 const vector nA = norm[0];
1686 label pBsurfaceHit = -1;
1688 geometryToConformTo_.findSurfaceNearest
1701 allGeometry_[pBsurfaceHit].getNormal
1707 const vector nB = norm[0];
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
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;
1792 label startIndex = existingSurfacePtLocations_.size();
1794 existingSurfacePtLocations_.append(pt);
1796 label endIndex = existingSurfacePtLocations_.size();
1798 return surfacePtLocationTreePtr_().insert(startIndex, endIndex);
1807 label startIndex = existingEdgeLocations_.size();
1809 existingEdgeLocations_.append(pt);
1811 label endIndex = existingEdgeLocations_.size();
1813 return edgeLocationTreePtr_().insert(startIndex, endIndex);
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);
1851 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1854 = edgeLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
1866 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1868 info = edgeLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
1881 pointIsNearSurfaceLocation(pt, info);
1893 const scalar exclusionRangeSqr = surfacePtExclusionDistanceSqr(pt);
1895 info = surfacePtLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
1909 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1911 bool closeToFeatureEdge =
1912 pointIsNearFeatureEdgeLocation(pt, nearestEdgeHit);
1914 if (closeToFeatureEdge)
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;
1976 const DynamicList<Foam::point>& existingEdgeLocations
1979 treeBoundBox overallBb
1981 geometryToConformTo_.globalBounds().extend(rndGen_, 1
e-4)
1984 overallBb.min() -=
Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
1985 overallBb.max() +=
Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
1987 edgeLocationTreePtr_.reset
1989 new dynamicIndexedOctree<dynamicTreeDataPoint>
1991 dynamicTreeDataPoint(existingEdgeLocations),
2003 const DynamicList<Foam::point>& existingSurfacePtLocations
2006 treeBoundBox overallBb
2008 geometryToConformTo_.globalBounds().extend(rndGen_, 1
e-4)
2011 overallBb.min() -=
Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
2012 overallBb.max() +=
Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
2014 surfacePtLocationTreePtr_.reset
2016 new dynamicIndexedOctree<dynamicTreeDataPoint>
2018 dynamicTreeDataPoint(existingSurfacePtLocations),
2031 const pointIndexHitAndFeatureDynList& surfaceIntersections,
2032 scalar surfacePtReplaceDistCoeffSqr,
2033 scalar edgeSearchDistCoeffSqr,
2034 pointIndexHitAndFeatureDynList& surfaceHits,
2035 pointIndexHitAndFeatureDynList& featureEdgeHits,
2036 DynamicList<label>& surfaceToTreeShape,
2037 DynamicList<label>& edgeToTreeShape,
2038 Map<scalar>& surfacePtToEdgePtDist,
2042 const scalar cellSize = targetCellSize(vit);
2043 const scalar cellSizeSqr =
sqr(cellSize);
2045 forAll(surfaceIntersections, sI)
2047 pointIndexHitAndFeature surfHitI = surfaceIntersections[sI];
2049 bool keepSurfacePoint =
true;
2051 if (!surfHitI.first().hit())
2056 const Foam::point& surfPt = surfHitI.first().hitPoint();
2058 bool isNearFeaturePt = nearFeaturePt(surfPt);
2060 bool isNearFeatureEdge = surfacePtNearFeatureEdge(surfPt);
2062 bool isNearSurfacePt = nearSurfacePoint(surfHitI);
2064 if (isNearFeaturePt || isNearSurfacePt || isNearFeatureEdge)
2066 keepSurfacePoint =
false;
2073 const scalar searchRadiusSqr = edgeSearchDistCoeffSqr*cellSizeSqr;
2075 geometryToConformTo_.findAllNearestEdges
2083 forAll(edHitsByFeature, i)
2085 const label featureHit = featuresHit[i];
2100 && !decomposition().positionOnThisProcessor(edPt)
2107 if (!nearFeaturePt(edPt))
2112 < surfacePtReplaceDistCoeffSqr*cellSizeSqr
2122 keepSurfacePoint =
false;
2134 !nearFeatureEdgeLocation(edHit, nearestEdgeHit)
2137 appendToEdgeLocationTree(edPt);
2139 edgeToTreeShape.append
2141 existingEdgeLocations_.size() - 1
2146 featureEdgeHits.append
2148 pointIndexHitAndFeature(edHit, featureHit)
2154 surfacePtToEdgePtDist.insert
2156 existingEdgeLocations_.size() - 1,
2162 label hitIndex = nearestEdgeHit.index();
2169 < surfacePtToEdgePtDist[hitIndex]
2172 featureEdgeHits[hitIndex] =
2173 pointIndexHitAndFeature(edHit, featureHit);
2175 existingEdgeLocations_[hitIndex] =
2177 surfacePtToEdgePtDist[hitIndex] =
2183 edgeLocationTreePtr_().remove(hitIndex);
2184 edgeLocationTreePtr_().insert
2196 if (keepSurfacePoint)
2198 surfaceHits.append(surfHitI);
2199 appendToSurfacePtTree(surfPt);
2200 surfaceToTreeShape.append(existingSurfacePtLocations_.size() - 1);
2214 Info<<
nl <<
"Storing surface conformation" <<
endl;
2216 surfaceConformationVertices_.clear();
2219 DynamicList<Vb> tempSurfaceVertices(number_of_vertices()/10);
2223 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
2224 vit != finite_vertices_end();
2233 && vit->boundaryPoint()
2234 && !vit->featurePoint()
2235 && !vit->constrained()
2238 tempSurfaceVertices.append
2251 tempSurfaceVertices.shrink();
2253 surfaceConformationVertices_.transfer(tempSurfaceVertices);
2258 label(surfaceConformationVertices_.size()),
2261 <<
" vertices" <<
nl <<
endl;
2267 Info<<
nl <<
"Reinserting stored surface conformation" <<
endl;
2269 Map<label> oldToNewIndices =
2270 insertPointPairs(surfaceConformationVertices_,
true,
true);
2272 ptPairs_.reIndex(oldToNewIndices);
2274 PackedBoolList selectedElems(surfaceConformationVertices_.size(),
true);
2276 forAll(surfaceConformationVertices_, vI)
2278 Vb& v = surfaceConformationVertices_[vI];
2283 if (iter != oldToNewIndices.end())
2285 const label newIndex = iter();
2293 selectedElems[vI] =
false;
2298 inplaceSubset<PackedBoolList, List<Vb> >
2301 surfaceConformationVertices_