52 const bitSet& internalOrCoupledFaces
61 solveScalar sumMagClosedBoundary = 0;
65 if (!internalOrCoupledFaces.
size() || !internalOrCoupledFaces[facei])
68 sumMagClosedBoundary +=
mag(areas[facei]);
81 Info<<
" ***Boundary openness " << openness
82 <<
" possible hole in boundary description."
91 Info<<
" Boundary openness " << openness <<
" OK."
114 label nErrorClosed = 0;
118 const cell& curCell = c[cI];
131 if (nErrorClosed > 0)
135 Info<<
" ***Cells with invalid face labels found, number of cells "
136 << nErrorClosed <<
endl;
156 scalar maxOpennessCell =
max(openness);
158 scalar maxAspectRatio =
max(aspectRatio);
163 if (openness[celli] > closedThreshold_)
173 if (aspectRatio[celli] > aspectThreshold_)
177 aspectSetPtr->
insert(celli);
195 Info<<
" ***Open cells found, max cell openness: "
196 << maxOpennessCell <<
", number of open cells " << nOpen
207 Info<<
" ***High aspect ratio cells found, Max aspect ratio: "
209 <<
", number of cells " << nAspect
218 Info<<
" Max cell openness = " << maxOpennessCell <<
" OK." <<
nl
219 <<
" Max aspect ratio = " << maxAspectRatio <<
" OK."
231 const bool detailedReport,
239 scalar minArea = GREAT;
240 scalar maxArea = -GREAT;
242 forAll(magFaceAreas, facei)
244 if (magFaceAreas[facei] < VSMALL)
252 if (isInternalFace(facei))
254 Pout<<
"Zero or negative face area detected for "
255 <<
"internal face "<< facei <<
" between cells "
256 << faceOwner()[facei] <<
" and "
257 << faceNeighbour()[facei]
258 <<
". Face area magnitude = " << magFaceAreas[facei]
263 Pout<<
"Zero or negative face area detected for "
264 <<
"boundary face " << facei <<
" next to cell "
265 << faceOwner()[facei] <<
". Face area magnitude = "
266 << magFaceAreas[facei] <<
endl;
271 minArea =
min(minArea, magFaceAreas[facei]);
272 maxArea =
max(maxArea, magFaceAreas[facei]);
278 if (minArea < VSMALL)
282 Info<<
" ***Zero or negative face area detected. "
283 "Minimum area: " << minArea <<
endl;
291 Info<<
" Minimum face area = " << minArea
292 <<
". Maximum face area = " << maxArea
293 <<
". Face area magnitudes OK." <<
endl;
304 const bool detailedReport,
310 scalar minVolume = GREAT;
311 scalar maxVolume = -GREAT;
313 label nNegVolCells = 0;
317 if (vols[celli] < VSMALL)
325 Pout<<
"Zero or negative cell volume detected for cell "
326 << celli <<
". Volume = " << vols[celli] <<
endl;
332 minVolume =
min(minVolume, vols[celli]);
333 maxVolume =
max(maxVolume, vols[celli]);
340 if (minVolume < VSMALL)
344 Info<<
" ***Zero or negative cell volume detected. "
345 <<
"Minimum negative volume: " << minVolume
346 <<
", Number of negative volume cells: " << nNegVolCells
355 Info<<
" Min volume = " << minVolume
356 <<
". Max volume = " << maxVolume
357 <<
". Total volume = " <<
gSum(vols)
358 <<
". Cell volumes OK." <<
endl;
384 const scalar severeNonorthogonalityThreshold =
387 scalar minDDotS =
min(ortho);
389 scalar sumDDotS =
sum(ortho);
391 label severeNonOrth = 0;
393 label errorNonOrth = 0;
398 if (ortho[facei] < severeNonorthogonalityThreshold)
400 if (ortho[facei] > SMALL)
428 label neiSize = ortho.
size();
435 Info<<
" Mesh non-orthogonality Max: "
442 if (severeNonOrth > 0)
444 Info<<
" *Number of severely non-orthogonal faces: "
445 << severeNonOrth <<
"." <<
endl;
449 if (errorNonOrth > 0)
453 Info<<
" ***Number of non-orthogonality errors: "
454 << errorNonOrth <<
"." <<
endl;
462 Info<<
" Non-orthogonality check OK." <<
endl;
474 const bool detailedReport,
475 const scalar minPyrVol,
498 label nErrorPyrs = 0;
502 if (ownPyrVol[facei] < minPyrVol)
510 Pout<<
"Negative pyramid volume: " << ownPyrVol[facei]
511 <<
" for face " << facei <<
" " <<
f[facei]
512 <<
" and owner cell: " << own[facei] <<
endl
513 <<
"Owner cell vertex labels: "
514 <<
cells()[own[facei]].labels(faces())
521 if (isInternalFace(facei))
523 if (neiPyrVol[facei] < minPyrVol)
531 Pout<<
"Negative pyramid volume: " << neiPyrVol[facei]
532 <<
" for face " << facei <<
" " <<
f[facei]
533 <<
" and neighbour cell: " << nei[facei] <<
nl
534 <<
"Neighbour cell vertex labels: "
535 <<
cells()[nei[facei]].labels(faces())
549 Info<<
" ***Error in face pyramids: "
550 << nErrorPyrs <<
" faces are incorrectly oriented."
559 Info<<
" Face pyramids OK." <<
endl;
591 scalar maxSkew =
max(skewness);
598 if (skewness[facei] > skewThreshold_)
616 Info<<
" ***Max skewness = " << maxSkew
617 <<
", " << nWarnSkew <<
" highly skew faces detected"
618 " which may impair the quality of the results"
627 Info<<
" Max skewness = " << maxSkew <<
" OK." <<
endl;
645 if (maxDeg < -SMALL || maxDeg > 180+SMALL)
648 <<
"maxDeg should be [0..180] but is now " << maxDeg
664 scalar maxEdgeSin =
max(faceAngles);
670 if (faceAngles[facei] > SMALL)
686 scalar maxConcaveDegr =
691 Info<<
" *There are " << nConcave
692 <<
" faces with concave angles between consecutive"
693 <<
" edges. Max concave angle = " << maxConcaveDegr
694 <<
" degrees." <<
endl;
702 Info<<
" All angles in faces OK." <<
endl;
715 const scalar warnFlatness,
721 if (warnFlatness < 0 || warnFlatness > 1)
724 <<
"warnFlatness should be [0..1] but is " << warnFlatness
741 scalar minFlatness = GREAT;
742 scalar sumFlatness = 0;
746 forAll(faceFlatness, facei)
748 if (fcs[facei].size() > 3 && magAreas[facei] > VSMALL)
750 sumFlatness += faceFlatness[facei];
753 minFlatness =
min(minFlatness, faceFlatness[facei]);
755 if (faceFlatness[facei] < warnFlatness)
778 Info<<
" Face flatness (1 = flat, 0 = butterfly) : min = "
779 << minFlatness <<
" average = " << sumFlatness / nSummed
789 Info<<
" *There are " << nWarped
790 <<
" faces with ratio between projected and actual area < "
791 << warnFlatness <<
endl;
793 Info<<
" Minimum ratio (minimum flatness, maximum warpage) = "
794 << minFlatness <<
endl;
802 Info<<
" All face flatness OK." <<
endl;
822 label nConcaveCells = 0;
826 const cell& cFaces = c[celli];
828 bool concave =
false;
837 label fI = cFaces[i];
839 const point& fC = fCentres[fI];
843 fN /=
max(
mag(fN), VSMALL);
847 if (fOwner[fI] != celli)
859 label fJ = cFaces[j];
861 const point& pt = fCentres[fJ];
871 pC /=
max(
mag(pC), VSMALL);
873 if ((pC & fN) > -planarCosAngle_)
895 if (nConcaveCells > 0)
899 Info<<
" ***Concave cells (using face planes) found,"
900 <<
" number of cells: " << nConcaveCells <<
endl;
908 Info<<
" Concave cell check OK." <<
endl;
929 label internal = nInternalFaces();
934 label nMultipleCells =
false;
938 for (label facei = 0; facei < internal; facei++)
940 if (own[facei] >= nei[facei])
964 label facei = curFaces[i];
966 if (facei >= nInternalFaces())
973 label nbrCelli = nei[facei];
975 if (nbrCelli == celli)
977 nbrCelli = own[facei];
980 if (celli < nbrCelli)
1000 label prevCell = nbr[0];
1001 label prevFace = curFaces[nbr.
indices()[0]];
1003 bool hasMultipleFaces =
false;
1005 for (label i = 1; i < nbr.
size(); i++)
1007 label thisCell = nbr[i];
1008 label thisFace = curFaces[nbr.
indices()[i]];
1015 if (thisCell == prevCell)
1017 hasMultipleFaces =
true;
1021 setPtr->
insert(prevFace);
1022 setPtr->
insert(thisFace);
1025 else if (thisFace < prevFace)
1031 setPtr->
insert(thisFace);
1035 prevCell = thisCell;
1036 prevFace = thisFace;
1039 if (hasMultipleFaces)
1048 if ((debug || report) && nMultipleCells > 0)
1050 Info<<
" <<Found " << nMultipleCells
1051 <<
" neighbouring cells with multiple inbetween faces." <<
endl;
1056 if (debug || report)
1058 Info<<
" ***Faces not in upper triangular order." <<
endl;
1064 if (debug || report)
1066 Info<<
" Upper triangular ordering OK." <<
endl;
1081 label nOpenCells = 0;
1090 const edgeList cellEdges = c[celli].edges(
f);
1096 edgeList curFaceEdges =
f[curFaces[facei]].edges();
1098 forAll(curFaceEdges, faceEdgeI)
1100 const edge& curEdge = curFaceEdges[faceEdgeI];
1102 forAll(cellEdges, cellEdgeI)
1104 if (cellEdges[cellEdgeI] == curEdge)
1106 edgeUsage[cellEdgeI]++;
1114 label nSingleEdges = 0;
1118 if (edgeUsage[edgeI] == 1)
1120 singleEdges[nSingleEdges] = cellEdges[edgeI];
1123 else if (edgeUsage[edgeI] != 2)
1132 if (nSingleEdges > 0)
1147 if (debug || report)
1149 Info<<
" ***Open cells found, number of cells: " << nOpenCells
1150 <<
". This problem may be fixable using the zipUpMesh utility."
1157 if (debug || report)
1159 Info<<
" Topological cell zip-up check OK." <<
endl;
1177 label nErrorFaces = 0;
1181 const face& curFace =
f[fI];
1198 bool inserted = facePoints.
insert(curFace[fp]);
1214 if (nErrorFaces > 0)
1216 if (debug || report)
1218 Info<<
" Faces with invalid vertex labels found, "
1219 <<
" number of faces: " << nErrorFaces <<
endl;
1225 if (debug || report)
1227 Info<<
" Face vertices OK." <<
endl;
1242 label nFaceErrors = 0;
1243 label nCellErrors = 0;
1249 if (pf[pointi].empty())
1279 if (nFaceErrors > 0 || nCellErrors > 0)
1281 if (debug || report)
1283 Info<<
" ***Unused points found in the mesh, "
1284 "number unused by faces: " << nFaceErrors
1285 <<
" number unused by cells: " << nCellErrors
1292 if (debug || report)
1305 label& nBaffleFaces,
1313 label nbFacei = iter.key();
1314 label nCommon = iter();
1316 const face& curFace = faces()[facei];
1317 const face& nbFace = faces()[nbFacei];
1319 if (nCommon == nbFace.
size() || nCommon == curFace.
size())
1321 if (nbFace.
size() != curFace.
size())
1353 label nbFacei = iter.key();
1354 label nCommon = iter();
1356 const face& curFace = faces()[facei];
1357 const face& nbFace = faces()[nbFacei];
1362 && nCommon != nbFace.
size()
1363 && nCommon != curFace.
size()
1369 label nb = nbFace.
find(curFace[fp]);
1383 label fpPlus1 = curFace.
fcIndex(fp);
1384 label fpMin1 = curFace.
rcIndex(fp);
1387 label nbPlus1 = nbFace.
fcIndex(nb);
1388 label nbMin1 = nbFace.
rcIndex(nb);
1395 if (nbFace[nbPlus1] == curFace[fpPlus1])
1400 else if (nbFace[nbPlus1] == curFace[fpMin1])
1405 else if (nbFace[nbMin1] == curFace[fpMin1])
1425 if (curFp >= curFace.
size())
1431 curFp = curFace.
size()-1;
1436 if (curNb >= nbFace.
size())
1442 curNb = nbFace.
size()-1;
1444 }
while (curFace[curFp] == nbFace[curNb]);
1453 for (label commonI = 0; commonI < nCommon; commonI++)
1457 if (curFp >= curFace.
size())
1463 curFp = curFace.
size()-1;
1468 if (curNb >= nbFace.
size())
1474 curNb = nbFace.
size()-1;
1477 if (curFace[curFp] != nbFace[curNb])
1513 label nBaffleFaces = 0;
1514 label nErrorDuplicate = 0;
1515 label nErrorOrder = 0;
1518 for (label facei = 0; facei <
nFaces(); facei++)
1520 const face& curFace = faces()[facei];
1524 nCommonPoints.
clear();
1528 label pointi = curFace[fp];
1534 label nbFacei = nbs[nbI];
1536 if (facei < nbFacei)
1539 ++(nCommonPoints(nbFacei, 0));
1547 if (checkDuplicateFaces(facei, nCommonPoints, nBaffleFaces, setPtr))
1553 if (checkCommonOrder(facei, nCommonPoints, setPtr))
1565 Info<<
" Number of identical duplicate faces (baffle faces): "
1566 << nBaffleFaces <<
endl;
1569 if (nErrorDuplicate > 0 || nErrorOrder > 0)
1572 if (nErrorDuplicate > 0)
1574 Info<<
" <<Number of duplicate (not baffle) faces found: "
1576 <<
". This might indicate a problem." <<
endl;
1579 if (nErrorOrder > 0)
1581 Info<<
" <<Number of faces with non-consecutive shared points: "
1582 << nErrorOrder <<
". This might indicate a problem." <<
endl;
1588 if (debug || report)
1590 Info<<
" Face-face connectivity OK." <<
endl;
1601 return checkClosedBoundary(faceAreas(), report,
bitSet());
1613 return checkClosedCells
1631 return checkFaceAreas
1647 return checkCellVolumes
1663 return checkFaceOrthogonality
1676 const scalar minPyrVol,
1680 return checkFacePyramids
1698 return checkFaceSkewness
1713 const scalar maxDeg,
1717 return checkFaceAngles
1731 const scalar warnFlatness,
1735 return checkFaceFlatness
1753 return checkConcaveCells
1765 label nFailedChecks = 0;
1767 if (checkPoints(report)) ++nFailedChecks;
1768 if (checkUpperTriangular(report)) ++nFailedChecks;
1769 if (checkCellsZipUp(report)) ++nFailedChecks;
1770 if (checkFaceVertices(report)) ++nFailedChecks;
1771 if (checkFaceFaces(report)) ++nFailedChecks;
1775 if (debug || report)
1777 Info<<
" Failed " << nFailedChecks
1778 <<
" mesh topology checks." <<
endl;
1784 if (debug || report)
1786 Info<<
" Mesh topology OK." <<
endl;
1795 label nFailedChecks = 0;
1797 if (checkClosedBoundary(report)) ++nFailedChecks;
1798 if (checkClosedCells(report)) ++nFailedChecks;
1799 if (checkFaceAreas(report)) ++nFailedChecks;
1800 if (checkCellVolumes(report)) ++nFailedChecks;
1801 if (checkFaceOrthogonality(report)) ++nFailedChecks;
1802 if (checkFacePyramids(report)) ++nFailedChecks;
1803 if (checkFaceSkewness(report)) ++nFailedChecks;
1807 if (debug || report)
1809 Info<<
" Failed " << nFailedChecks
1810 <<
" mesh geometry checks." <<
endl;
1816 if (debug || report)
1818 Info<<
" Mesh geometry OK." <<
endl;
1833 if (debug || report)
1835 Info<<
" Failed " << nFailedChecks
1836 <<
" mesh checks." <<
endl;
1842 if (debug || report)
1853 scalar prev = closedThreshold_;
1854 closedThreshold_ = val;
1862 scalar prev = aspectThreshold_;
1863 aspectThreshold_ = val;
1871 scalar prev = nonOrthThreshold_;
1872 nonOrthThreshold_ = val;
1880 scalar prev = skewThreshold_;
1881 skewThreshold_ = val;
Various functions to operate on Lists.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
void clear()
Clear all entries from table.
A HashTable to objects of type <T> with a label key.
label size() const noexcept
Number of entries.
A list that is sorted upon construction or when explicitly requested with the sort() method.
const labelList & indices() const noexcept
Return the list of sorted indices. Updated every sort.
void sort()
Forward (stable) sort the list (if changed after construction).
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
label rcIndex(const label i) const noexcept
void size(const label n)
Older name for setAddressableSize.
label fcIndex(const label i) const noexcept
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
tmp< scalarField > faceSkewness() const
Return face skewness.
A cell is defined as a list of faces with extra functionality.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Class to handle errors and exceptions in a simple, consistent stream-based manner.
A face is a list of labels corresponding to mesh vertices.
void checkMesh() const
Debug: Check coupled mesh for correctness.
Smooth ATC in cells having a point to a set of patches supplied by type.
bool checkFaceOrthogonality(const vectorField &fAreas, const vectorField &cellCtrs, const bool report, labelHashSet *setPtr) const
Check for non-orthogonality.
static scalar planarCosAngle_
Threshold where faces are considered coplanar.
virtual bool checkCellsZipUp(const bool report=false, labelHashSet *setPtr=nullptr) const
Check cell zip-up.
static scalar skewThreshold_
Skewness warning threshold.
static scalar setNonOrthThreshold(const scalar)
Set the non-orthogonality warning threshold in degrees.
label nInternalFaces() const noexcept
Number of internal faces.
virtual bool checkTopology(const bool report=false) const
Check mesh topology for correctness.
bool checkCommonOrder(const label, const Map< label > &, labelHashSet *) const
Check that shared points are in consecutive order.
virtual bool checkUpperTriangular(const bool report=false, labelHashSet *setPtr=nullptr) const
Check face ordering.
static scalar closedThreshold_
Static data to control mesh checking.
bool checkFacePyramids(const pointField &points, const vectorField &ctrs, const bool report, const bool detailedReport, const scalar minPyrVol, labelHashSet *setPtr) const
Check face pyramid volume.
bool checkFaceAngles(const pointField &points, const vectorField &faceAreas, const bool report, const scalar maxDeg, labelHashSet *setPtr) const
Check face angles.
static scalar nonOrthThreshold_
Non-orthogonality warning threshold in deg.
static scalar setSkewThreshold(const scalar)
Set the skewness warning threshold as percentage.
static scalar setAspectThreshold(const scalar)
Set the aspect ratio warning threshold.
bool checkDuplicateFaces(const label, const Map< label > &, label &nBaffleFaces, labelHashSet *) const
Check if all points on face are shared with another face.
virtual bool checkFaceVertices(const bool report=false, labelHashSet *setPtr=nullptr) const
Check uniqueness of face vertices.
bool checkCellVolumes(const scalarField &vols, const bool report, const bool detailedReport, labelHashSet *setPtr) const
Check for negative cell volumes.
bool checkFaceAreas(const vectorField &faceAreas, const bool report, const bool detailedReport, labelHashSet *setPtr) const
Check for negative face areas.
virtual bool checkPoints(const bool report=false, labelHashSet *setPtr=nullptr) const
Check for unused points.
bool checkFaceSkewness(const pointField &points, const vectorField &fCtrs, const vectorField &fAreas, const vectorField &cellCtrs, const bool report, labelHashSet *setPtr) const
Check face skewness.
bool checkFaceFlatness(const pointField &points, const vectorField &faceCentres, const vectorField &faceAreas, const bool report, const scalar warnFlatness, labelHashSet *setPtr) const
Check face warpage.
static scalar aspectThreshold_
Aspect ratio warning threshold.
bool checkConcaveCells(const vectorField &fAreas, const pointField &fCentres, const bool report, labelHashSet *setPtr) const
Check for concave cells by the planes of faces.
bool checkClosedBoundary(const vectorField &areas, const bool report, const bitSet &internalOrCoupledFaces) const
Check boundary for closedness.
virtual bool checkFaceFaces(const bool report=false, labelHashSet *setPtr=nullptr) const
Check face-face connectivity.
bool checkClosedCells(const vectorField &faceAreas, const scalarField &cellVolumes, const bool report, labelHashSet *setPtr, labelHashSet *aspectSetPtr, const Vector< label > &meshD) const
Check cells for closedness.
virtual bool checkGeometry(const bool report=false) const
Check mesh geometry (& implicitly topology) for correctness.
static scalar setClosedThreshold(const scalar)
Set the closedness ratio warning threshold.
A class for managing temporary objects.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
#define DebugInFunction
Report an information message using Foam::Info.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedScalar asin(const dimensionedScalar &ds)
Type gSum(const FieldField< Field, Type > &f)
label checkGeometry(const polyMesh &mesh, const bool allGeometry, autoPtr< surfaceWriter > &surfWriter, autoPtr< coordSetWriter > &setWriter)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
dimensionedScalar sin(const dimensionedScalar &ds)
label checkTopology(const polyMesh &mesh, const bool allTopology, const bool allGeometry, autoPtr< surfaceWriter > &surfWriter, autoPtr< coordSetWriter > &setWriter)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
constexpr scalar radToDeg() noexcept
Multiplication factor for radians to degrees conversion.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
static constexpr const zero Zero
Global zero (0)
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar acos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Unit conversion functions.