44 void Foam::snappyLayerDriver::sumWeights
46 const PackedBoolList& isMasterEdge,
59 if (isMasterEdge.get(meshEdges[edgeI]) == 1)
61 const edge& e = edges[edgeI];
71 - pts[meshPoints[e[0]]]
74 scalar eWeight = 1.0/eMag;
76 invSumWeight[e[0]] += eWeight;
77 invSumWeight[e[1]] += eWeight;
90 forAll(invSumWeight, pointi)
92 scalar w = invSumWeight[pointi];
96 invSumWeight[pointi] = 1.0/w;
103 void Foam::snappyLayerDriver::smoothField
105 const motionSmoother& meshMover,
106 const PackedBoolList& isMasterPoint,
107 const PackedBoolList& isMasterEdge,
110 const label nSmoothDisp,
116 const labelList& meshPoints = pp.meshPoints();
129 Info<<
"shrinkMeshDistance : Smoothing field ..." <<
endl;
131 for (
label iter = 0; iter < nSmoothDisp; iter++)
155 average[pointi] < field[pointi]
156 &&
average[pointi] >= fieldMin[pointi]
159 field[pointi] =
average[pointi];
164 if ((iter % 10) == 0)
173 Info<<
" Iteration " << iter <<
" residual " << resid <<
endl;
270 void Foam::snappyLayerDriver::smoothPatchNormals
272 const motionSmoother& meshMover,
273 const PackedBoolList& isMasterPoint,
274 const PackedBoolList& isMasterEdge,
276 const label nSmoothDisp,
282 const labelList& meshPoints = pp.meshPoints();
297 Info<<
"shrinkMeshDistance : Smoothing normals ..." <<
endl;
299 for (
label iter = 0; iter < nSmoothDisp; iter++)
315 if ((iter % 10) == 0)
324 Info<<
" Iteration " << iter <<
" residual " << resid <<
endl;
332 normals[pointi] =
average[pointi];
333 normals[pointi] /=
mag(normals[pointi]) + vSmall;
340 void Foam::snappyLayerDriver::smoothNormals
342 const label nSmoothDisp,
343 const PackedBoolList& isMasterPoint,
344 const PackedBoolList& isMasterEdge,
350 Info<<
"shrinkMeshDistance : Smoothing normals in interior ..." <<
endl;
352 const fvMesh& mesh = meshRefiner_.
mesh();
353 const edgeList& edges = mesh.edges();
356 PackedBoolList isFixedPoint(mesh.nPoints());
361 label meshPointi = fixedPoints[i];
362 isFixedPoint.set(meshPointi, 1);
386 for (
label iter = 0; iter < nSmoothDisp; iter++)
402 if ((iter % 10) == 0)
410 Info<<
" Iteration " << iter <<
" residual " << resid <<
endl;
417 if (isFixedPoint.get(pointi) == 0)
421 normals[pointi] =
average[pointi];
422 normals[pointi] /=
mag(normals[pointi]) + vSmall;
431 bool Foam::snappyLayerDriver::isMaxEdge
433 const List<pointData>& pointWallDist,
438 const fvMesh& mesh = meshRefiner_.
mesh();
443 const edge& e = mesh.edges()[edgeI];
445 vector v0(points[e[0]] - pointWallDist[e[0]].origin());
446 scalar magV0(
mag(v0));
453 vector v1(points[e[1]] - pointWallDist[e[1]].origin());
454 scalar magV1(
mag(v1));
479 if ((pointWallDist[e[0]].v() & pointWallDist[e[1]].v()) < minCos)
492 void Foam::snappyLayerDriver::handleFeatureAngleLayerTerminations
495 const PackedBoolList& isMasterPoint,
498 List<extrudeMode>& extrudeStatus,
504 const fvMesh& mesh = meshRefiner_.
mesh();
509 boolList extrudedFaces(pp.size(),
true);
511 forAll(pp.localFaces(), facei)
513 const face& f = pp.localFaces()[facei];
519 extrudedFaces[facei] =
false;
534 List<List<point>> edgeFaceNormals(pp.nEdges());
535 List<List<bool>> edgeFaceExtrude(pp.nEdges());
539 const labelList& meshPoints = pp.meshPoints();
543 const labelList& eFaces = edgeFaces[edgeI];
545 edgeFaceNormals[edgeI].
setSize(eFaces.size());
546 edgeFaceExtrude[edgeI].setSize(eFaces.size());
549 label facei = eFaces[i];
550 edgeFaceNormals[edgeI][i] = faceNormals[facei];
551 edgeFaceExtrude[edgeI][i] = extrudedFaces[facei];
560 globalMeshData::ListPlusEqOp<List<point>>(),
569 globalMeshData::ListPlusEqOp<List<bool>>(),
574 forAll(edgeFaceNormals, edgeI)
576 const List<point>& eFaceNormals = edgeFaceNormals[edgeI];
577 const List<bool>& eFaceExtrude = edgeFaceExtrude[edgeI];
579 if (eFaceNormals.size() == 2)
581 const edge& e = pp.edges()[edgeI];
591 if (!eFaceExtrude[0] || !eFaceExtrude[1])
593 const vector& n0 = eFaceNormals[0];
594 const vector& n1 = eFaceNormals[1];
596 if ((n0 & n1) < minCos)
609 if (isMasterPoint[meshPoints[v0]])
625 if (isMasterPoint[meshPoints[v1]])
644 void Foam::snappyLayerDriver::findIsolatedRegions
646 const scalar minCosLayerTermination,
647 const PackedBoolList& isMasterPoint,
648 const PackedBoolList& isMasterEdge,
652 List<extrudeMode>& extrudeStatus,
657 const fvMesh& mesh = meshRefiner_.
mesh();
659 Info<<
"shrinkMeshDistance : Removing isolated regions ..." <<
endl;
662 label nPointCounter = 0;
668 handleFeatureAngleLayerTerminations
670 minCosLayerTermination,
681 syncPatchDisplacement
694 boolList extrudedFaces(pp.size(),
true);
695 forAll(pp.localFaces(), facei)
697 const face& f = pp.localFaces()[facei];
702 extrudedFaces[facei] =
false;
710 boolList keptPoints(pp.nPoints(),
false);
711 forAll(keptPoints, patchPointi)
713 const labelList& pFaces = pointFaces[patchPointi];
717 label facei = pFaces[i];
718 if (extrudedFaces[facei])
720 keptPoints[patchPointi] =
true;
737 forAll(keptPoints, patchPointi)
739 if (!keptPoints[patchPointi])
775 if (isMasterEdge.get(meshEdges[edgeI]) == 1)
777 const edge& e = edges[edgeI];
784 isolatedPoint[v0] += 1;
788 isolatedPoint[v1] += 1;
806 const face& f = pp.localFaces()[facei];
810 if (isolatedPoint[f[fp]] > 2)
816 bool allPointsExtruded =
true;
823 allPointsExtruded =
false;
828 if (allPointsExtruded)
850 reduce(nPointCounter, sumOp<label>());
851 Info<<
"Number isolated points extrusion stopped : "<< nPointCounter
862 void Foam::snappyLayerDriver::medialAxisSmoothingInfo
864 const motionSmoother& meshMover,
865 const label nSmoothNormals,
866 const label nSmoothSurfaceNormals,
867 const scalar minMedialAxisAngleCos,
868 const scalar featureAngle,
877 Info<<
"medialAxisSmoothingInfo :" 878 <<
" Calculate distance to Medial Axis ..." <<
endl;
880 const polyMesh& mesh = meshMover.mesh();
884 const labelList& meshPoints = pp.meshPoints();
912 UIndirectList<point>(meshPointNormals, meshPoints) = pointNormals;
928 nSmoothSurfaceNormals,
936 UIndirectList<point>(meshPointNormals, meshPoints) = pointNormals;
939 "smoothed pointNormals",
949 List<pointData> pointWallDist(mesh.nPoints());
952 int dummyTrackData = 0;
958 List<pointData> wallInfo(meshPoints.size());
960 forAll(meshPoints, patchPointi)
962 label pointi = meshPoints[patchPointi];
963 wallInfo[patchPointi] = pointData
968 pointNormals[patchPointi]
973 List<pointData> edgeWallDist(mesh.nEdges());
974 PointEdgeWave<pointData> wallDistCalc
981 mesh.globalData().nTotalPoints(),
988 wallDistCalc.getUnsetPoints(),
995 <<
"Walking did not visit all points." <<
nl 996 <<
" Did not visit " << nUnvisit
997 <<
" out of " << mesh.globalData().nTotalPoints()
998 <<
" points. This is not necessarily a problem" <<
nl 999 <<
" and might be due to faceZones splitting of part" 1000 <<
" of the domain." <<
nl <<
endl;
1012 forAll(pointWallDist, pointi)
1014 origin[pointi] = pointWallDist[pointi].origin();
1015 distSqr[pointi] = pointWallDist[pointi].distSqr();
1017 passiveV[pointi] = pointWallDist[pointi].v();
1029 List<pointData> pointMedialDist(mesh.nPoints());
1030 List<pointData> edgeMedialDist(mesh.nEdges());
1033 DynamicList<pointData> maxInfo(meshPoints.size());
1034 DynamicList<label> maxPoints(meshPoints.size());
1038 const edgeList& edges = mesh.edges();
1042 const edge& e = edges[edgeI];
1046 !pointWallDist[e[0]].valid(dummyTrackData)
1047 || !pointWallDist[e[1]].valid(dummyTrackData)
1052 else if (isMaxEdge(pointWallDist, edgeI, minMedialAxisAngleCos))
1059 vector eVec = e.vec(points);
1060 scalar eMag =
mag(eVec);
1066 const point& p0 = points[e[0]];
1067 const point& p1 = points[e[1]];
1068 scalar dist0 = (p0-pointWallDist[e[0]].origin()) & eVec;
1069 scalar dist1 = (pointWallDist[e[1]].origin()-p1) & eVec;
1070 scalar s = 0.5*(dist1+eMag+dist0);
1077 else if (s >= dist0+eMag)
1083 medialAxisPt = p0+(s-dist0)*eVec;
1088 label pointi = e[ep];
1090 if (!pointMedialDist[pointi].valid(dummyTrackData))
1092 maxPoints.append(pointi);
1098 magSqr(points[pointi]-medialAxisPt),
1103 pointMedialDist[pointi] = maxInfo.last();
1112 const polyBoundaryMesh& patches = mesh.boundaryMesh();
1119 const polyPatch& pp = patches[
patchi];
1121 meshMover.displacement().boundaryField()[
patchi];
1126 && !isA<emptyPolyPatch>(pp)
1127 && !adaptPatches.found(
patchi)
1130 const labelList& meshPoints = pp.meshPoints();
1137 if (pvf.fixesValue())
1140 Info<<
"Inserting all points on patch " << pp.name()
1145 label pointi = meshPoints[i];
1146 if (!pointMedialDist[pointi].valid(dummyTrackData))
1148 maxPoints.append(pointi);
1159 pointMedialDist[pointi] = maxInfo.last();
1170 Info<<
"Inserting points on patch " << pp.name()
1171 <<
" if angle to nearest layer patch > " 1172 << featureAngle <<
" degrees." <<
endl;
1179 label pointi = meshPoints[i];
1183 pointWallDist[pointi].valid(dummyTrackData)
1184 && !pointMedialDist[pointi].valid(dummyTrackData)
1190 -pointWallDist[pointi].v()
1193 if (cosAngle > featureAngleCos)
1198 maxPoints.append(pointi);
1209 pointMedialDist[pointi] = maxInfo.last();
1226 PointEdgeWave<pointData> medialDistCalc
1234 mesh.globalData().nTotalPoints(),
1239 forAll(pointMedialDist, pointi)
1241 medialDist[pointi] =
Foam::sqrt(pointMedialDist[pointi].distSqr());
1242 medialVec[pointi] = pointMedialDist[pointi].origin();
1256 if (!pointWallDist[i].valid(dummyTrackData))
1258 dispVec[i] =
vector(1, 0, 0);
1262 dispVec[i] = pointWallDist[i].v();
1282 forAll(medialRatio, pointi)
1284 if (!pointWallDist[pointi].valid(dummyTrackData))
1286 medialRatio[pointi] = 0.0;
1290 scalar wDist2 = pointWallDist[pointi].distSqr();
1291 scalar mDist = medialDist[pointi];
1293 if (wDist2 <
sqr(small) && mDist < small)
1300 medialRatio[pointi] = 0.0;
1304 medialRatio[pointi] = mDist / (
Foam::sqrt(wDist2) + mDist);
1319 Info<<
"medialAxisSmoothingInfo :" 1320 <<
" Writing:" <<
nl 1321 <<
" " << dispVec.name()
1322 <<
" : normalised direction of nearest displacement" <<
nl 1323 <<
" " << medialDist.name()
1324 <<
" : distance to medial axis" <<
nl 1325 <<
" " << medialVec.name()
1326 <<
" : nearest point on medial axis" <<
nl 1327 <<
" " << medialRatio.name()
1328 <<
" : ratio of medial distance to wall distance" <<
nl 1339 mesh.time().path()/meshRefiner_.
timeName()
1344 medialRatio.write();
1349 void Foam::snappyLayerDriver::shrinkMeshMedialDistance
1351 motionSmoother& meshMover,
1352 const dictionary& meshQualityDict,
1353 const List<labelPair>& baffles,
1354 const label nSmoothPatchThickness,
1355 const label nSmoothDisplacement,
1356 const scalar maxThicknessToMedialRatio,
1357 const label nAllowableErrors,
1359 const scalar minCosLayerTermination,
1369 List<extrudeMode>& extrudeStatus,
1374 Info<<
"shrinkMeshMedialDistance : Smoothing using Medial Axis ..." <<
endl;
1376 const polyMesh& mesh = meshMover.mesh();
1379 const labelList& meshPoints = pp.meshPoints();
1397 thickness =
mag(patchDisp);
1399 forAll(thickness, patchPointi)
1401 if (extrudeStatus[patchPointi] ==
NOEXTRUDE)
1403 thickness[patchPointi] = 0.0;
1407 label numThicknessRatioExclude = 0;
1412 autoPtr<OBJstream> str;
1420 /
"thicknessRatioExcludePoints_" 1425 Info<<
"Writing points with too large an extrusion distance to " 1426 << str().name() <<
endl;
1429 autoPtr<OBJstream> medialVecStr;
1437 /
"thicknessRatioExcludeMedialVec_" 1442 Info<<
"Writing points with too large an extrusion distance to " 1443 << medialVecStr().name() <<
endl;
1446 forAll(meshPoints, patchPointi)
1448 if (extrudeStatus[patchPointi] !=
NOEXTRUDE)
1450 label pointi = meshPoints[patchPointi];
1454 scalar mDist = medialDist[pointi];
1455 scalar thicknessRatio = thickness[patchPointi]/(mDist+vSmall);
1460 patchDisp[patchPointi]
1461 / (
mag(patchDisp[patchPointi]) + vSmall);
1462 vector mVec = mesh.points()[pointi]-medialVec[pointi];
1463 mVec /=
mag(mVec)+vSmall;
1464 thicknessRatio *= (n&mVec);
1466 if (thicknessRatio > maxThicknessToMedialRatio)
1471 Pout<<
"truncating displacement at " 1472 << mesh.points()[pointi]
1473 <<
" from " << thickness[patchPointi]
1477 minThickness[patchPointi]
1478 +thickness[patchPointi]
1480 <<
" medial direction:" << mVec
1481 <<
" extrusion direction:" << n
1482 <<
" with thicknessRatio:" << thicknessRatio
1486 thickness[patchPointi] =
1487 0.5*(minThickness[patchPointi]+thickness[patchPointi]);
1489 patchDisp[patchPointi] = thickness[patchPointi]*
n;
1491 if (isMasterPoint[pointi])
1493 numThicknessRatioExclude++;
1498 const point& pt = mesh.points()[pointi];
1499 str().write(
linePointRef(pt, pt+patchDisp[patchPointi]));
1501 if (medialVecStr.valid())
1503 const point& pt = mesh.points()[pointi];
1504 medialVecStr().write
1517 reduce(numThicknessRatioExclude, sumOp<label>());
1518 Info<<
"shrinkMeshMedialDistance : Reduce layer thickness at " 1519 << numThicknessRatioExclude
1520 <<
" nodes where thickness to medial axis distance is large " <<
endl;
1527 minCosLayerTermination,
1540 forAll(thickness, patchPointi)
1542 if (extrudeStatus[patchPointi] ==
NOEXTRUDE)
1544 thickness[patchPointi] = 0.0;
1557 nSmoothPatchThickness,
1563 int dummyTrackData = 0;
1565 List<pointData> pointWallDist(mesh.nPoints());
1572 List<pointData> edgeWallDist(mesh.nEdges());
1573 labelList wallPoints(meshPoints.size());
1576 List<pointData> wallInfo(meshPoints.size());
1578 forAll(meshPoints, patchPointi)
1580 label pointi = meshPoints[patchPointi];
1581 wallPoints[patchPointi] = pointi;
1582 wallInfo[patchPointi] = pointData
1586 thickness[patchPointi],
1592 PointEdgeWave<pointData> wallDistCalc
1599 mesh.globalData().nTotalPoints(),
1607 forAll(displacement, pointi)
1609 if (!pointWallDist[pointi].valid(dummyTrackData))
1611 displacement[pointi] =
Zero;
1620 displacement[pointi] =
1621 -medialRatio[pointi]
1622 *pointWallDist[pointi].s()
1630 if (nSmoothDisplacement > 0)
1644 Info<<
"shrinkMeshDistance : Smoothing displacement ..." <<
endl;
1646 const scalar lambda = 0.33;
1647 const scalar mu = -0.34;
1650 for (
label iter = 0; iter < nSmoothDisplacement; iter++)
1667 if (medialRatio[i] > small && medialRatio[i] < 1-small)
1670 (1-
lambda)*displacement[i]
1691 if (medialRatio[i] > small && medialRatio[i] < 1-small)
1693 displacement[i] = (1-
mu)*displacement[i]+mu*
average[i];
1699 if ((iter % 10) == 0)
1707 Info<<
" Iteration " << iter <<
" residual " << resid
1715 meshMover.setDisplacementPatchFields();
1726 forAll(pointWallDist, pointi)
1728 pWallDist[pointi] = pointWallDist[pointi].s();
1738 "displacement BEFORE",
1743 meshMover.correctBoundaryConditions(displacement);
1746 "displacement AFTER",
1755 const_cast<Time&
>(mesh.time())++;
1756 Info<<
"Writing wanted-displacement mesh (possibly illegal) to " 1762 meshMover.movePoints();
1780 mesh.time().path()/meshRefiner_.
timeName()
1784 medialRatio.write();
1789 meshMover.movePoints();
1796 Info<<
"shrinkMeshMedialDistance : Moving mesh ..." <<
endl;
1797 scalar oldErrorReduction = -1;
1799 for (
label iter = 0; iter < 2*nSnap ; iter++)
1801 Info<<
"Iteration " << iter <<
endl;
1804 Info<<
"Displacement scaling for error reduction set to 0." <<
endl;
1805 oldErrorReduction = meshMover.setErrorReduction(0.0);
1814 meshMover.paramDict(),
1821 Info<<
"shrinkMeshMedialDistance : Successfully moved mesh" <<
endl;
1826 if (oldErrorReduction >= 0)
1828 meshMover.setErrorReduction(oldErrorReduction);
1831 Info<<
"shrinkMeshMedialDistance : Finished moving mesh ..." <<
endl;
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
List< labelList > labelListList
A List of labelList.
#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.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Unit conversion functions.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
static writeType writeLevel()
Get/set write level.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Vector< scalar > vector
A scalar version of the templated Vector.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
const fvMesh & mesh() const
Reference to mesh.
virtual const pointField & points() const
Return raw points.
List< bool > boolList
Bool container classes.
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
vectorField pointField
pointField is a vectorField.
dimensionedScalar cos(const dimensionedScalar &ds)
line< point, const point & > linePointRef
Line using referred points.
Do not extrude. No layers added.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
List< label > labelList
A List of labels.
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
word timeName() const
Replacement for Time::timeName() : return oldInstance (if.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
volScalarField scalarField(fieldObject, mesh)
pointPatchField< vector > pointPatchVectorField
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
void setInstance(const fileName &)
Set the instance for mesh files.
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
void setSize(const label)
Reset size of List.
static T gAverage(const PackedBoolList &isMasterElem, const UList< T > &values)
Helper: calculate average.
vector point
Point is a vector.
const dimensionedScalar & mu
Atomic mass unit.
#define WarningInFunction
Report a warning using Foam::Warning.
static void testSyncPointList(const string &msg, const polyMesh &mesh, const List< scalar > &fld)
prefixOSstream Pout(cout, "Pout")
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
GeometricField< vector, pointPatchField, pointMesh > pointVectorField
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)
bool write() const
Write mesh and all data.