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])
184 if (
mag(n0 & n1) < 0.1)
186 candidateFaces.append(face0);
189 else if (cellLevel[cell1] > cellLevel[cell0])
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 for (
const label facei : mesh_.cells()[iter.key()])
482 facePatch[facei] == -1
487 facePatch[facei] = nearestAdaptPatch[facei];
488 faceZone[facei] = nearestAdaptZone[facei];
502 Info<<
"markFacesOnProblemCells : Marked " 504 <<
" additional internal faces to be converted into baffles" 507 <<
" cells edge-connected to lower level cells." <<
endl;
511 cellSet problemCellSet(mesh_,
"problemCells", problemCells.toc());
512 problemCellSet.instance() =
timeName();
513 Pout<<
"Writing " << problemCellSet.size()
514 <<
" cells that are edge connected to coarser cell to set " 515 << problemCellSet.objectPath() <<
endl;
516 problemCellSet.
write();
548 const scalar volFraction =
549 motionDict.getOrDefault<scalar>(
"minVolCollapseRatio", -1);
551 const bool checkCollapse = (volFraction > 0);
553 scalar maxNonOrtho = -1;
561 minArea = get<scalar>(motionDict,
"minArea", dryRun_);
562 maxNonOrtho = get<scalar>(motionDict,
"maxNonOrtho", dryRun_);
564 Info<<
"markFacesOnProblemCells :" 565 <<
" Deleting all-anchor surface cells only if" 566 <<
" snapping them violates mesh quality constraints:" <<
nl 567 <<
" snapped/original cell volume < " << volFraction <<
nl 568 <<
" face area < " << minArea <<
nl 569 <<
" non-orthogonality > " << maxNonOrtho <<
nl 573 autoPtr<indirectPrimitivePatch> ppPtr
585 List<pointIndexHit> hitInfo;
587 surfaces_.findNearest
597 newPoints = mesh_.points();
601 if (hitInfo[i].hit())
603 newPoints[meshPoints[i]] = hitInfo[i].point();
609 const_cast<Time&
>(mesh_.time())++;
611 mesh_.movePoints(newPoints);
621 writeType(writeLevel() | WRITEMESH),
622 mesh_.time().path()/
"newPoints" 624 mesh_.movePoints(oldPoints);
645 bitSet hasSevenBoundaryAnchorPoints(mesh_.nCells());
650 DynamicList<label> dynFEdges;
651 DynamicList<label> dynCPoints;
656 const labelList& cPoints = mesh_.cellPoints(celli, pSet, dynCPoints);
660 label nBoundaryAnchors = 0;
662 label nonBoundaryAnchor = -1;
664 for (
const label pointi : cPoints)
666 if (pointLevel[pointi] <= cellLevel[celli])
669 if (isBoundaryPoint[pointi])
676 nonBoundaryAnchor = pointi;
685 if (nBoundaryAnchors == 8)
687 const cell& cFaces = mesh_.cells()[celli];
710 && !isCollapsedCell(newPoints, volFraction, celli)
725 for (
const label facei : cFaces)
731 facePatch[facei] == -1
736 facePatch[facei] = nearestAdaptPatch[facei];
737 faceZone[facei] = nearestAdaptZone[facei];
753 else if (nBoundaryAnchors == 7 && nonBoundaryAnchor != -1)
756 hasSevenBoundaryAnchorPoints.set(celli);
757 nonBoundaryAnchors.insert(nonBoundaryAnchor);
766 DynamicList<label> dynPCells;
768 for (
const label pointi : nonBoundaryAnchors)
770 const labelList& pCells = mesh_.pointCells(pointi, dynPCells);
775 for (
const label celli : pCells)
777 if (hasSevenBoundaryAnchorPoints.test(celli))
786 for (
const label celli : pCells)
788 if (hasSevenBoundaryAnchorPoints.test(celli))
793 && !isCollapsedCell(newPoints, volFraction, celli)
808 const cell& cFaces = mesh_.cells()[celli];
810 for (
const label facei : cFaces)
816 facePatch[facei] == -1
821 facePatch[facei] = nearestAdaptPatch[facei];
822 faceZone[facei] = nearestAdaptZone[facei];
870 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
872 if (facePatch[facei] == -1)
874 const labelList& fEdges = mesh_.faceEdges(facei, dynFEdges);
875 label nFaceBoundaryEdges = 0;
877 for (
const label edgei : fEdges)
879 if (isBoundaryEdge[edgei])
881 ++nFaceBoundaryEdges;
885 if (nFaceBoundaryEdges == fEdges.size())
909 facePatch[facei] = nearestAdaptPatch[facei];
910 faceZone[facei] = nearestAdaptZone[facei];
924 label facei =
pp.start();
928 if (facePatch[facei] == -1)
930 const labelList& fEdges = mesh_.faceEdges(facei, dynFEdges);
931 label nFaceBoundaryEdges = 0;
933 for (
const label edgei : fEdges)
935 if (isBoundaryEdge[edgei])
937 ++nFaceBoundaryEdges;
941 if (nFaceBoundaryEdges == fEdges.size())
965 facePatch[facei] = nearestAdaptPatch[facei];
966 faceZone[facei] = nearestAdaptZone[facei];
967 if (isMasterFace.test(facei))
1002 Info<<
"markFacesOnProblemCells : marked " 1004 <<
" additional internal faces to be converted into baffles." 1009 Info<<
"markFacesOnProblemCells : prevented " 1011 <<
" internal faces from getting converted into baffles." 1019 void Foam::meshRefinement::markFacesOnProblemCellsGeometric
1021 const snapParameters& snapParams,
1022 const dictionary& motionDict,
1034 const labelList adaptPatchIDs(meshedPatches());
1037 autoPtr<indirectPrimitivePatch> ppPtr
1055 Info<<
"Constructing mesh displacer ..." <<
endl;
1056 Info<<
"Using mesh parameters " << motionDict <<
nl <<
endl;
1060 motionSmoother meshMover
1071 Info<<
"Checking initial mesh ..." <<
endl;
1087 Info<<
"Detected " << nInitErrors <<
" illegal faces" 1088 <<
" (concave, zero area or negative cell pyramid volume)" 1092 Info<<
"Checked initial mesh in = " 1093 << mesh_.time().cpuTimeIncrement() <<
" s\n" <<
nl <<
endl;
1109 snappySnapDriver::calcNearestSurface
1111 snapParams.strictRegionSnap(),
1113 globalToMasterPatch,
1127 newPoints[meshPoints[i]] += disp[i];
1134 minMagSqrEqOp<point>(),
1135 vector(GREAT, GREAT, GREAT)
1138 mesh_.movePoints(newPoints);
1141 mesh_.moving(
false);
1149 nearestPatch(meshedPatches(), nearestAdaptPatch, nearestAdaptZone);
1153 facePatch.setSize(mesh_.nFaces());
1155 faceZone.setSize(mesh_.nFaces());
1158 label nBaffleFaces = 0;
1161 faceSet wrongFaces(mesh_,
"wrongFaces", 100);
1169 label nWrongFaces = 0;
1201 scalar minArea(get<scalar>(motionDict,
"minArea", dryRun_));
1202 if (minArea > -SMALL)
1220 Info<<
" faces with area < " 1221 <<
setw(5) << minArea
1223 << nNewWrongFaces-nWrongFaces <<
endl;
1225 nWrongFaces = nNewWrongFaces;
1228 scalar minDet(get<scalar>(motionDict,
"minDeterminant", dryRun_));
1248 Info<<
" faces on cells with determinant < " 1249 <<
setw(5) << minDet <<
" : " 1250 << nNewWrongFaces-nWrongFaces <<
endl;
1252 nWrongFaces = nNewWrongFaces;
1257 for (
const label facei : wrongFaces)
1259 const label patchi = mesh_.boundaryMesh().whichPatch(facei);
1261 if (patchi == -1 || mesh_.boundaryMesh()[patchi].coupled())
1263 facePatch[facei] = nearestAdaptPatch[facei];
1264 faceZone[facei] = nearestAdaptZone[facei];
1277 mesh_.movePoints(oldPoints);
1280 mesh_.moving(
false);
1283 Info<<
"markFacesOnProblemCellsGeometric : marked " 1285 <<
" additional internal and coupled faces" 1286 <<
" 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.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
const Field< point_type > & localPoints() const
Return pointField of points in patch.
const labelListList & faceEdges() const
static const pointMesh & New(const polyMesh &mesh, Args &&... args)
Get existing or create MeshObject registered with typeName.
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)
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.
List< labelList > labelListList
List of labelList.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
#define forAll(list, i)
Loop across all elements in list.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
const Field< point_type > & faceNormals() const
Return face unit normals for patch.
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), works like std::iota() but returning a...
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
const labelListList & edgeFaces() const
Return edge-face addressing.
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.
label whichPatch(const label meshFacei) const
Return patch index for a given mesh face index. Uses binary search.
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.
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
forAllConstIters(mixture.phases(), phase)
A HashTable to objects of type <T> with a label key.