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.getOrDefault<scalar>(
"slipFeatureAngle", (0.5*featureAngle))
187 const label nSmoothNormals
189 meshRefinement::get<label>(coeffDict,
"nSmoothNormals", dryRun_)
193 const label nMedialAxisIter = coeffDict.getOrDefault<label>
199 const bool disableWallEdges = coeffDict.getOrDefault<
bool>
224 const bitSet isPatchMasterPoint
232 const bitSet isPatchMasterEdge
247 fieldSmoother_.smoothPatchNormals
249 nSmoothSurfaceNormals,
261 List<PointData<vector>> pointWallDist(
mesh().
nPoints());
264 int dummyTrackData = 0;
270 List<PointData<vector>> wallInfo(meshPoints.size());
272 forAll(meshPoints, patchPointI)
274 label pointI = meshPoints[patchPointI];
275 wallInfo[patchPointI] = PointData<vector>
279 pointNormals[patchPointI]
284 List<PointData<vector>> edgeWallDist(
mesh().nEdges());
285 PointEdgeWave<PointData<vector>> wallDistCalc
295 wallDistCalc.iterate(nMedialAxisIter);
299 wallDistCalc.nUnvisitedPoints(),
305 if (nMedialAxisIter > 0)
308 <<
" : Limited walk to " << nMedialAxisIter
309 <<
" steps. Not visited " << nUnvisit
311 <<
" points" <<
endl;
316 <<
"Walking did not visit all points." <<
nl
317 <<
" Did not visit " << nUnvisit
319 <<
" points. This is not necessarily a problem" <<
nl
320 <<
" and might be due to faceZones splitting of part"
321 <<
" of the domain." <<
nl <<
endl;
330 List<pointEdgePoint> pointMedialDist(
mesh().
nPoints());
331 List<pointEdgePoint> edgeMedialDist(
mesh().nEdges());
334 DynamicList<pointEdgePoint> maxInfo(meshPoints.size());
335 DynamicList<label> maxPoints(meshPoints.size());
343 const edge&
e = edges[edgeI];
347 !pointWallDist[
e[0]].valid(dummyTrackData)
348 || !pointWallDist[
e[1]].valid(dummyTrackData)
354 label pointI =
e[ep];
356 if (!pointMedialDist[pointI].valid(dummyTrackData))
358 maxPoints.append(pointI);
367 pointMedialDist[pointI] = maxInfo.last();
378 minMedialAxisAngleCos,
389 scalar eMag =
mag(eVec);
396 const point& origin0 = pointWallDist[
e[0]].origin();
398 const point& origin1 = pointWallDist[
e[1]].origin();
399 scalar dist0 = (
p0-origin0) & eVec;
400 scalar dist1 = (origin1-p1) & eVec;
401 scalar
s = 0.5*(dist1+eMag+dist0);
413 else if (
s >= dist0+eMag)
424 medialAxisPt =
p0+(
s-dist0)*eVec;
431 label pointI =
e[ep];
433 if (!pointMedialDist[pointI].valid(dummyTrackData))
435 maxPoints.append(pointI);
444 pointMedialDist[pointI] = maxInfo.last();
461 const polyPatch& pp =
patches[patchI];
468 && !isA<emptyPolyPatch>(pp)
469 && !adaptPatches.found(patchI)
472 const labelList& meshPoints = pp.meshPoints();
479 if (pvf.fixesValue())
483 <<
" : Inserting all points on patch " << pp.name()
488 label pointI = meshPoints[i];
489 if (!pointMedialDist[pointI].valid(dummyTrackData))
491 maxPoints.append(pointI);
500 pointMedialDist[pointI] = maxInfo.last();
512 <<
" : Inserting points on patch " << pp.name()
513 <<
" if angle to nearest layer patch > "
514 << slipFeatureAngle <<
" degrees." <<
endl;
527 label pointI = meshPoints[i];
531 pointWallDist[pointI].valid(dummyTrackData)
532 && !pointMedialDist[pointI].valid(dummyTrackData)
538 -pointWallDist[pointI].data()
541 if (cosAngle > slipFeatureAngleCos)
546 maxPoints.append(pointI);
555 pointMedialDist[pointI] = maxInfo.last();
575 OBJstream str(
mesh().time().timePath()/
"medialSurfacePoints.obj");
580 "medialSurfacePoints",
585 <<
" : Writing estimated medial surface:" <<
nl <<
incrIndent
590 for (
const auto& info : maxInfo)
592 str.
write(info.origin());
599 PointEdgeWave<pointEdgePoint> medialDistCalc
610 medialDistCalc.iterate(2*nMedialAxisIter);
614 forAll(pointMedialDist, pointI)
616 if (pointMedialDist[pointI].valid(dummyTrackData))
620 pointMedialDist[pointI].distSqr()
622 medialVec_[pointI] = pointMedialDist[pointI].origin();
627 medialDist_[pointI] = 0.0;
628 medialVec_[pointI] =
point(1, 0, 0);
636 if (!pointWallDist[i].valid(dummyTrackData))
638 dispVec_[i] =
vector(1, 0, 0);
642 dispVec_[i] = pointWallDist[i].data();
647 fieldSmoother_.smoothNormals
657 forAll(medialRatio_, pointI)
659 if (!pointWallDist[pointI].valid(dummyTrackData))
661 medialRatio_[pointI] = 0.0;
665 scalar wDist2 = pointWallDist[pointI].distSqr();
666 scalar mDist = medialDist_[pointI];
668 if (wDist2 <
sqr(SMALL) && mDist < SMALL)
675 medialRatio_[pointI] = 0.0;
679 medialRatio_[pointI] = mDist / (
Foam::sqrt(wDist2) + mDist);
689 <<
indent <<
"ratio of medial distance to wall distance : "
690 << medialRatio_.
name() <<
nl
691 <<
indent <<
"distance to nearest medial axis : "
692 << medialDist_.
name() <<
nl
693 <<
indent <<
"nearest medial axis location : "
694 << medialVec_.
name() <<
nl
695 <<
indent <<
"normal at nearest wall : "
700 medialRatio_.write();
707 bool Foam::medialAxisMeshMover::unmarkExtrusion
709 const label patchPointI,
711 List<snappyLayerDriver::extrudeMode>& extrudeStatus
717 patchDisp[patchPointI] =
Zero;
723 patchDisp[patchPointI] =
Zero;
731 void Foam::medialAxisMeshMover::syncPatchDisplacement
735 List<snappyLayerDriver::extrudeMode>& extrudeStatus
739 const labelList& meshPoints = pp.meshPoints();
741 label nChangedTotal = 0;
753 minMagSqrEqOp<vector>(),
760 if (
mag(patchDisp[i]) < minThickness[i])
762 if (unmarkExtrusion(i, patchDisp, extrudeStatus))
833 nChangedTotal += nChanged;
849 void Foam::medialAxisMeshMover::
850 handleFeatureAngleLayerTerminations
853 const bitSet& isPatchMasterPoint,
855 List<snappyLayerDriver::extrudeMode>& extrudeStatus,
865 boolList extrudedFaces(pp.size(),
true);
867 forAll(pp.localFaces(), faceI)
869 const face&
f = pp.localFaces()[faceI];
875 extrudedFaces[faceI] =
false;
890 List<List<point>> edgeFaceNormals(pp.nEdges());
891 List<List<bool>> edgeFaceExtrude(pp.nEdges());
898 const labelList& eFaces = edgeFaces[edgeI];
900 edgeFaceNormals[edgeI].
setSize(eFaces.size());
901 edgeFaceExtrude[edgeI].setSize(eFaces.size());
904 label faceI = eFaces[i];
906 edgeFaceExtrude[edgeI][i] = extrudedFaces[faceI];
915 ListOps::appendEqOp<point>(),
924 ListOps::appendEqOp<bool>(),
929 forAll(edgeFaceNormals, edgeI)
931 const List<point>& eFaceNormals = edgeFaceNormals[edgeI];
932 const List<bool>& eFaceExtrude = edgeFaceExtrude[edgeI];
934 if (eFaceNormals.size() == 2)
936 const edge&
e = pp.edges()[edgeI];
946 if (!eFaceExtrude[0] || !eFaceExtrude[1])
948 const vector& n0 = eFaceNormals[0];
949 const vector& n1 = eFaceNormals[1];
951 if ((n0 & n1) < minCos)
953 if (unmarkExtrusion(v0, patchDisp, extrudeStatus))
955 if (isPatchMasterPoint[v0])
960 if (unmarkExtrusion(v1, patchDisp, extrudeStatus))
962 if (isPatchMasterPoint[v1])
982 void Foam::medialAxisMeshMover::findIsolatedRegions
984 const scalar minCosLayerTermination,
985 const bool detectExtrusionIsland,
986 const bitSet& isPatchMasterPoint,
987 const bitSet& isPatchMasterEdge,
990 List<snappyLayerDriver::extrudeMode>& extrudeStatus,
996 const labelList& meshPoints = pp.meshPoints();
999 Info<< typeName <<
" : Removing isolated regions ..." <<
nl
1000 <<
indent <<
"- if partially extruded faces make angle < "
1002 if (detectExtrusionIsland)
1004 Info<<
indent <<
"- if exclusively surrounded by non-extruded points"
1009 Info<<
indent <<
"- if exclusively surrounded by non-extruded faces"
1014 label nPointCounter = 0;
1020 if (minCosLayerTermination > -1)
1022 handleFeatureAngleLayerTerminations
1024 minCosLayerTermination,
1033 syncPatchDisplacement(minThickness, patchDisp, extrudeStatus);
1044 boolList keptPoints(pp.nPoints(),
false);
1046 if (detectExtrusionIsland)
1055 const face&
f = pp.localFaces()[faceI];
1061 if (islandPoint[faceI] == -1)
1064 islandPoint[faceI] =
f[fp];
1066 else if (islandPoint[faceI] != -2)
1069 islandPoint[faceI] = -2;
1081 forAll(pointFaces, patchPointI)
1090 if (islandPoint[faceI] != patchPointI)
1092 keptPoints[patchPointI] =
true;
1105 boolList extrudedFaces(pp.size(),
true);
1106 forAll(pp.localFaces(), faceI)
1108 const face&
f = pp.localFaces()[faceI];
1113 extrudedFaces[faceI] =
false;
1121 forAll(keptPoints, patchPointI)
1128 if (extrudedFaces[faceI])
1130 keptPoints[patchPointI] =
true;
1149 forAll(keptPoints, patchPointI)
1151 if (!keptPoints[patchPointI])
1153 if (unmarkExtrusion(patchPointI, patchDisp, extrudeStatus))
1168 const edgeList& edges = pp.edges();
1178 if (isPatchMasterEdge[edgeI])
1180 const edge&
e = edges[edgeI];
1187 isolatedPoint[v0] += 1;
1191 isolatedPoint[v1] += 1;
1209 const face&
f = pp.localFaces()[faceI];
1210 bool failed =
false;
1213 if (isolatedPoint[
f[fp]] > 2)
1219 bool allPointsExtruded =
true;
1226 allPointsExtruded =
false;
1231 if (allPointsExtruded)
1252 reduce(nPointCounter, sumOp<label>());
1254 <<
" : Number of isolated points extrusion stopped : "<< nPointCounter
1261 Foam::medialAxisMeshMover::medialAxisMeshMover
1270 adaptPatchIDs_(getFixedValueBCs(pointDisplacement)),
1271 adaptPatchPtr_(getPatch(
mesh(), adaptPatchIDs_)),
1277 pointDisplacement.time().timeName(),
1278 pointDisplacement.db(),
1298 fieldSmoother_(
mesh()),
1304 pointDisplacement.time().timeName(),
1305 pointDisplacement.db(),
1318 pointDisplacement.time().timeName(),
1319 pointDisplacement.db(),
1332 pointDisplacement.time().timeName(),
1333 pointDisplacement.db(),
1346 pointDisplacement.time().timeName(),
1347 pointDisplacement.db(),
1368 void Foam::medialAxisMeshMover::calculateDisplacement
1376 Info<< typeName <<
" : Smoothing using Medial Axis ..." <<
endl;
1386 const label nSmoothDisplacement =
1390 const scalar maxThicknessToMedialRatio =
1391 coeffDict.
get<scalar>(
"maxThicknessToMedialRatio");
1394 const scalar featureAngle = coeffDict.
get<scalar>(
"featureAngle");
1397 scalar layerTerminationAngle = coeffDict.
getOrDefault<scalar>
1399 "layerTerminationAngle",
1405 const label nSmoothPatchThickness =
1406 coeffDict.
get<label>(
"nSmoothThickness");
1409 const label nMedialAxisIter = coeffDict.
getOrDefault<label>
1416 const bool detectExtrusionIsland = coeffDict.
getOrDefault
1418 "detectExtrusionIsland",
1437 const bitSet isPatchMasterPoint
1445 const bitSet isPatchMasterEdge
1457 forAll(thickness, patchPointI)
1461 thickness[patchPointI] = 0.0;
1465 label numThicknessRatioExclude = 0;
1470 autoPtr<OBJstream> str;
1478 /
"thicknessRatioExcludePoints_"
1484 <<
" : Writing points with too large an extrusion distance to "
1485 << str().name() <<
endl;
1488 autoPtr<OBJstream> medialVecStr;
1496 /
"thicknessRatioExcludeMedialVec_"
1502 <<
" : Writing medial axis vectors on points with too large"
1503 <<
" an extrusion distance to " << medialVecStr().name() <<
endl;
1506 forAll(meshPoints, patchPointI)
1510 label pointI = meshPoints[patchPointI];
1514 scalar mDist = medialDist_[pointI];
1515 scalar thicknessRatio = thickness[patchPointI]/(mDist+VSMALL);
1526 thicknessRatio *= (
n & mVec);
1528 if (thicknessRatio > maxThicknessToMedialRatio)
1533 Pout<<
"truncating displacement at "
1535 <<
" from " << thickness[patchPointI]
1539 minThickness[patchPointI]
1540 +thickness[patchPointI]
1542 <<
" medial direction:" << mVec
1543 <<
" extrusion direction:" <<
n
1544 <<
" with thicknessRatio:" << thicknessRatio
1548 thickness[patchPointI] =
1549 0.5*(minThickness[patchPointI]+thickness[patchPointI]);
1551 patchDisp[patchPointI] = thickness[patchPointI]*
n;
1553 if (isPatchMasterPoint[patchPointI])
1555 numThicknessRatioExclude++;
1561 str().write(
linePointRef(pt, pt+patchDisp[patchPointI]));
1566 medialVecStr().write
1579 reduce(numThicknessRatioExclude, sumOp<label>());
1580 Info<< typeName <<
" : Reducing layer thickness at "
1581 << numThicknessRatioExclude
1582 <<
" nodes where thickness to medial axis distance is large " <<
endl;
1589 minCosLayerTermination,
1590 detectExtrusionIsland,
1602 forAll(thickness, patchPointI)
1606 thickness[patchPointI] = 0.0;
1615 fieldSmoother_.minSmoothField
1617 nSmoothPatchThickness,
1627 int dummyTrackData = 0;
1629 List<PointData<scalar>> pointWallDist(
mesh().
nPoints());
1636 List<PointData<scalar>> edgeWallDist(
mesh().nEdges());
1637 labelList wallPoints(meshPoints.size());
1640 List<PointData<scalar>> wallInfo(meshPoints.size());
1642 forAll(meshPoints, patchPointI)
1644 label pointI = meshPoints[patchPointI];
1645 wallPoints[patchPointI] = pointI;
1646 wallInfo[patchPointI] = PointData<scalar>
1650 thickness[patchPointI]
1655 PointEdgeWave<PointData<scalar>> wallDistCalc
1665 wallDistCalc.iterate(nMedialAxisIter);
1670 pointField& displacement = pointDisplacement_;
1672 forAll(displacement, pointI)
1674 if (!pointWallDist[pointI].valid(dummyTrackData))
1676 displacement[pointI] =
Zero;
1685 displacement[pointI] =
1686 -medialRatio_[pointI]
1687 *pointWallDist[pointI].data()
1694 if (nSmoothDisplacement > 0)
1696 bitSet isToBeSmoothed(displacement.size(),
false);
1700 if (medialRatio_[i] > SMALL && medialRatio_[i] < 1-SMALL)
1702 isToBeSmoothed.set(i);
1706 fieldSmoother_.smoothLambdaMuDisplacement
1708 nSmoothDisplacement,
1718 bool Foam::medialAxisMeshMover::shrinkMesh
1720 const dictionary& meshQualityDict,
1721 const label nAllowableErrors,
1726 const label nSnap = meshQualityDict.get<label>(
"nRelaxIter");
1731 meshMover_.setDisplacementPatchFields();
1733 Info<< typeName <<
" : Moving mesh ..." <<
endl;
1734 scalar oldErrorReduction = -1;
1736 bool meshOk =
false;
1738 for (label iter = 0; iter < 2*nSnap ; iter++)
1741 <<
" : Iteration " << iter <<
endl;
1745 <<
" : Displacement scaling for error reduction set to 0."
1747 oldErrorReduction = meshMover_.setErrorReduction(0.0);
1752 meshMover_.scaleMesh
1756 meshMover_.paramDict(),
1763 Info<< typeName <<
" : Successfully moved mesh" <<
endl;
1769 if (oldErrorReduction >= 0)
1771 meshMover_.setErrorReduction(oldErrorReduction);
1774 Info<< typeName <<
" : Finished moving mesh ..." <<
endl;
1783 const label nAllowableErrors,
1791 const word minThicknessName = moveDict.
get<
word>(
"minThicknessName");
1798 if (minThicknessName ==
"none")
1804 (minThicknessName ==
"none")
1812 pointDisplacement_.internalField(),
1821 forAll(extrudeStatus, pointI)
1823 if (
mag(patchDisp[pointI]) <= minThickness[pointI]+SMALL)
1831 calculateDisplacement(moveDict, minThickness, extrudeStatus, patchDisp);
1848 meshMover_.movePoints();
1851 meshMover_.correct();