61 Foam::label Foam::autoLayerDriver::mergePatchFacesUndo
64 const scalar concaveCos,
65 const dictionary& motionDict
71 combineFaces faceCombiner(mesh,
true);
74 labelHashSet boundaryCells(mesh.nFaces()-mesh.nInternalFaces());
79 const polyBoundaryMesh&
patches = mesh.boundaryMesh();
83 label patchI = patchIDs[i];
85 const polyPatch& patch = patches[patchI];
91 boundaryCells.insert(mesh.faceOwner()[patch.start()+i]);
100 faceCombiner.getMergeSets
108 label nFaceSets =
returnReduce(allFaceSets.size(), sumOp<label>());
110 Info<<
"Merging " << nFaceSets <<
" sets of faces." <<
nl <<
endl;
116 faceSet allSets(mesh,
"allFaceSets", allFaceSets.size());
119 forAll(allFaceSets[setI], i)
121 allSets.insert(allFaceSets[setI][i]);
124 Pout<<
"Writing all faces to be merged to set "
125 << allSets.objectPath() <<
endl;
131 polyTopoChange meshMod(mesh);
134 faceCombiner.setRefinement(allFaceSets, meshMod);
139 faceCombiner.savedPointLabels(),
145 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh,
false,
true);
148 mesh.updateMesh(map);
151 if (map().hasMotionPoints())
153 mesh.movePoints(map().preMotionPoints());
166 faceCombiner.updateMesh(map);
171 for (label iteration = 0; iteration < 100; iteration++)
174 <<
"Undo iteration " << iteration <<
nl
175 <<
"----------------" <<
endl;
185 mesh.nFaces()-mesh.nInternalFaces()
223 Pout<<
"Writing all faces in error to faceSet "
224 << errorFaces.objectPath() <<
nl <<
endl;
232 DynamicList<label> mastersToRestore(allFaceSets.size());
236 label masterFaceI = faceCombiner.masterFace()[setI];
238 if (masterFaceI != -1)
240 label masterCellII = mesh.faceOwner()[masterFaceI];
242 const cell& cFaces = mesh.cells()[masterCellII];
246 if (errorFaces.found(cFaces[i]))
248 mastersToRestore.append(masterFaceI);
254 mastersToRestore.shrink();
258 mastersToRestore.size(),
262 Info<<
"Masters that need to be restored:"
267 faceSet restoreSet(mesh,
"mastersToRestore", mastersToRestore);
268 Pout<<
"Writing all " << mastersToRestore.size()
269 <<
" masterfaces to be restored to set "
270 << restoreSet.objectPath() <<
endl;
285 polyTopoChange meshMod(mesh);
289 Map<label> restoredPoints(0);
290 Map<label> restoredFaces(0);
291 Map<label> restoredCells(0);
293 faceCombiner.setUnrefinement
303 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh,
false,
true);
306 mesh.updateMesh(map);
309 if (map().hasMotionPoints())
311 mesh.movePoints(map().preMotionPoints());
324 faceCombiner.updateMesh(map);
346 Pout<<
"Writing merged-faces mesh to time "
353 Info<<
"No faces merged ..." <<
endl;
363 removePoints& pointRemover,
367 fvMesh& mesh = meshRefiner_.mesh();
370 polyTopoChange meshMod(mesh);
372 pointRemover.setRefinement(pointCanBeDeleted, meshMod);
375 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh,
false,
true);
378 mesh.updateMesh(map);
381 if (map().hasMotionPoints())
383 mesh.movePoints(map().preMotionPoints());
391 if (meshRefiner_.overwrite())
393 mesh.setInstance(meshRefiner_.oldInstance());
396 pointRemover.updateMesh(map);
397 meshRefiner_.updateMesh(map,
labelList(0));
406 removePoints& pointRemover,
410 fvMesh& mesh = meshRefiner_.mesh();
413 polyTopoChange meshMod(mesh);
417 pointRemover.getUnrefimentSet
425 pointRemover.setUnrefinement
433 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh,
false,
true);
436 mesh.updateMesh(map);
439 if (map().hasMotionPoints())
441 mesh.movePoints(map().preMotionPoints());
449 if (meshRefiner_.overwrite())
451 mesh.setInstance(meshRefiner_.oldInstance());
454 pointRemover.updateMesh(map);
455 meshRefiner_.updateMesh(map,
labelList(0));
469 const fvMesh& mesh = meshRefiner_.mesh();
472 boolList selected(mesh.nFaces(),
false);
476 label faceI = candidateFaces[i];
478 if (
set.found(faceI))
480 selected[faceI] =
true;
493 return selectedFaces;
503 const fvMesh& mesh = meshRefiner_.mesh();
505 boolList selected(mesh.nFaces(),
false);
509 label faceI = iter.key();
511 label own = mesh.faceOwner()[faceI];
513 const cell& ownFaces = mesh.cells()[own];
514 forAll(ownFaces, ownFaceI)
516 selected[ownFaces[ownFaceI]] =
true;
519 if (mesh.isInternalFace(faceI))
521 label nbr = mesh.faceNeighbour()[faceI];
523 const cell& nbrFaces = mesh.cells()[nbr];
524 forAll(nbrFaces, nbrFaceI)
526 selected[nbrFaces[nbrFaceI]] =
true;
543 Foam::label Foam::autoLayerDriver::mergeEdgesUndo
546 const dictionary& motionDict
549 fvMesh& mesh = meshRefiner_.mesh();
552 <<
"Merging all points on surface that" <<
nl
553 <<
"- are used by only two boundary faces and" <<
nl
554 <<
"- make an angle with a cosine of more than " << minCos
558 removePoints pointRemover(mesh,
true);
562 label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
566 Info<<
"Removing " << nRemove
567 <<
" straight edge points ..." <<
nl <<
endl;
572 doRemovePoints(pointRemover, pointCanBeDeleted);
575 for (label iteration = 0; iteration < 100; iteration++)
578 <<
"Undo iteration " << iteration <<
nl
579 <<
"----------------" <<
endl;
589 mesh.nFaces()-mesh.nInternalFaces()
626 Pout<<
"**Writing all faces in error to faceSet "
627 << errorFaces.objectPath() <<
nl <<
endl;
635 pointRemover.savedFaceLabels(),
640 label n =
returnReduce(masterErrorFaces.size(), sumOp<label>());
642 Info<<
"Detected " << n
643 <<
" error faces on boundaries that have been merged."
644 <<
" These will be restored to their original faces." <<
nl
653 <<
" error faces in mesh."
654 <<
" Restoring neighbours of faces in error." <<
nl
665 doRestorePoints(pointRemover, expandedErrorFaces);
671 doRestorePoints(pointRemover, masterErrorFaces);
676 Pout<<
"Writing merged-edges mesh to time "
677 << meshRefiner_.timeName() <<
nl <<
endl;
683 Info<<
"No straight edges simplified and no points removed ..." <<
endl;
691 void Foam::autoLayerDriver::dumpDisplacement
693 const fileName& prefix,
696 const List<extrudeMode>& extrudeStatus
699 OFstream dispStr(prefix +
"_disp.obj");
700 Info<<
"Writing all displacements to " << dispStr.name() <<
nl <<
endl;
704 forAll(patchDisp, patchPointI)
706 const point& pt = pp.localPoints()[patchPointI];
711 dispStr <<
"l " << vertI-1 <<
' ' << vertI <<
nl;
715 OFstream illStr(prefix +
"_illegal.obj");
716 Info<<
"Writing invalid displacements to " << illStr.name() <<
nl <<
endl;
720 forAll(patchDisp, patchPointI)
722 if (extrudeStatus[patchPointI] != EXTRUDE)
724 const point& pt = pp.localPoints()[patchPointI];
729 illStr <<
"l " << vertI-1 <<
' ' << vertI <<
nl;
737 void Foam::autoLayerDriver::checkManifold
740 pointSet& nonManifoldPoints
744 fp.checkPointManifold(
false, &nonManifoldPoints);
751 const labelList& eFaces = edgeFaces[edgeI];
753 if (eFaces.size() > 2)
755 const edge&
e = fp.edges()[edgeI];
757 nonManifoldPoints.insert(fp.meshPoints()[e[0]]);
758 nonManifoldPoints.insert(fp.meshPoints()[e[1]]);
764 void Foam::autoLayerDriver::checkMeshManifold()
const
766 const fvMesh& mesh = meshRefiner_.mesh();
768 Info<<
nl <<
"Checking mesh manifoldness ..." <<
endl;
771 labelList outsideFaces(mesh.nFaces() - mesh.nInternalFaces());
773 for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
775 outsideFaces[faceI - mesh.nInternalFaces()] = faceI;
778 pointSet nonManifoldPoints
790 IndirectList<face>(mesh.faces(), outsideFaces),
796 label nNonManif =
returnReduce(nonManifoldPoints.size(), sumOp<label>());
800 Info<<
"Outside of mesh is multiply connected across edges or"
802 <<
"This is not a fatal error but might cause some unexpected"
803 <<
" behaviour." <<
nl
804 <<
"Writing " << nNonManif
805 <<
" points where this happens to pointSet "
806 << nonManifoldPoints.name() <<
endl;
808 nonManifoldPoints.
write();
816 bool Foam::autoLayerDriver::unmarkExtrusion
818 const label patchPointI,
821 List<extrudeMode>& extrudeStatus
824 if (extrudeStatus[patchPointI] == EXTRUDE)
826 extrudeStatus[patchPointI] = NOEXTRUDE;
827 patchNLayers[patchPointI] = 0;
831 else if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
833 extrudeStatus[patchPointI] = NOEXTRUDE;
834 patchNLayers[patchPointI] = 0;
846 bool Foam::autoLayerDriver::unmarkExtrusion
848 const face& localFace,
851 List<extrudeMode>& extrudeStatus
854 bool unextruded =
false;
877 void Foam::autoLayerDriver::handleNonManifolds
883 List<extrudeMode>& extrudeStatus
886 const fvMesh& mesh = meshRefiner_.mesh();
888 Info<<
nl <<
"Handling non-manifold points ..." <<
endl;
891 Info<<
nl <<
"Checking patch manifoldness ..." <<
endl;
893 pointSet nonManifoldPoints(mesh,
"nonManifoldPoints", pp.nPoints());
896 checkManifold(pp, nonManifoldPoints);
898 label nNonManif =
returnReduce(nonManifoldPoints.size(), sumOp<label>());
900 Info<<
"Outside of local patch is multiply connected across edges or"
901 <<
" points at " << nNonManif <<
" points." <<
endl;
904 const labelList& meshPoints = pp.meshPoints();
909 const labelHashSet sharedEdgeSet(mesh.globalData().sharedEdgeLabels());
913 if (sharedEdgeSet.found(meshEdges[edgeI]))
915 const edge& e = mesh.edges()[meshEdges[edgeI]];
917 Pout<<
"Disabling extrusion at edge "
918 << mesh.points()[e[0]] << mesh.points()[e[1]]
919 <<
" since it is non-manifold across coupled patches."
922 nonManifoldPoints.insert(e[0]);
923 nonManifoldPoints.insert(e[1]);
935 nNonManif =
returnReduce(nonManifoldPoints.size(), sumOp<label>());
939 Info<<
"Outside of patches is multiply connected across edges or"
940 <<
" points at " << nNonManif <<
" points." <<
nl
941 <<
"Writing " << nNonManif
942 <<
" points where this happens to pointSet "
943 << nonManifoldPoints.name() <<
nl
944 <<
"and setting layer thickness to zero on these points."
947 nonManifoldPoints.
write();
949 forAll(meshPoints, patchPointI)
951 if (nonManifoldPoints.found(meshPoints[patchPointI]))
964 Info<<
"Set displacement to zero for all " << nNonManif
965 <<
" non-manifold points" <<
endl;
970 void Foam::autoLayerDriver::handleFeatureAngle
977 List<extrudeMode>& extrudeStatus
980 const fvMesh& mesh = meshRefiner_.mesh();
982 Info<<
nl <<
"Handling feature edges ..." <<
endl;
984 if (minCos < 1-SMALL)
993 const labelList& eFaces = pp.edgeFaces()[edgeI];
995 label meshEdgeI = meshEdges[edgeI];
1001 edgeNormal[meshEdgeI],
1002 pp.faceNormals()[eFaces[i]]
1017 autoPtr<OFstream> str;
1020 str.reset(
new OFstream(mesh.time().path()/
"featureEdges.obj"));
1021 Info<<
"Writing feature edges to " << str().name() <<
endl;
1030 const labelList& eFaces = pp.edgeFaces()[edgeI];
1032 label meshEdgeI = meshEdges[edgeI];
1034 const vector& n = edgeNormal[meshEdgeI];
1038 scalar
cos = n & pp.faceNormals()[eFaces[0]];
1042 const edge& e = pp.edges()[edgeI];
1067 str()<<
"l " << vertI-1 <<
' ' << vertI <<
nl;
1073 Info<<
"Set displacement to zero for points on "
1075 <<
" feature edges" <<
endl;
1084 void Foam::autoLayerDriver::handleWarpedFaces
1087 const scalar faceRatio,
1088 const scalar edge0Len,
1092 List<extrudeMode>& extrudeStatus
1095 const fvMesh& mesh = meshRefiner_.mesh();
1097 Info<<
nl <<
"Handling cells with warped patch faces ..." <<
nl;
1101 label nWarpedFaces = 0;
1105 const face&
f = pp[i];
1109 label faceI = pp.addressing()[i];
1111 label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
1112 scalar edgeLen = edge0Len/(1<<ownLevel);
1115 const point& fc = mesh.faceCentres()[faceI];
1116 const vector& fn = pp.faceNormals()[i];
1122 vector n = points[f[fp]] - fc;
1123 vProj[fp] = (n & fn);
1127 scalar minVal =
min(vProj);
1128 scalar maxVal =
max(vProj);
1130 if ((maxVal - minVal) > faceRatio * edgeLen)
1149 Info<<
"Set displacement to zero on "
1151 <<
" warped faces since layer would be > " << faceRatio
1152 <<
" of the size of the bounding box." <<
endl;
1260 void Foam::autoLayerDriver::setNumLayers
1267 List<extrudeMode>& extrudeStatus
1270 const fvMesh& mesh = meshRefiner_.mesh();
1272 Info<< nl <<
"Handling points with inconsistent layer specification ..."
1277 labelList maxLayers(patchNLayers.size(), labelMin);
1278 labelList minLayers(patchNLayers.size(), labelMax);
1282 label patchI = patchIDs[i];
1284 const labelList& meshPoints = mesh.boundaryMesh()[patchI].meshPoints();
1286 label wantedLayers = patchToNLayers[patchI];
1288 forAll(meshPoints, patchPointI)
1290 label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
1292 maxLayers[ppPointI] =
max(wantedLayers, maxLayers[ppPointI]);
1293 minLayers[ppPointI] =
min(wantedLayers, minLayers[ppPointI]);
1323 if (maxLayers[i] == labelMin || minLayers[i] == labelMax)
1326 <<
"Patchpoint:" << i <<
" coord:" << pp.localPoints()[i]
1327 <<
" maxLayers:" << maxLayers
1328 <<
" minLayers:" << minLayers
1331 else if (maxLayers[i] == minLayers[i])
1334 patchNLayers[i] = maxLayers[i];
1352 patchNLayers[i] = maxLayers[i];
1365 void Foam::autoLayerDriver::growNoExtrusion
1370 List<extrudeMode>& extrudeStatus
1373 Info<< nl <<
"Growing non-extrusion points by one layer ..." <<
endl;
1375 List<extrudeMode> grownExtrudeStatus(extrudeStatus);
1377 const faceList& localFaces = pp.localFaces();
1381 forAll(localFaces, faceI)
1383 const face& f = localFaces[faceI];
1385 bool hasSqueeze =
false;
1388 if (extrudeStatus[f[fp]] == NOEXTRUDE)
1402 extrudeStatus[f[fp]] == NOEXTRUDE
1403 && grownExtrudeStatus[f[fp]] != NOEXTRUDE
1406 grownExtrudeStatus[f[fp]] = NOEXTRUDE;
1413 extrudeStatus.transfer(grownExtrudeStatus);
1415 forAll(extrudeStatus, patchPointI)
1417 if (extrudeStatus[patchPointI] == NOEXTRUDE)
1420 patchNLayers[patchPointI] = 0;
1424 reduce(nGrown, sumOp<label>());
1426 Info<<
"Set displacement to zero for an additional " << nGrown
1427 <<
" points." <<
endl;
1431 void Foam::autoLayerDriver::calculateLayerThickness
1437 const bool relativeSizes,
1443 const scalar edge0Len,
1450 const fvMesh& mesh = meshRefiner_.mesh();
1451 const polyBoundaryMesh& patches = mesh.boundaryMesh();
1458 expansionRatio.setSize(pp.nPoints());
1459 expansionRatio = GREAT;
1460 thickness.setSize(pp.nPoints());
1462 minThickness.setSize(pp.nPoints());
1463 minThickness = GREAT;
1467 label patchI = patchIDs[i];
1469 const labelList& meshPoints = patches[patchI].meshPoints();
1471 forAll(meshPoints, patchPointI)
1473 label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
1475 expansionRatio[ppPointI] =
min
1477 expansionRatio[ppPointI],
1478 patchExpansionRatio[patchI]
1480 thickness[ppPointI] =
min
1482 thickness[ppPointI],
1483 patchFinalLayerThickness[patchI]
1485 minThickness[ppPointI] =
min
1487 minThickness[ppPointI],
1488 patchMinThickness[patchI]
1532 if (
min(patchMinThickness) < 0 ||
max(patchMinThickness) > 2)
1535 <<
"Thickness should be factor of local undistorted cell size."
1536 <<
" Valid values are [0..2]." << nl
1537 <<
" minThickness:" << patchMinThickness
1545 labelList maxPointLevel(pp.nPoints(), labelMin);
1549 label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];
1551 const face& f = pp.localFaces()[i];
1555 maxPointLevel[f[fp]] =
max(maxPointLevel[f[fp]], ownLevel);
1570 forAll(maxPointLevel, pointI)
1573 scalar edgeLen = edge0Len/(1<<maxPointLevel[pointI]);
1574 thickness[pointI] *= edgeLen;
1575 minThickness[pointI] *= edgeLen;
1584 forAll(thickness, pointI)
1588 if (expansionRatio[pointI] == 1)
1590 thickness[pointI] *= patchNLayers[pointI];
1595 scalar invExpansion = 1.0 / expansionRatio[pointI];
1596 label nLay = patchNLayers[pointI];
1597 thickness[pointI] *=
1598 (1.0 -
pow(invExpansion, nLay))
1599 / (1.0 - invExpansion);
1610 void Foam::autoLayerDriver::syncPatchDisplacement
1612 const motionSmoother& meshMover,
1616 List<extrudeMode>& extrudeStatus
1619 const fvMesh& mesh = meshRefiner_.mesh();
1620 const labelList& meshPoints = meshMover.patch().meshPoints();
1622 label nChangedTotal = 0;
1642 if (
mag(patchDisp[i]) < minThickness[i])
1660 labelList syncPatchNLayers(patchNLayers);
1673 forAll(syncPatchNLayers, i)
1675 if (syncPatchNLayers[i] != patchNLayers[i])
1704 forAll(syncPatchNLayers, i)
1706 if (syncPatchNLayers[i] != patchNLayers[i])
1723 nChangedTotal += nChanged;
1731 Info<<
"Prevented extrusion on "
1733 <<
" coupled patch points during syncPatchDisplacement." <<
endl;
1742 void Foam::autoLayerDriver::getPatchDisplacement
1744 const motionSmoother& meshMover,
1749 List<extrudeMode>& extrudeStatus
1752 Info<< nl <<
"Determining displacement for added points"
1753 <<
" according to pointNormal ..." <<
endl;
1755 const fvMesh& mesh = meshRefiner_.mesh();
1757 const vectorField& faceNormals = pp.faceNormals();
1759 const pointField& localPoints = pp.localPoints();
1760 const labelList& meshPoints = pp.meshPoints();
1769 forAll(faceNormals, faceI)
1771 const face& f = pp.localFaces()[faceI];
1775 pointNormals[f[fp]] += faceNormals[faceI];
1776 nPointFaces[f[fp]] ++;
1802 pointNormals[i] /= nPointFaces[i];
1811 patchDisp = thickness*pointNormals;
1814 forAll(pointNormals, patchPointI)
1816 label meshPointI = pp.meshPoints()[patchPointI];
1818 if (extrudeStatus[patchPointI] == NOEXTRUDE)
1821 patchNLayers[patchPointI] = 0;
1827 const vector& n = pointNormals[patchPointI];
1831 Pout<<
"No valid normal for point " << meshPointI
1832 <<
' ' << pp.points()[meshPointI]
1833 <<
"; setting displacement to " << patchDisp[patchPointI]
1836 extrudeStatus[patchPointI] = EXTRUDEREMOVE;
1842 forAll(extrudeStatus, patchPointI)
1844 if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
1849 const labelList& pEdges = pp.pointEdges()[patchPointI];
1853 label edgeI = pEdges[i];
1855 label otherPointI = pp.edges()[edgeI].otherVertex(patchPointI);
1857 if (extrudeStatus[otherPointI] != NOEXTRUDE)
1859 avg += localPoints[otherPointI] + patchDisp[otherPointI];
1866 Pout<<
"Displacement at illegal point "
1867 << localPoints[patchPointI]
1868 <<
" set to " << (avg / nPoints - localPoints[patchPointI])
1871 patchDisp[patchPointI] =
1873 - localPoints[patchPointI];
1879 syncPatchDisplacement
1895 Foam::label Foam::autoLayerDriver::truncateDisplacement
1897 const motionSmoother& meshMover,
1899 const faceSet& illegalPatchFaces,
1902 List<extrudeMode>& extrudeStatus
1905 const polyMesh& mesh = meshMover.mesh();
1910 const Map<label>& meshPointMap = pp.meshPointMap();
1914 label faceI = iter.key();
1916 if (mesh.isInternalFace(faceI))
1919 <<
"Faceset " << illegalPatchFaces.name()
1920 <<
" contains internal face " << faceI << nl
1924 const face& f = mesh.faces()[faceI];
1929 if (meshPointMap.found(f[fp]))
1931 label patchPointI = meshPointMap[f[fp]];
1933 if (extrudeStatus[patchPointI] != NOEXTRUDE)
1948 forAll(patchDisp, patchPointI)
1950 if (
mag(patchDisp[patchPointI]) < minThickness[patchPointI])
1966 else if (extrudeStatus[patchPointI] == NOEXTRUDE)
1970 patchNLayers[patchPointI] = 0;
1975 const faceList& localFaces = pp.localFaces();
1979 syncPatchDisplacement
1996 const face& localF = localFaces[i];
2001 extrudeMode prevMode = extrudeStatus[localF.prevLabel(0)];
2005 extrudeMode fpMode = extrudeStatus[localF[fp]];
2007 if (prevMode == NOEXTRUDE && fpMode != NOEXTRUDE)
2034 reduce(nPinched, sumOp<label>());
2036 Info<<
"truncateDisplacement : Unextruded " << nPinched
2037 <<
" faces due to non-consecutive vertices being extruded." <<
endl;
2043 label nDiffering = 0;
2087 if (nPinched+nDiffering == 0)
2099 void Foam::autoLayerDriver::setupLayerInfoTruncation
2101 const motionSmoother& meshMover,
2103 const List<extrudeMode>& extrudeStatus,
2104 const label nBufferCellsNoExtrude,
2109 Info<< nl <<
"Setting up information for layer truncation ..." <<
endl;
2112 const polyMesh& mesh = meshMover.mesh();
2114 if (nBufferCellsNoExtrude < 0)
2116 Info<< nl <<
"Performing no layer truncation."
2117 <<
" nBufferCellsNoExtrude set to less than 0 ..." <<
endl;
2120 forAll(pp.localFaces(), patchFaceI)
2122 const face& f = pp.localFaces()[patchFaceI];
2126 if (patchNLayers[f[fp]] > 0)
2128 nPatchFaceLayers[patchFaceI] = patchNLayers[f[fp]];
2133 nPatchPointLayers = patchNLayers;
2140 forAll(pp.localFaces(), patchFaceI)
2142 const face& f = pp.localFaces()[patchFaceI];
2147 bool noExtrude =
false;
2152 if (extrudeStatus[f[fp]] == NOEXTRUDE)
2156 mLevel =
max(mLevel, patchNLayers[f[fp]]);
2166 nPatchFaceLayers[patchFaceI] = 1;
2167 maxLevel[patchFaceI] = mLevel;
2171 maxLevel[patchFaceI] = mLevel;
2183 label nLevels =
gMax(patchNLayers);
2186 for (label ilevel = 1; ilevel < nLevels; ilevel++)
2192 nBuffer = nBufferCellsNoExtrude - 1;
2196 nBuffer = nBufferCellsNoExtrude;
2199 for (label ibuffer = 0; ibuffer < nBuffer + 1; ibuffer++)
2201 labelList tempCounter(nPatchFaceLayers);
2203 boolList foundNeighbour(pp.nPoints(),
false);
2205 forAll(pp.meshPoints(), patchPointI)
2207 forAll(pointFaces[patchPointI], pointFaceI)
2209 label faceI = pointFaces[patchPointI][pointFaceI];
2213 nPatchFaceLayers[faceI] != -1
2214 && maxLevel[faceI] > 0
2217 foundNeighbour[patchPointI] =
true;
2233 forAll(pp.meshPoints(), patchPointI)
2235 if (foundNeighbour[patchPointI])
2237 forAll(pointFaces[patchPointI], pointFaceI)
2239 label faceI = pointFaces[patchPointI][pointFaceI];
2242 nPatchFaceLayers[faceI] == -1
2243 && maxLevel[faceI] > 0
2244 && ilevel < maxLevel[faceI]
2247 tempCounter[faceI] = ilevel;
2252 nPatchFaceLayers = tempCounter;
2256 forAll(pp.localFaces(), patchFaceI)
2258 if (nPatchFaceLayers[patchFaceI] == -1)
2260 nPatchFaceLayers[patchFaceI] = maxLevel[patchFaceI];
2264 forAll(pp.meshPoints(), patchPointI)
2266 if (extrudeStatus[patchPointI] != NOEXTRUDE)
2268 forAll(pointFaces[patchPointI], pointFaceI)
2270 label face = pointFaces[patchPointI][pointFaceI];
2271 nPatchPointLayers[patchPointI] =
max
2273 nPatchPointLayers[patchPointI],
2274 nPatchFaceLayers[face]
2280 nPatchPointLayers[patchPointI] = 0;
2297 bool Foam::autoLayerDriver::cellsUseFace
2299 const polyMesh& mesh,
2306 const cell& cFaces = mesh.cells()[cellLabels[i]];
2310 if (faces.found(cFaces[cFaceI]))
2323 Foam::label Foam::autoLayerDriver::checkAndUnmark
2325 const addPatchCellLayer& addLayer,
2326 const dictionary& meshQualityDict,
2328 const fvMesh& newMesh,
2332 List<extrudeMode>& extrudeStatus
2336 Info<< nl <<
"Checking mesh with layer ..." <<
endl;
2337 faceSet wrongFaces(newMesh,
"wrongFaces", newMesh.nFaces()/1000);
2341 <<
" (concave, zero area or negative cell pyramid volume)"
2355 addLayer.layerFaces()
2360 forAll(addedCells, oldPatchFaceI)
2364 const labelList& fCells = addedCells[oldPatchFaceI];
2366 if (cellsUseFace(newMesh, fCells, wrongFaces))
2373 pp.localFaces()[oldPatchFaceI],
2390 Foam::label Foam::autoLayerDriver::countExtrusion
2393 const List<extrudeMode>& extrudeStatus
2397 label nExtruded = 0;
2399 const faceList& localFaces = pp.localFaces();
2403 const face& localFace = localFaces[i];
2407 if (extrudeStatus[localFace[fp]] != NOEXTRUDE)
2421 void Foam::autoLayerDriver::getLayerCellsFaces
2423 const polyMesh& mesh,
2424 const addPatchCellLayer& addLayer,
2429 flaggedCells.setSize(mesh.nCells());
2430 flaggedCells =
false;
2431 flaggedFaces.setSize(mesh.nFaces());
2432 flaggedFaces =
false;
2440 forAll(addedCells, oldPatchFaceI)
2442 const labelList& added = addedCells[oldPatchFaceI];
2446 flaggedCells[added[i]] =
true;
2450 forAll(layerFaces, oldPatchFaceI)
2452 const labelList& layer = layerFaces[oldPatchFaceI];
2456 for (label i = 1; i < layer.size()-1; i++)
2458 flaggedFaces[layer[i]] =
true;
2469 meshRefiner_(meshRefiner)
2475 void Foam::autoLayerDriver::mergePatchFacesUndo
2494 <<
"Merging all faces of a cell" << nl
2495 <<
"---------------------------" << nl
2496 <<
" - which are on the same patch" << nl
2497 <<
" - which make an angle < " << layerParams.
featureAngle()
2500 <<
" (cos:" << minCos <<
')' << nl
2501 <<
" - as long as the resulting face doesn't become concave"
2504 <<
" (0=straight, 180=fully concave)" << nl
2507 label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict);
2509 nChanged += mergeEdgesUndo(minCos, motionDict);
2518 const label nAllowableErrors,
2523 fvMesh& mesh = meshRefiner_.mesh();
2535 Info<<
"Constructing mesh displacer ..." <<
endl;
2574 meshMover().adaptPatchIDs(),
2616 const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
2617 const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
2645 for (label i = 0; i < layerParams.
nGrow(); i++)
2659 const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
2660 const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
2666 calculateLayerThickness
2669 meshMover().adaptPatchIDs(),
2691 label maxPatchNameLen = 0;
2692 forAll(meshMover().adaptPatchIDs(), i)
2694 label patchI = meshMover().adaptPatchIDs()[i];
2695 word patchName = patches[patchI].
name();
2696 maxPatchNameLen =
max(maxPatchNameLen, label(patchName.size()));
2700 <<
setf(ios_base::left) <<
setw(maxPatchNameLen) <<
"patch"
2701 <<
setw(0) <<
" faces layers avg thickness[m]" << nl
2702 <<
setf(ios_base::left) <<
setw(maxPatchNameLen) <<
" "
2703 <<
setw(0) <<
" near-wall overall" << nl
2704 <<
setf(ios_base::left) <<
setw(maxPatchNameLen) <<
"-----"
2705 <<
setw(0) <<
" ----- ------ --------- -------" <<
endl;
2707 forAll(meshMover().adaptPatchIDs(), i)
2709 label patchI = meshMover().adaptPatchIDs()[i];
2711 const labelList& meshPoints = patches[patchI].meshPoints();
2715 scalar sumThickness = 0;
2716 scalar sumNearWallThickness = 0;
2718 forAll(meshPoints, patchPointI)
2720 label ppPointI = pp().meshPointMap()[meshPoints[patchPointI]];
2724 sumThickness += thickness[ppPointI];
2726 label nLay = patchNLayers[ppPointI];
2729 if (expansionRatio[ppPointI] == 1)
2731 sumNearWallThickness += thickness[ppPointI]/nLay;
2736 (1.0-
pow(expansionRatio[ppPointI], nLay))
2737 / (1.0-expansionRatio[ppPointI]);
2738 sumNearWallThickness += thickness[ppPointI]/s;
2746 scalar avgThickness = 0;
2747 scalar avgNearWallThickness = 0;
2756 avgNearWallThickness =
2766 <<
" " <<
setw(8) << avgNearWallThickness
2767 <<
" " <<
setw(8) << avgThickness
2784 meshRefiner_.timeName(),
2790 meshMover().pMesh(),
2799 meshRefiner_.timeName(),
2805 meshMover().pMesh(),
2814 meshRefiner_.timeName(),
2820 meshMover().pMesh(),
2830 medialAxisSmoothingInfo
2853 for (label iteration = 0; iteration < layerParams.
nLayerIter(); iteration++)
2856 <<
"Layer addition iteration " << iteration << nl
2857 <<
"--------------------------" <<
endl;
2865 : motionDict.
subDict(
"relaxed")
2870 Info<<
"Switched to relaxed meshQuality constraints." <<
endl;
2876 syncPatchDisplacement
2886 getPatchDisplacement
2913 shrinkMeshMedialDistance
2921 layerParams.
nSnap(),
2937 patchDisp = oldPatchPos - pp().localPoints();
2942 faceSet dummySet(mesh,
"wrongPatchFaces", 0);
2943 truncateDisplacement
2966 Info<<
"Writing shrunk mesh to " << meshRefiner_.timeName() <<
endl;
2985 labelList nPatchFaceLayers(pp().localFaces().size(),-1);
2986 setupLayerInfoTruncation
3002 if (patchNLayers[i] > 0)
3004 if (expansionRatio[i] == 1.0)
3006 firstDisp[i] = patchDisp[i]/nPatchPointLayers[i];
3010 label nLay = nPatchPointLayers[i];
3012 pow(expansionRatio[i], nLay - 1)
3013 * (1.0 - expansionRatio[i])
3014 / (1.0 -
pow(expansionRatio[i], nLay));
3015 firstDisp[i] = h*patchDisp[i];
3020 scalarField invExpansionRatio = 1.0 / expansionRatio;
3040 savedMeshMod = meshMod;
3062 fvMesh& newMesh = newMeshPtr();
3067 if (meshRefiner_.overwrite())
3094 Info<<
"Writing layer mesh to " << meshRefiner_.timeName() <<
endl;
3104 <<
" added cells to cellSet "
3106 addedCellSet.
write();
3116 <<
" faces inside added layer to faceSet "
3118 layerFacesSet.
write();
3122 label nTotChanged = checkAndUnmark
3134 Info<<
"Extruding " << countExtrusion(pp, extrudeStatus)
3136 <<
" faces. Removed extrusion at " << nTotChanged <<
" faces."
3139 if (nTotChanged == 0)
3145 meshMover().movePoints(oldPoints);
3146 meshMover().correct();
3163 if (map().hasMotionPoints())
3173 if (meshRefiner_.overwrite())
3178 meshRefiner_.updateMesh(map,
labelList(0));
3187 <<
"Doing final balancing" << nl
3188 <<
"---------------------" << nl
3207 map().distributeCellData(flaggedCells);
3208 map().distributeFaceData(flaggedFaces);
3218 <<
" added cells to cellSet "
3220 addedCellSet.
write();
3225 <<
" faces inside added layer to faceSet "
3227 layerFacesSet.
write();
3236 const bool preBalance,
3241 const fvMesh& mesh = meshRefiner_.mesh();
3244 <<
"Shrinking and layer addition phase" << nl
3245 <<
"----------------------------------" << nl
3248 Info<<
"Using mesh parameters " << motionDict << nl <<
endl;
3251 mergePatchFacesUndo(layerParams, motionDict);
3258 label nFacesWithLayers = 0;
3259 forAll(numLayers, patchI)
3261 if (numLayers[patchI] > 0)
3271 Info<< nl <<
"No layers to generate ..." <<
endl;
3276 checkMeshManifold();
3279 Info<<
"Checking initial mesh ..." <<
endl;
3288 Info<<
"Detected " << nInitErrors <<
" illegal faces"
3289 <<
" (concave, zero area or negative cell pyramid volume)"
3297 <<
"Doing initial balancing" << nl
3298 <<
"-----------------------" << nl
3302 forAll(numLayers, patchI)
3304 if (numLayers[patchI] > 0)
3309 cellWeights[pp.
faceCells()[i]] += numLayers[patchI];