40 void Foam::autoLayerDriver::sumWeights
53 if (isMasterEdge.get(meshEdges[edgeI]) == 1)
55 const edge&
e = edges[edgeI];
59 invSumWeight[e[0]] += eWeight;
60 invSumWeight[e[1]] += eWeight;
74 forAll(invSumWeight, pointI)
76 scalar w = invSumWeight[pointI];
80 invSumWeight[pointI] = 1.0/w;
87 void Foam::autoLayerDriver::smoothField
89 const motionSmoother& meshMover,
93 const label nSmoothDisp,
99 const labelList& meshPoints = pp.meshPoints();
112 Info<<
"shrinkMeshDistance : Smoothing field ..." <<
endl;
114 for (label iter = 0; iter < nSmoothDisp; iter++)
138 average[pointI] < field[pointI]
139 &&
average[pointI] >= fieldMin[pointI]
142 field[pointI] =
average[pointI];
147 if ((iter % 10) == 0)
149 Info<<
" Iteration " << iter <<
" residual "
159 void Foam::autoLayerDriver::smoothPatchNormals
161 const motionSmoother& meshMover,
164 const label nSmoothDisp,
170 const labelList& meshPoints = pp.meshPoints();
185 Info<<
"shrinkMeshDistance : Smoothing normals ..." <<
endl;
187 for (label iter = 0; iter < nSmoothDisp; iter++)
203 if ((iter % 10) == 0)
205 Info<<
" Iteration " << iter <<
" residual "
216 normals[pointI] =
average[pointI];
217 normals[pointI] /=
mag(normals[pointI]) + VSMALL;
224 void Foam::autoLayerDriver::smoothNormals
226 const label nSmoothDisp,
233 Info<<
"shrinkMeshDistance : Smoothing normals ..." <<
endl;
235 const fvMesh&
mesh = meshRefiner_.mesh();
236 const edgeList& edges = mesh.edges();
244 label meshPointI = fixedPoints[i];
245 isFixedPoint.set(meshPointI, 1);
264 Info<<
"shrinkMeshDistance : Smoothing normals in interior ..." <<
endl;
266 for (label iter = 0; iter < nSmoothDisp; iter++)
282 if ((iter % 10) == 0)
284 Info<<
" Iteration " << iter <<
" residual "
294 if (isFixedPoint.get(pointI) == 0)
298 normals[pointI] =
average[pointI];
299 normals[pointI] /=
mag(normals[pointI]) + VSMALL;
308 bool Foam::autoLayerDriver::isMaxEdge
310 const List<pointData>& pointWallDist,
315 const fvMesh& mesh = meshRefiner_.mesh();
320 const edge&
e = mesh.edges()[edgeI];
322 vector v0(points[e[0]] - pointWallDist[e[0]].origin());
323 scalar magV0(
mag(v0));
330 vector v1(points[e[1]] - pointWallDist[e[1]].origin());
331 scalar magV1(
mag(
v1));
342 if ((v0 &
v1) < minCos)
355 void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
359 List<extrudeMode>& extrudeStatus,
368 boolList extrudedFaces(pp.size(),
true);
370 forAll(pp.localFaces(), faceI)
372 const face&
f = pp.localFaces()[faceI];
376 if (extrudeStatus[f[fp]] == NOEXTRUDE)
378 extrudedFaces[faceI] =
false;
389 forAll(pp.edgeFaces(), edgeI)
391 const labelList& eFaces = pp.edgeFaces()[edgeI];
393 if (eFaces.size() == 2)
395 const edge& e = pp.edges()[edgeI];
401 extrudeStatus[v0] != NOEXTRUDE
402 || extrudeStatus[v1] != NOEXTRUDE
405 if (!extrudedFaces[eFaces[0]] || !extrudedFaces[eFaces[1]])
407 const vector& n0 = pp.faceNormals()[eFaces[0]];
408 const vector& n1 = pp.faceNormals()[eFaces[1]];
410 if ((n0 & n1) < minCos)
448 void Foam::autoLayerDriver::findIsolatedRegions
453 const scalar minCosLayerTermination,
455 List<extrudeMode>& extrudeStatus,
460 const fvMesh& mesh = meshRefiner_.mesh();
462 Info<<
"shrinkMeshDistance : Removing isolated regions ..." <<
endl;
465 label nPointCounter = 0;
471 handleFeatureAngleLayerTerminations
474 minCosLayerTermination,
486 boolList extrudedFaces(pp.size(),
true);
487 forAll(pp.localFaces(), faceI)
489 const face& f = pp.localFaces()[faceI];
492 if (extrudeStatus[f[fp]] == NOEXTRUDE)
494 extrudedFaces[faceI] =
false;
502 boolList keptPoints(pp.nPoints(),
false);
503 forAll(keptPoints, patchPointI)
509 label faceI = pFaces[i];
510 if (extrudedFaces[faceI])
512 keptPoints[patchPointI] =
true;
530 forAll(keptPoints, patchPointI)
532 if (!keptPoints[patchPointI])
568 if (isMasterEdge.get(meshEdges[edgeI]) == 1)
570 const edge& e = edges[edgeI];
575 if (extrudeStatus[v1] != NOEXTRUDE)
577 isolatedPoint[v0] += 1;
579 if (extrudeStatus[v0] != NOEXTRUDE)
581 isolatedPoint[
v1] += 1;
600 const face& f = pp.localFaces()[faceI];
604 if (isolatedPoint[f[fp]] > 2)
610 bool allPointsExtruded =
true;
615 if (extrudeStatus[f[fp]] == NOEXTRUDE)
617 allPointsExtruded =
false;
622 if (allPointsExtruded)
645 reduce(nPointCounter, sumOp<label>());
646 Info<<
"Number isolated points extrusion stopped : "<< nPointCounter<<
endl;
656 void Foam::autoLayerDriver::medialAxisSmoothingInfo
658 const motionSmoother& meshMover,
659 const label nSmoothNormals,
660 const label nSmoothSurfaceNormals,
661 const scalar minMedianAxisAngleCos,
669 Info<<
"medialAxisSmoothingInfo :"
670 <<
" Calculate distance to Medial Axis ..." <<
endl;
672 const polyMesh& mesh = meshMover.mesh();
677 const labelList& meshPoints = pp.meshPoints();
687 forAll(meshEdges, patchEdgeI)
689 const edge& e = pp.edges()[patchEdgeI];
691 label v0 = pp.meshPoints()[e[0]];
692 label v1 = pp.meshPoints()[e[1]];
696 mesh.pointEdges()[v0],
710 forAll(faceNormals, faceI)
712 const face& f = pp.localFaces()[faceI];
716 pointNormals[f[fp]] += faceNormals[faceI];
717 nPointFaces[f[fp]] ++;
743 pointNormals[i] /= nPointFaces[i];
753 nSmoothSurfaceNormals,
762 List<pointData> pointWallDist(mesh.nPoints());
768 List<pointData> wallInfo(meshPoints.size());
770 forAll(meshPoints, patchPointI)
772 label pointI = meshPoints[patchPointI];
773 wallInfo[patchPointI] = pointData
778 pointNormals[patchPointI]
783 List<pointData> edgeWallDist(mesh.nEdges());
784 PointEdgeWave<pointData> wallDistCalc
791 mesh.globalData().nTotalPoints()
798 List<pointData> pointMedialDist(mesh.nPoints());
799 List<pointData> edgeMedialDist(mesh.nEdges());
802 DynamicList<pointData> maxInfo(meshPoints.size());
803 DynamicList<label> maxPoints(meshPoints.size());
807 const edgeList& edges = mesh.edges();
811 if (isMaxEdge(pointWallDist, edgeI, minMedianAxisAngleCos))
815 const edge& e = edges[edgeI];
819 label pointI = e[ep];
821 if (!pointMedialDist[pointI].valid())
823 maxPoints.append(pointI);
834 pointMedialDist[pointI] = maxInfo[maxInfo.size()-1];
842 const polyBoundaryMesh&
patches = mesh.boundaryMesh();
848 const polyPatch& pp = patches[patchI];
853 && !isA<emptyPolyPatch>(pp)
854 && !adaptPatches.found(patchI)
857 Info<<
"Inserting points on patch " << pp.name() <<
endl;
859 const labelList& meshPoints = pp.meshPoints();
863 label pointI = meshPoints[i];
865 if (!pointMedialDist[pointI].valid())
867 maxPoints.append(pointI);
878 pointMedialDist[pointI] = maxInfo[maxInfo.size()-1];
888 PointEdgeWave<pointData> medialDistCalc
896 mesh.globalData().nTotalPoints()
900 forAll(pointMedialDist, pointI)
902 medialDist[pointI] =
Foam::sqrt(pointMedialDist[pointI].distSqr());
909 dispVec[i] = pointWallDist[i].v();
913 smoothNormals(nSmoothNormals, isMasterEdge, meshPoints, dispVec);
916 forAll(medialRatio, pointI)
918 scalar wDist2 = pointWallDist[pointI].distSqr();
919 scalar mDist = medialDist[pointI];
921 if (wDist2 <
sqr(SMALL) && mDist < SMALL)
923 medialRatio[pointI] = 0.0;
927 medialRatio[pointI] = mDist / (
Foam::sqrt(wDist2) + mDist);
933 Info<<
"medialAxisSmoothingInfo :"
935 <<
" " << dispVec.name()
936 <<
" : normalised direction of nearest displacement" <<
nl
937 <<
" " << medialDist.name()
938 <<
" : distance to medial axis" <<
nl
939 <<
" " << medialRatio.name()
940 <<
" : ratio of medial distance to wall distance" <<
nl
949 void Foam::autoLayerDriver::shrinkMeshMedialDistance
951 motionSmoother& meshMover,
952 const dictionary& meshQualityDict,
953 const label nSmoothThickness,
954 const scalar maxThicknessToMedialRatio,
955 const label nAllowableErrors,
957 const scalar minCosLayerTermination,
966 List<extrudeMode>& extrudeStatus,
971 Info<<
"shrinkMeshMedialDistance : Smoothing using Medial Axis ..." <<
endl;
973 const polyMesh& mesh = meshMover.mesh();
976 const labelList& meshPoints = pp.meshPoints();
983 forAll(meshEdges, patchEdgeI)
985 const edge& e = pp.edges()[patchEdgeI];
987 label v0 = pp.meshPoints()[e[0]];
988 label v1 = pp.meshPoints()[e[1]];
992 mesh.pointEdges()[v0],
1001 thickness =
mag(patchDisp);
1003 forAll(thickness, patchPointI)
1005 if (extrudeStatus[patchPointI] == NOEXTRUDE)
1007 thickness[patchPointI] = 0.0;
1011 label numThicknessRatioExclude = 0;
1014 forAll(meshPoints, patchPointI)
1016 if (extrudeStatus[patchPointI] != NOEXTRUDE)
1018 label pointI = meshPoints[patchPointI];
1020 scalar mDist = medialDist[pointI];
1022 scalar thicknessRatio = thickness[patchPointI]/(mDist+VSMALL);
1024 if (thicknessRatio > maxThicknessToMedialRatio)
1027 thickness[patchPointI] =
1028 0.5*(minThickness[patchPointI]+thickness[patchPointI]);
1030 patchDisp[patchPointI] =
1031 thickness[patchPointI]
1032 * patchDisp[patchPointI]
1033 / (
mag(patchDisp[patchPointI]) + VSMALL);
1034 numThicknessRatioExclude++;
1039 reduce(numThicknessRatioExclude, sumOp<label>());
1040 Info<<
"shrinkMeshMedialDistance : Reduce layer thickness at "
1041 << numThicknessRatioExclude
1042 <<
" nodes where thickness to medial axis distance is large " <<
endl;
1051 minCosLayerTermination,
1069 List<pointData> pointWallDist(mesh.nPoints());
1076 List<pointData> edgeWallDist(mesh.nEdges());
1077 labelList wallPoints(meshPoints.size());
1080 List<pointData> wallInfo(meshPoints.size());
1082 forAll(meshPoints, patchPointI)
1084 label pointI = meshPoints[patchPointI];
1085 wallPoints[patchPointI] = pointI;
1086 wallInfo[patchPointI] = pointData
1090 thickness[patchPointI],
1096 PointEdgeWave<pointData> wallDistCalc
1103 mesh.globalData().nTotalPoints()
1110 forAll(displacement, pointI)
1117 displacement[pointI] =
1118 -medialRatio[pointI]
1119 *pointWallDist[pointI].s()
1126 Info<<
"shrinkMeshMedialDistance : Moving mesh ..." <<
endl;
1127 scalar oldErrorReduction = -1;
1129 for (label iter = 0; iter < 2*nSnap ; iter++)
1131 Info<<
"Iteration " << iter <<
endl;
1134 Info<<
"Displacement scaling for error reduction set to 0." <<
endl;
1135 oldErrorReduction = meshMover.setErrorReduction(0.0);
1144 meshMover.paramDict(),
1151 Info<<
"shrinkMeshMedialDistance : Successfully moved mesh" <<
endl;
1156 if (oldErrorReduction >= 0)
1158 meshMover.setErrorReduction(oldErrorReduction);
1161 Info<<
"shrinkMeshMedialDistance : Finished moving mesh ..." <<
endl;