50 externalDisplacementMeshMover,
61 bool Foam::medialAxisMeshMover::isMaxEdge
63 const List<PointData<vector>>& pointWallDist,
66 const bool disableWallEdges
76 scalar magV0(
mag(v0));
83 scalar magV1(
mag(v1));
117 if ((pointWallDist[
e[0]].data() & pointWallDist[
e[1]].data()) < minCos)
126 void Foam::medialAxisMeshMover::update(
const dictionary& coeffDict)
129 <<
" : Calculating distance to Medial Axis ..." <<
endl;
134 const labelList& meshPoints = pp.meshPoints();
141 const label nSmoothSurfaceNormals
143 meshRefinement::get<label>
146 "nSmoothSurfaceNormals",
154 scalar minMedialAxisAngle(
Zero);
157 !coeffDict.readCompat
159 "minMedialAxisAngle",
160 {{
"minMedianAxisAngle", 1712 }},
168 <<
"Entry '" <<
"minMedialAxisAngle"
169 <<
"' not found in dictionary " << coeffDict.name() <<
endl;
175 const scalar featureAngle
177 meshRefinement::get<scalar>(coeffDict,
"featureAngle", dryRun_)
181 const scalar slipFeatureAngle =
183 coeffDict.found(
"slipFeatureAngle")
184 ? coeffDict.get<scalar>(
"slipFeatureAngle")
189 const label nSmoothNormals
191 meshRefinement::get<label>(coeffDict,
"nSmoothNormals", dryRun_)
195 const label nMedialAxisIter = coeffDict.lookupOrDefault<
label>
201 const bool disableWallEdges = coeffDict.lookupOrDefault<
bool>
226 const bitSet isPatchMasterPoint
234 const bitSet isPatchMasterEdge
249 fieldSmoother_.smoothPatchNormals
251 nSmoothSurfaceNormals,
263 List<PointData<vector>> pointWallDist(
mesh().
nPoints());
266 int dummyTrackData = 0;
272 List<PointData<vector>> wallInfo(meshPoints.size());
274 forAll(meshPoints, patchPointI)
276 label pointI = meshPoints[patchPointI];
277 wallInfo[patchPointI] = PointData<vector>
281 pointNormals[patchPointI]
286 List<PointData<vector>> edgeWallDist(
mesh().nEdges());
287 PointEdgeWave<PointData<vector>> wallDistCalc
297 wallDistCalc.iterate(nMedialAxisIter);
301 wallDistCalc.nUnvisitedPoints(),
307 if (nMedialAxisIter > 0)
310 <<
" : Limited walk to " << nMedialAxisIter
311 <<
" steps. Not visited " << nUnvisit
313 <<
" points" <<
endl;
318 <<
"Walking did not visit all points." <<
nl
319 <<
" Did not visit " << nUnvisit
321 <<
" points. This is not necessarily a problem" <<
nl
322 <<
" and might be due to faceZones splitting of part"
323 <<
" of the domain." <<
nl <<
endl;
332 List<pointEdgePoint> pointMedialDist(
mesh().
nPoints());
333 List<pointEdgePoint> edgeMedialDist(
mesh().nEdges());
336 DynamicList<pointEdgePoint> maxInfo(meshPoints.size());
337 DynamicList<label> maxPoints(meshPoints.size());
345 const edge&
e = edges[edgeI];
349 !pointWallDist[
e[0]].valid(dummyTrackData)
350 || !pointWallDist[
e[1]].valid(dummyTrackData)
358 if (!pointMedialDist[pointI].valid(dummyTrackData))
360 maxPoints.append(pointI);
369 pointMedialDist[pointI] = maxInfo.last();
380 minMedialAxisAngleCos,
391 scalar eMag =
mag(eVec);
398 const point& origin0 = pointWallDist[
e[0]].origin();
400 const point& origin1 = pointWallDist[
e[1]].origin();
401 scalar dist0 = (
p0-origin0) & eVec;
402 scalar dist1 = (origin1-p1) & eVec;
403 scalar
s = 0.5*(dist1+eMag+dist0);
415 else if (
s >= dist0+eMag)
426 medialAxisPt =
p0+(
s-dist0)*eVec;
435 if (!pointMedialDist[pointI].valid(dummyTrackData))
437 maxPoints.append(pointI);
446 pointMedialDist[pointI] = maxInfo.last();
463 const polyPatch& pp =
patches[patchI];
470 && !isA<emptyPolyPatch>(pp)
471 && !adaptPatches.found(patchI)
474 const labelList& meshPoints = pp.meshPoints();
481 if (pvf.fixesValue())
485 <<
" : Inserting all points on patch " << pp.name()
490 label pointI = meshPoints[i];
491 if (!pointMedialDist[pointI].valid(dummyTrackData))
493 maxPoints.append(pointI);
502 pointMedialDist[pointI] = maxInfo.last();
514 <<
" : Inserting points on patch " << pp.name()
515 <<
" if angle to nearest layer patch > "
516 << slipFeatureAngle <<
" degrees." <<
endl;
529 label pointI = meshPoints[i];
533 pointWallDist[pointI].valid(dummyTrackData)
534 && !pointMedialDist[pointI].valid(dummyTrackData)
540 -pointWallDist[pointI].data()
543 if (cosAngle > slipFeatureAngleCos)
548 maxPoints.append(pointI);
557 pointMedialDist[pointI] = maxInfo.last();
577 OBJstream str(
mesh().time().timePath()/
"medialSurfacePoints.obj");
582 "medialSurfacePoints",
587 <<
" : Writing estimated medial surface:" <<
nl <<
incrIndent
592 for (
const auto& info : maxInfo)
594 str.
write(info.origin());
601 PointEdgeWave<pointEdgePoint> medialDistCalc
612 medialDistCalc.iterate(2*nMedialAxisIter);
616 forAll(pointMedialDist, pointI)
618 if (pointMedialDist[pointI].valid(dummyTrackData))
622 pointMedialDist[pointI].distSqr()
624 medialVec_[pointI] = pointMedialDist[pointI].origin();
629 medialDist_[pointI] = 0.0;
630 medialVec_[pointI] =
point(1, 0, 0);
638 if (!pointWallDist[i].valid(dummyTrackData))
640 dispVec_[i] =
vector(1, 0, 0);
644 dispVec_[i] = pointWallDist[i].data();
649 fieldSmoother_.smoothNormals
659 forAll(medialRatio_, pointI)
661 if (!pointWallDist[pointI].valid(dummyTrackData))
663 medialRatio_[pointI] = 0.0;
667 scalar wDist2 = pointWallDist[pointI].distSqr();
668 scalar mDist = medialDist_[pointI];
670 if (wDist2 <
sqr(SMALL) && mDist < SMALL)
677 medialRatio_[pointI] = 0.0;
681 medialRatio_[pointI] = mDist / (
Foam::sqrt(wDist2) + mDist);
691 <<
indent <<
"ratio of medial distance to wall distance : "
692 << medialRatio_.
name() <<
nl
693 <<
indent <<
"distance to nearest medial axis : "
694 << medialDist_.
name() <<
nl
695 <<
indent <<
"nearest medial axis location : "
696 << medialVec_.
name() <<
nl
697 <<
indent <<
"normal at nearest wall : "
702 medialRatio_.write();
709 bool Foam::medialAxisMeshMover::unmarkExtrusion
711 const label patchPointI,
713 List<snappyLayerDriver::extrudeMode>& extrudeStatus
719 patchDisp[patchPointI] =
Zero;
725 patchDisp[patchPointI] =
Zero;
733 void Foam::medialAxisMeshMover::syncPatchDisplacement
737 List<snappyLayerDriver::extrudeMode>& extrudeStatus
741 const labelList& meshPoints = pp.meshPoints();
743 label nChangedTotal = 0;
755 minMagSqrEqOp<vector>(),
762 if (
mag(patchDisp[i]) < minThickness[i])
764 if (unmarkExtrusion(i, patchDisp, extrudeStatus))
835 nChangedTotal += nChanged;
851 void Foam::medialAxisMeshMover::
852 handleFeatureAngleLayerTerminations
855 const bitSet& isPatchMasterPoint,
857 List<snappyLayerDriver::extrudeMode>& extrudeStatus,
867 boolList extrudedFaces(pp.size(),
true);
869 forAll(pp.localFaces(), faceI)
871 const face&
f = pp.localFaces()[faceI];
877 extrudedFaces[faceI] =
false;
892 List<List<point>> edgeFaceNormals(pp.nEdges());
893 List<List<bool>> edgeFaceExtrude(pp.nEdges());
900 const labelList& eFaces = edgeFaces[edgeI];
902 edgeFaceNormals[edgeI].
setSize(eFaces.size());
903 edgeFaceExtrude[edgeI].setSize(eFaces.size());
906 label faceI = eFaces[i];
908 edgeFaceExtrude[edgeI][i] = extrudedFaces[faceI];
917 globalMeshData::ListPlusEqOp<List<point>>(),
926 globalMeshData::ListPlusEqOp<List<bool>>(),
931 forAll(edgeFaceNormals, edgeI)
933 const List<point>& eFaceNormals = edgeFaceNormals[edgeI];
934 const List<bool>& eFaceExtrude = edgeFaceExtrude[edgeI];
936 if (eFaceNormals.size() == 2)
938 const edge&
e = pp.edges()[edgeI];
948 if (!eFaceExtrude[0] || !eFaceExtrude[1])
950 const vector& n0 = eFaceNormals[0];
951 const vector& n1 = eFaceNormals[1];
953 if ((n0 & n1) < minCos)
955 if (unmarkExtrusion(v0, patchDisp, extrudeStatus))
957 if (isPatchMasterPoint[v0])
962 if (unmarkExtrusion(v1, patchDisp, extrudeStatus))
964 if (isPatchMasterPoint[v1])
984 void Foam::medialAxisMeshMover::findIsolatedRegions
986 const scalar minCosLayerTermination,
987 const bool detectExtrusionIsland,
988 const bitSet& isPatchMasterPoint,
989 const bitSet& isPatchMasterEdge,
992 List<snappyLayerDriver::extrudeMode>& extrudeStatus,
998 const labelList& meshPoints = pp.meshPoints();
1001 Info<< typeName <<
" : Removing isolated regions ..." <<
nl
1002 <<
indent <<
"- if partially extruded faces make angle < "
1004 if (detectExtrusionIsland)
1006 Info<<
indent <<
"- if exclusively surrounded by non-extruded points"
1011 Info<<
indent <<
"- if exclusively surrounded by non-extruded faces"
1016 label nPointCounter = 0;
1022 if (minCosLayerTermination > -1)
1024 handleFeatureAngleLayerTerminations
1026 minCosLayerTermination,
1035 syncPatchDisplacement(minThickness, patchDisp, extrudeStatus);
1046 boolList keptPoints(pp.nPoints(),
false);
1048 if (detectExtrusionIsland)
1057 const face&
f = pp.localFaces()[faceI];
1063 if (islandPoint[faceI] == -1)
1066 islandPoint[faceI] =
f[fp];
1068 else if (islandPoint[faceI] != -2)
1071 islandPoint[faceI] = -2;
1083 forAll(pointFaces, patchPointI)
1092 if (islandPoint[faceI] != patchPointI)
1094 keptPoints[patchPointI] =
true;
1107 boolList extrudedFaces(pp.size(),
true);
1108 forAll(pp.localFaces(), faceI)
1110 const face&
f = pp.localFaces()[faceI];
1115 extrudedFaces[faceI] =
false;
1123 forAll(keptPoints, patchPointI)
1130 if (extrudedFaces[faceI])
1132 keptPoints[patchPointI] =
true;
1151 forAll(keptPoints, patchPointI)
1153 if (!keptPoints[patchPointI])
1155 if (unmarkExtrusion(patchPointI, patchDisp, extrudeStatus))
1170 const edgeList& edges = pp.edges();
1180 if (isPatchMasterEdge[edgeI])
1182 const edge&
e = edges[edgeI];
1189 isolatedPoint[v0] += 1;
1193 isolatedPoint[v1] += 1;
1211 const face&
f = pp.localFaces()[faceI];
1212 bool failed =
false;
1215 if (isolatedPoint[
f[fp]] > 2)
1221 bool allPointsExtruded =
true;
1228 allPointsExtruded =
false;
1233 if (allPointsExtruded)
1254 reduce(nPointCounter, sumOp<label>());
1256 <<
" : Number of isolated points extrusion stopped : "<< nPointCounter
1263 Foam::medialAxisMeshMover::medialAxisMeshMover
1272 adaptPatchIDs_(getFixedValueBCs(pointDisplacement)),
1273 adaptPatchPtr_(getPatch(
mesh(), adaptPatchIDs_)),
1279 pointDisplacement.time().timeName(),
1280 pointDisplacement.db(),
1300 fieldSmoother_(
mesh()),
1306 pointDisplacement.time().timeName(),
1307 pointDisplacement.db(),
1320 pointDisplacement.time().timeName(),
1321 pointDisplacement.db(),
1334 pointDisplacement.time().timeName(),
1335 pointDisplacement.db(),
1348 pointDisplacement.time().timeName(),
1349 pointDisplacement.db(),
1370 void Foam::medialAxisMeshMover::calculateDisplacement
1378 Info<< typeName <<
" : Smoothing using Medial Axis ..." <<
endl;
1388 const label nSmoothDisplacement =
1392 const scalar maxThicknessToMedialRatio =
1393 coeffDict.
get<scalar>(
"maxThicknessToMedialRatio");
1396 const scalar featureAngle = coeffDict.
get<scalar>(
"featureAngle");
1401 "layerTerminationAngle",
1407 const label nSmoothPatchThickness =
1408 coeffDict.
get<
label>(
"nSmoothThickness");
1420 "detectExtrusionIsland",
1439 const bitSet isPatchMasterPoint
1447 const bitSet isPatchMasterEdge
1459 forAll(thickness, patchPointI)
1463 thickness[patchPointI] = 0.0;
1467 label numThicknessRatioExclude = 0;
1472 autoPtr<OBJstream> str;
1480 /
"thicknessRatioExcludePoints_"
1486 <<
" : Writing points with too large an extrusion distance to "
1487 << str().name() <<
endl;
1490 autoPtr<OBJstream> medialVecStr;
1498 /
"thicknessRatioExcludeMedialVec_"
1504 <<
" : Writing medial axis vectors on points with too large"
1505 <<
" an extrusion distance to " << medialVecStr().name() <<
endl;
1508 forAll(meshPoints, patchPointI)
1512 label pointI = meshPoints[patchPointI];
1516 scalar mDist = medialDist_[pointI];
1517 scalar thicknessRatio = thickness[patchPointI]/(mDist+VSMALL);
1528 thicknessRatio *= (
n & mVec);
1530 if (thicknessRatio > maxThicknessToMedialRatio)
1535 Pout<<
"truncating displacement at "
1537 <<
" from " << thickness[patchPointI]
1541 minThickness[patchPointI]
1542 +thickness[patchPointI]
1544 <<
" medial direction:" << mVec
1545 <<
" extrusion direction:" <<
n
1546 <<
" with thicknessRatio:" << thicknessRatio
1550 thickness[patchPointI] =
1551 0.5*(minThickness[patchPointI]+thickness[patchPointI]);
1553 patchDisp[patchPointI] = thickness[patchPointI]*
n;
1555 if (isPatchMasterPoint[patchPointI])
1557 numThicknessRatioExclude++;
1563 str().write(
linePointRef(pt, pt+patchDisp[patchPointI]));
1565 if (medialVecStr.valid())
1568 medialVecStr().write
1581 reduce(numThicknessRatioExclude, sumOp<label>());
1582 Info<< typeName <<
" : Reducing layer thickness at "
1583 << numThicknessRatioExclude
1584 <<
" nodes where thickness to medial axis distance is large " <<
endl;
1591 minCosLayerTermination,
1592 detectExtrusionIsland,
1604 forAll(thickness, patchPointI)
1608 thickness[patchPointI] = 0.0;
1617 fieldSmoother_.minSmoothField
1619 nSmoothPatchThickness,
1629 int dummyTrackData = 0;
1631 List<PointData<scalar>> pointWallDist(
mesh().
nPoints());
1638 List<PointData<scalar>> edgeWallDist(
mesh().nEdges());
1639 labelList wallPoints(meshPoints.size());
1642 List<PointData<scalar>> wallInfo(meshPoints.size());
1644 forAll(meshPoints, patchPointI)
1646 label pointI = meshPoints[patchPointI];
1647 wallPoints[patchPointI] = pointI;
1648 wallInfo[patchPointI] = PointData<scalar>
1652 thickness[patchPointI]
1657 PointEdgeWave<PointData<scalar>> wallDistCalc
1667 wallDistCalc.iterate(nMedialAxisIter);
1672 pointField& displacement = pointDisplacement_;
1674 forAll(displacement, pointI)
1676 if (!pointWallDist[pointI].valid(dummyTrackData))
1678 displacement[pointI] =
Zero;
1687 displacement[pointI] =
1688 -medialRatio_[pointI]
1689 *pointWallDist[pointI].data()
1696 if (nSmoothDisplacement > 0)
1698 bitSet isToBeSmoothed(displacement.size(),
false);
1702 if (medialRatio_[i] > SMALL && medialRatio_[i] < 1-SMALL)
1704 isToBeSmoothed.set(i);
1708 fieldSmoother_.smoothLambdaMuDisplacement
1710 nSmoothDisplacement,
1720 bool Foam::medialAxisMeshMover::shrinkMesh
1722 const dictionary& meshQualityDict,
1723 const label nAllowableErrors,
1728 const label nSnap = meshQualityDict.get<
label>(
"nRelaxIter");
1733 meshMover_.setDisplacementPatchFields();
1735 Info<< typeName <<
" : Moving mesh ..." <<
endl;
1736 scalar oldErrorReduction = -1;
1738 bool meshOk =
false;
1740 for (
label iter = 0; iter < 2*nSnap ; iter++)
1743 <<
" : Iteration " << iter <<
endl;
1747 <<
" : Displacement scaling for error reduction set to 0."
1749 oldErrorReduction = meshMover_.setErrorReduction(0.0);
1754 meshMover_.scaleMesh
1758 meshMover_.paramDict(),
1765 Info<< typeName <<
" : Successfully moved mesh" <<
endl;
1771 if (oldErrorReduction >= 0)
1773 meshMover_.setErrorReduction(oldErrorReduction);
1776 Info<< typeName <<
" : Finished moving mesh ..." <<
endl;
1785 const label nAllowableErrors,
1793 const word minThicknessName = moveDict.
get<
word>(
"minThicknessName");
1800 if (minThicknessName ==
"none")
1806 (minThicknessName ==
"none")
1814 pointDisplacement_.internalField(),
1823 forAll(extrudeStatus, pointI)
1825 if (
mag(patchDisp[pointI]) <= minThickness[pointI]+SMALL)
1833 calculateDisplacement(moveDict, minThickness, extrudeStatus, patchDisp);
1850 meshMover_.movePoints();
1853 meshMover_.correct();