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)
383 void Foam::meshRefinement::markFacesOnProblemCells
385 const dictionary& motionDict,
386 const bool removeEdgeConnectedCells,
394 const labelList& cellLevel = meshCutter_.cellLevel();
395 const labelList& pointLevel = meshCutter_.pointLevel();
396 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
401 boolList isBoundaryPoint(mesh_.nPoints(),
false);
402 boolList isBoundaryEdge(mesh_.nEdges(),
false);
403 boolList isBoundaryFace(mesh_.nFaces(),
false);
407 const labelList adaptPatchIDs(meshedPatches());
411 const polyPatch& pp =
patches[adaptPatchIDs[i]];
413 label facei = pp.
start();
433 nearestPatch(adaptPatchIDs, nearestAdaptPatch, nearestAdaptZone);
438 facePatch.setSize(mesh_.nFaces());
440 faceZone.setSize(mesh_.nFaces());
444 labelList neiLevel(mesh_.nBoundaryFaces());
446 calcNeighbourData(neiLevel, neiCc);
450 label nBaffleFaces = 0;
454 label nPrevented = 0;
456 if (removeEdgeConnectedCells &&
max(perpendicularAngle) >= 0)
458 Info<<
"markFacesOnProblemCells :" 459 <<
" Checking for edge-connected cells of highly differing sizes." 464 Map<label> problemCells
466 findEdgeConnectedProblemCells
476 const cell& cFaces = mesh_.cells()[iter.key()];
480 label facei = cFaces[i];
482 if (facePatch[facei] == -1 && mesh_.isInternalFace(facei))
484 facePatch[facei] = nearestAdaptPatch[facei];
485 faceZone[facei] = nearestAdaptZone[facei];
499 Info<<
"markFacesOnProblemCells : Marked " 501 <<
" additional internal faces to be converted into baffles" 504 <<
" cells edge-connected to lower level cells." <<
endl;
508 cellSet problemCellSet(mesh_,
"problemCells", problemCells.toc());
509 problemCellSet.instance() =
timeName();
510 Pout<<
"Writing " << problemCellSet.size()
511 <<
" cells that are edge connected to coarser cell to set " 512 << problemCellSet.objectPath() <<
endl;
513 problemCellSet.
write();
545 const scalar volFraction =
546 motionDict.getOrDefault<scalar>(
"minVolCollapseRatio", -1);
548 const bool checkCollapse = (volFraction > 0);
550 scalar maxNonOrtho = -1;
558 minArea = get<scalar>(motionDict,
"minArea", dryRun_);
559 maxNonOrtho = get<scalar>(motionDict,
"maxNonOrtho", dryRun_);
561 Info<<
"markFacesOnProblemCells :" 562 <<
" Deleting all-anchor surface cells only if" 563 <<
" snapping them violates mesh quality constraints:" <<
nl 564 <<
" snapped/original cell volume < " << volFraction <<
nl 565 <<
" face area < " << minArea <<
nl 566 <<
" non-orthogonality > " << maxNonOrtho <<
nl 570 autoPtr<indirectPrimitivePatch> ppPtr
579 const pointField& localPoints = pp.localPoints();
580 const labelList& meshPoints = pp.meshPoints();
582 List<pointIndexHit> hitInfo;
584 surfaces_.findNearest
594 newPoints = mesh_.points();
598 if (hitInfo[i].hit())
600 newPoints[meshPoints[i]] = hitInfo[i].point();
606 const_cast<Time&
>(mesh_.time())++;
608 mesh_.movePoints(newPoints);
618 writeType(writeLevel() | WRITEMESH),
619 mesh_.time().path()/
"newPoints" 621 mesh_.movePoints(oldPoints);
642 bitSet hasSevenBoundaryAnchorPoints(mesh_.nCells());
647 DynamicList<label> dynFEdges;
648 DynamicList<label> dynCPoints;
653 const labelList& cPoints = mesh_.cellPoints(celli, pSet, dynCPoints);
657 label nBoundaryAnchors = 0;
658 label nNonAnchorBoundary = 0;
659 label nonBoundaryAnchor = -1;
661 for (
const label pointi : cPoints)
663 if (pointLevel[pointi] <= cellLevel[celli])
666 if (isBoundaryPoint[pointi])
673 nonBoundaryAnchor = pointi;
676 else if (isBoundaryPoint[pointi])
678 nNonAnchorBoundary++;
682 if (nBoundaryAnchors == 8)
684 const cell& cFaces = mesh_.cells()[celli];
691 if (isBoundaryFace[cFaces[cFacei]])
707 && !isCollapsedCell(newPoints, volFraction, celli)
724 label facei = cFaces[cf];
728 facePatch[facei] == -1
729 && mesh_.isInternalFace(facei)
732 facePatch[facei] = nearestAdaptPatch[facei];
733 faceZone[facei] = nearestAdaptZone[facei];
749 else if (nBoundaryAnchors == 7 && nonBoundaryAnchor != -1)
752 hasSevenBoundaryAnchorPoints.set(celli);
753 nonBoundaryAnchors.insert(nonBoundaryAnchor);
762 DynamicList<label> dynPCells;
764 for (
const label pointi : nonBoundaryAnchors)
766 const labelList& pCells = mesh_.pointCells(pointi, dynPCells);
771 for (
const label celli : pCells)
773 if (hasSevenBoundaryAnchorPoints.test(celli))
782 for (
const label celli : pCells)
784 if (hasSevenBoundaryAnchorPoints.test(celli))
789 && !isCollapsedCell(newPoints, volFraction, celli)
804 const cell& cFaces = mesh_.cells()[celli];
806 for (
const label facei : cFaces)
810 facePatch[facei] == -1
811 && mesh_.isInternalFace(facei)
814 facePatch[facei] = nearestAdaptPatch[facei];
815 faceZone[facei] = nearestAdaptZone[facei];
863 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
865 if (facePatch[facei] == -1)
867 const labelList& fEdges = mesh_.faceEdges(facei, dynFEdges);
868 label nFaceBoundaryEdges = 0;
870 for (
const label edgei : fEdges)
872 if (isBoundaryEdge[edgei])
874 ++nFaceBoundaryEdges;
878 if (nFaceBoundaryEdges == fEdges.size())
902 facePatch[facei] = nearestAdaptPatch[facei];
903 faceZone[facei] = nearestAdaptZone[facei];
913 for (
const polyPatch& pp :
patches)
917 label facei = pp.start();
921 if (facePatch[facei] == -1)
923 const labelList& fEdges = mesh_.faceEdges(facei, dynFEdges);
924 label nFaceBoundaryEdges = 0;
926 for (
const label edgei : fEdges)
928 if (isBoundaryEdge[edgei])
930 ++nFaceBoundaryEdges;
934 if (nFaceBoundaryEdges == fEdges.size())
958 facePatch[facei] = nearestAdaptPatch[facei];
959 faceZone[facei] = nearestAdaptZone[facei];
960 if (isMasterFace.test(facei))
995 Info<<
"markFacesOnProblemCells : marked " 997 <<
" additional internal faces to be converted into baffles." 1002 Info<<
"markFacesOnProblemCells : prevented " 1004 <<
" internal faces from getting converted into baffles." 1012 void Foam::meshRefinement::markFacesOnProblemCellsGeometric
1014 const snapParameters& snapParams,
1015 const dictionary& motionDict,
1027 const labelList adaptPatchIDs(meshedPatches());
1030 autoPtr<indirectPrimitivePatch> ppPtr
1048 Info<<
"Constructing mesh displacer ..." <<
endl;
1049 Info<<
"Using mesh parameters " << motionDict <<
nl <<
endl;
1053 motionSmoother meshMover
1064 Info<<
"Checking initial mesh ..." <<
endl;
1080 Info<<
"Detected " << nInitErrors <<
" illegal faces" 1081 <<
" (concave, zero area or negative cell pyramid volume)" 1085 Info<<
"Checked initial mesh in = " 1086 << mesh_.time().cpuTimeIncrement() <<
" s\n" <<
nl <<
endl;
1102 snappySnapDriver::calcNearestSurface
1104 snapParams.strictRegionSnap(),
1106 globalToMasterPatch,
1115 const labelList& meshPoints = pp.meshPoints();
1120 newPoints[meshPoints[i]] += disp[i];
1127 minMagSqrEqOp<point>(),
1128 vector(GREAT, GREAT, GREAT)
1131 mesh_.movePoints(newPoints);
1134 mesh_.moving(
false);
1142 nearestPatch(meshedPatches(), nearestAdaptPatch, nearestAdaptZone);
1146 facePatch.setSize(mesh_.nFaces());
1148 faceZone.setSize(mesh_.nFaces());
1151 label nBaffleFaces = 0;
1154 faceSet wrongFaces(mesh_,
"wrongFaces", 100);
1162 label nWrongFaces = 0;
1194 scalar minArea(get<scalar>(motionDict,
"minArea", dryRun_));
1195 if (minArea > -SMALL)
1213 Info<<
" faces with area < " 1214 <<
setw(5) << minArea
1216 << nNewWrongFaces-nWrongFaces <<
endl;
1218 nWrongFaces = nNewWrongFaces;
1221 scalar minDet(get<scalar>(motionDict,
"minDeterminant", dryRun_));
1241 Info<<
" faces on cells with determinant < " 1242 <<
setw(5) << minDet <<
" : " 1243 << nNewWrongFaces-nWrongFaces <<
endl;
1245 nWrongFaces = nNewWrongFaces;
1250 for (
const label facei : wrongFaces)
1252 const label patchi = mesh_.boundaryMesh().whichPatch(facei);
1254 if (patchi == -1 || mesh_.boundaryMesh()[patchi].coupled())
1256 facePatch[facei] = nearestAdaptPatch[facei];
1257 faceZone[facei] = nearestAdaptZone[facei];
1270 mesh_.movePoints(oldPoints);
1273 mesh_.moving(
false);
1276 Info<<
"markFacesOnProblemCellsGeometric : marked " 1278 <<
" additional internal and coupled faces" 1279 <<
" to be converted into baffles." <<
endl;
static bool checkMesh(const bool report, const polyMesh &mesh, const dictionary &dict, labelHashSet &wrongFaces, const bool dryRun=false)
Check mesh with mesh settings in dict. Collects incorrect faces.
List< labelList > labelListList
A List of labelList.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
const labelListList & faceEdges() const
static const pointMesh & New(const polyMesh &mesh, Args &&... args)
Get existing or create a new MeshObject.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Unit conversion functions.
constexpr char nl
The newline '\n' character (0x0a)
label start() const noexcept
The start label of boundary faces in the polyMesh face list.
Ostream & endl(Ostream &os)
Add newline and flush stream.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field...
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
#define forAll(list, i)
Loop across all elements in list.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
void write(vtk::formatter &fmt, const Type &val, const label n=1)
Component-wise write of a value (N times)
vectorField pointField
pointField is a vectorField.
dimensionedScalar cos(const dimensionedScalar &ds)
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
static void preSmoothPatch(const meshRefinement &meshRefiner, const snapParameters &snapParams, const label nInitErrors, const List< labelPair > &baffles, motionSmoother &)
Smooth the mesh (patch and internal) to increase visibility.
virtual const faceList & faces() const
Return raw faces.
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
static scalarField calcSnapDistance(const fvMesh &mesh, const snapParameters &snapParams, const indirectPrimitivePatch &)
Calculate edge length per patch point.
dimensionedScalar sin(const dimensionedScalar &ds)
Istream and Ostream manipulators taking arguments.
static bool checkFaceArea(const bool report, const scalar minArea, const polyMesh &, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
Small faces.
int debug
Static debugging option.
static autoPtr< indirectPrimitivePatch > makePatch(const polyMesh &, const labelList &)
Create patch from set of patches.
vector point
Point is a vector.
const polyBoundaryMesh & patches
static labelList affectedCells(const polyMesh &, const labelList &changedFaces)
Helper function: get affected cells from faces.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< vector > vectorField
Specialisation of Field<T> for vector.
Omanip< int > setw(const int i)
List< label > labelList
A List of labels.
static tmp< pointVectorField > makeDisplacementField(const pointMesh &pMesh, const labelList &adaptPatchIDs)
Helper function to make a pointVectorField with correct.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
static bool checkCellDeterminant(const bool report, const scalar minDet, const polyMesh &, const vectorField &faceAreas, const labelList &checkFaces, const labelList &affectedCells, labelHashSet *setPtr)
Area of internal faces v.s. boundary faces.
List< bool > boolList
A List of bools.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
forAllConstIters(mixture.phases(), phase)
A HashTable to objects of type <T> with a label key.