49 void Foam::meshRefinement::markBoundaryFace
57 isBoundaryFace[facei] =
true;
61 for (
const label edgei : fEdges)
63 isBoundaryEdge[edgei] =
true;
66 const face&
f = mesh_.
faces()[facei];
68 for (
const label pointi :
f)
70 isBoundaryPoint[pointi] =
true;
75 void Foam::meshRefinement::findNearest
78 List<pointIndexHit>& nearestInfo,
87 fc[i] = mesh_.faceCentres()[meshFaces[i]];
102 nearestNormal.setSize(nearestInfo.size());
103 nearestRegion.setSize(nearestInfo.size());
105 forAll(allSurfaces, surfI)
107 DynamicList<pointIndexHit> localHits;
111 if (nearestSurface[i] == surfI)
113 localHits.append(nearestInfo[i]);
117 label geomI = surfaces_.surfaces()[surfI];
120 surfaces_.geometry()[geomI].getNormal(localHits, localNormals);
123 surfaces_.geometry()[geomI].getRegion(localHits, localRegion);
128 if (nearestSurface[i] == surfI)
130 nearestNormal[i] = localNormals[localI];
131 nearestRegion[i] = localRegion[localI];
146 autoPtr<indirectPrimitivePatch> ppPtr
160 DynamicList<label> candidateFaces(pp.size()/20);
164 const labelList& cellLevel = meshCutter_.cellLevel();
168 const labelList& eFaces = edgeFaces[edgeI];
170 if (eFaces.size() == 2)
172 label face0 = pp.addressing()[eFaces[0]];
173 label face1 = pp.addressing()[eFaces[1]];
175 label cell0 = mesh_.faceOwner()[face0];
176 label cell1 = mesh_.faceOwner()[face1];
178 if (cellLevel[cell0] > cellLevel[cell1])
181 const vector& n0 = pp.faceNormals()[eFaces[0]];
182 const vector& n1 = pp.faceNormals()[eFaces[1]];
184 if (
mag(n0 & n1) < 0.1)
186 candidateFaces.append(face0);
189 else if (cellLevel[cell1] > cellLevel[cell0])
192 const vector& n0 = pp.faceNormals()[eFaces[0]];
193 const vector& n1 = pp.faceNormals()[eFaces[1]];
195 if (
mag(n0 & n1) < 0.1)
197 candidateFaces.append(face1);
202 candidateFaces.shrink();
205 <<
" faces on edge-connected cells of differing level."
210 faceSet fSet(mesh_,
"edgeConnectedFaces", candidateFaces);
212 Pout<<
"Writing " << fSet.size()
213 <<
" with problematic topology to faceSet "
214 << fSet.objectPath() <<
endl;
222 List<pointIndexHit> nearestInfo;
239 Map<label> candidateCells(candidateFaces.size());
241 faceSet perpFaces(mesh_,
"perpendicularFaces", pp.size()/100);
245 label facei = candidateFaces[i];
249 label region = surfaces_.globalRegion
255 scalar angle =
degToRad(perpendicularAngle[region]);
261 perpFaces.insert(facei);
262 candidateCells.insert
264 mesh_.faceOwner()[facei],
265 globalToPatch[region]
274 Pout<<
"Writing " << perpFaces.size()
275 <<
" faces that are perpendicular to the surface to set "
276 << perpFaces.objectPath() <<
endl;
279 return candidateCells;
286 bool Foam::meshRefinement::isCollapsedFace
290 const scalar minFaceArea,
291 const scalar maxNonOrtho,
296 const scalar severeNonorthogonalityThreshold =
300 scalar magS =
mag(
s);
303 if (magS < minFaceArea)
309 const point& ownCc = mesh_.cellCentres()[mesh_.faceOwner()[facei]];
311 if (mesh_.isInternalFace(facei))
313 label nei = mesh_.faceNeighbour()[facei];
314 vector d = mesh_.cellCentres()[nei] - ownCc;
316 scalar dDotS = (d &
s)/(
mag(d)*magS + VSMALL);
318 if (dDotS < severeNonorthogonalityThreshold)
329 label patchi = mesh_.boundaryMesh().whichPatch(facei);
331 if (mesh_.boundaryMesh()[patchi].coupled())
333 vector d = neiCc[facei-mesh_.nInternalFaces()] - ownCc;
335 scalar dDotS = (d &
s)/(
mag(d)*magS + VSMALL);
337 if (dDotS < severeNonorthogonalityThreshold)
357 bool Foam::meshRefinement::isCollapsedCell
360 const scalar volFraction,
364 scalar vol = mesh_.cells()[celli].mag(
points, mesh_.faces());
366 if (vol/mesh_.cellVolumes()[celli] < volFraction)
385 const dictionary& motionDict,
386 const bool removeEdgeConnectedCells,
391 const labelList& cellLevel = meshCutter_.cellLevel();
392 const labelList& pointLevel = meshCutter_.pointLevel();
393 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
398 boolList isBoundaryPoint(mesh_.nPoints(),
false);
399 boolList isBoundaryEdge(mesh_.nEdges(),
false);
400 boolList isBoundaryFace(mesh_.nFaces(),
false);
404 labelList adaptPatchIDs(meshedPatches());
408 const polyPatch& pp =
patches[adaptPatchIDs[i]];
410 label facei = pp.
start();
428 const labelList nearestAdaptPatch(nearestPatch(adaptPatchIDs));
436 labelList neiLevel(mesh_.nBoundaryFaces());
438 calcNeighbourData(neiLevel, neiCc);
442 label nBaffleFaces = 0;
446 label nPrevented = 0;
448 if (removeEdgeConnectedCells &&
max(perpendicularAngle) >= 0)
450 Info<<
"markFacesOnProblemCells :"
451 <<
" Checking for edge-connected cells of highly differing sizes."
456 Map<label> problemCells
458 findEdgeConnectedProblemCells
468 const cell& cFaces = mesh_.cells()[iter.key()];
472 label facei = cFaces[i];
474 if (facePatch[facei] == -1 && mesh_.isInternalFace(facei))
476 facePatch[facei] = nearestAdaptPatch[facei];
490 Info<<
"markFacesOnProblemCells : Marked "
492 <<
" additional internal faces to be converted into baffles"
495 <<
" cells edge-connected to lower level cells." <<
endl;
499 cellSet problemCellSet(mesh_,
"problemCells", problemCells.toc());
500 problemCellSet.instance() =
timeName();
501 Pout<<
"Writing " << problemCellSet.size()
502 <<
" cells that are edge connected to coarser cell to set "
503 << problemCellSet.objectPath() <<
endl;
504 problemCellSet.
write();
536 const scalar volFraction =
537 motionDict.getOrDefault<scalar>(
"minVolCollapseRatio", -1);
539 const bool checkCollapse = (volFraction > 0);
541 scalar maxNonOrtho = -1;
549 minArea = get<scalar>(motionDict,
"minArea", dryRun_);
550 maxNonOrtho = get<scalar>(motionDict,
"maxNonOrtho", dryRun_);
552 Info<<
"markFacesOnProblemCells :"
553 <<
" Deleting all-anchor surface cells only if"
554 <<
" snapping them violates mesh quality constraints:" <<
nl
555 <<
" snapped/original cell volume < " << volFraction <<
nl
556 <<
" face area < " << minArea <<
nl
557 <<
" non-orthogonality > " << maxNonOrtho <<
nl
561 autoPtr<indirectPrimitivePatch> ppPtr
570 const pointField& localPoints = pp.localPoints();
571 const labelList& meshPoints = pp.meshPoints();
573 List<pointIndexHit> hitInfo;
575 surfaces_.findNearest
585 newPoints = mesh_.points();
589 if (hitInfo[i].hit())
591 newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
597 const_cast<Time&
>(mesh_.time())++;
599 mesh_.movePoints(newPoints);
605 writeType(writeLevel() | WRITEMESH),
606 mesh_.time().path()/
"newPoints"
608 mesh_.movePoints(oldPoints);
626 bitSet hasSevenBoundaryAnchorPoints(mesh_.nCells());
631 DynamicList<label> dynFEdges;
632 DynamicList<label> dynCPoints;
637 const labelList& cPoints = mesh_.cellPoints(celli, pSet, dynCPoints);
641 label nBoundaryAnchors = 0;
642 label nNonAnchorBoundary = 0;
643 label nonBoundaryAnchor = -1;
645 for (
const label pointi : cPoints)
647 if (pointLevel[pointi] <= cellLevel[celli])
650 if (isBoundaryPoint[pointi])
657 nonBoundaryAnchor = pointi;
660 else if (isBoundaryPoint[pointi])
662 nNonAnchorBoundary++;
666 if (nBoundaryAnchors == 8)
668 const cell& cFaces = mesh_.cells()[celli];
675 if (isBoundaryFace[cFaces[cFacei]])
691 && !isCollapsedCell(newPoints, volFraction, celli)
708 label facei = cFaces[cf];
712 facePatch[facei] == -1
713 && mesh_.isInternalFace(facei)
716 facePatch[facei] = nearestAdaptPatch[facei];
732 else if (nBoundaryAnchors == 7 && nonBoundaryAnchor != -1)
735 hasSevenBoundaryAnchorPoints.set(celli);
736 nonBoundaryAnchors.insert(nonBoundaryAnchor);
745 DynamicList<label> dynPCells;
747 for (
const label pointi : nonBoundaryAnchors)
749 const labelList& pCells = mesh_.pointCells(pointi, dynPCells);
754 for (
const label celli : pCells)
756 if (hasSevenBoundaryAnchorPoints.test(celli))
765 for (
const label celli : pCells)
767 if (hasSevenBoundaryAnchorPoints.test(celli))
772 && !isCollapsedCell(newPoints, volFraction, celli)
787 const cell& cFaces = mesh_.cells()[celli];
789 for (
const label facei : cFaces)
793 facePatch[facei] == -1
794 && mesh_.isInternalFace(facei)
797 facePatch[facei] = nearestAdaptPatch[facei];
845 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
847 if (facePatch[facei] == -1)
849 const labelList& fEdges = mesh_.faceEdges(facei, dynFEdges);
850 label nFaceBoundaryEdges = 0;
852 for (
const label edgei : fEdges)
854 if (isBoundaryEdge[edgei])
856 ++nFaceBoundaryEdges;
860 if (nFaceBoundaryEdges == fEdges.size())
884 facePatch[facei] = nearestAdaptPatch[facei];
894 for (
const polyPatch& pp :
patches)
898 label facei = pp.start();
902 if (facePatch[facei] == -1)
904 const labelList& fEdges = mesh_.faceEdges(facei, dynFEdges);
905 label nFaceBoundaryEdges = 0;
907 for (
const label edgei : fEdges)
909 if (isBoundaryEdge[edgei])
911 ++nFaceBoundaryEdges;
915 if (nFaceBoundaryEdges == fEdges.size())
939 facePatch[facei] = nearestAdaptPatch[facei];
940 if (isMasterFace.test(facei))
969 Info<<
"markFacesOnProblemCells : marked "
971 <<
" additional internal faces to be converted into baffles."
976 Info<<
"markFacesOnProblemCells : prevented "
978 <<
" internal faces from getting converted into baffles."
990 const snapParameters& snapParams,
991 const dictionary& motionDict,
1000 labelList adaptPatchIDs(meshedPatches());
1003 autoPtr<indirectPrimitivePatch> ppPtr
1021 Info<<
"Constructing mesh displacer ..." <<
endl;
1022 Info<<
"Using mesh parameters " << motionDict <<
nl <<
endl;
1026 motionSmoother meshMover
1037 Info<<
"Checking initial mesh ..." <<
endl;
1053 Info<<
"Detected " << nInitErrors <<
" illegal faces"
1054 <<
" (concave, zero area or negative cell pyramid volume)"
1058 Info<<
"Checked initial mesh in = "
1059 << mesh_.time().cpuTimeIncrement() <<
" s\n" <<
nl <<
endl;
1075 snappySnapDriver::calcNearestSurface
1077 snapParams.strictRegionSnap(),
1079 globalToMasterPatch,
1088 const labelList& meshPoints = pp.meshPoints();
1093 newPoints[meshPoints[i]] += disp[i];
1100 minMagSqrEqOp<point>(),
1101 vector(GREAT, GREAT, GREAT)
1104 mesh_.movePoints(newPoints);
1109 const labelList nearestAdaptPatch(nearestPatch(meshedPatches()));
1113 labelList facePatch(mesh_.nFaces(), -1);
1115 label nBaffleFaces = 0;
1118 faceSet wrongFaces(mesh_,
"wrongFaces", 100);
1126 label nWrongFaces = 0;
1158 scalar minArea(get<scalar>(motionDict,
"minArea", dryRun_));
1159 if (minArea > -SMALL)
1177 Info<<
" faces with area < "
1178 <<
setw(5) << minArea
1180 << nNewWrongFaces-nWrongFaces <<
endl;
1182 nWrongFaces = nNewWrongFaces;
1185 scalar minDet(get<scalar>(motionDict,
"minDeterminant", dryRun_));
1205 Info<<
" faces on cells with determinant < "
1206 <<
setw(5) << minDet <<
" : "
1207 << nNewWrongFaces-nWrongFaces <<
endl;
1209 nWrongFaces = nNewWrongFaces;
1214 for (
const label facei : wrongFaces)
1216 const label patchi = mesh_.boundaryMesh().whichPatch(facei);
1218 if (patchi == -1 || mesh_.boundaryMesh()[patchi].coupled())
1220 facePatch[facei] = nearestAdaptPatch[facei];
1233 mesh_.movePoints(oldPoints);
1236 Info<<
"markFacesOnProblemCellsGeometric : marked "
1238 <<
" additional internal and coupled faces"
1239 <<
" to be converted into baffles." <<
endl;