Go to the documentation of this file.
52 const bitSet& internalOrCoupledFaces
61 solveScalar sumMagClosedBoundary = 0;
63 for (label facei = nInternalFaces(); facei < areas.size(); facei++)
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];
120 if (
min(curCell) < 0 ||
max(curCell) > nFaces())
131 if (nErrorClosed > 0)
135 Info<<
" ***Cells with invalid face labels found, number of cells "
136 << nErrorClosed <<
endl;
145 primitiveMeshTools::cellClosedness
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,
488 primitiveMeshTools::facePyramidVolume
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]++;
1113 edgeList singleEdges(cellEdges.size());
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;
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;
int debug
Static debugging option.
bool checkFaceFlatness(const pointField &points, const vectorField &faceCentres, const vectorField &faceAreas, const bool report, const scalar warnFlatness, labelHashSet *setPtr) const
Check face warpage.
bool checkFaceSkewness(const pointField &points, const vectorField &fCtrs, const vectorField &fAreas, const vectorField &cellCtrs, const bool report, labelHashSet *setPtr) const
Check face skewness.
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 checkFaceAreas(const vectorField &faceAreas, const bool report, const bool detailedReport, labelHashSet *setPtr) const
Check for negative face areas.
virtual bool checkFaceFaces(const bool report=false, labelHashSet *setPtr=nullptr) const
Check face-face connectivity.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
virtual bool checkMesh(const bool report=false) const
Check mesh for correctness. Returns false for no error.
A class for managing temporary objects.
static constexpr const zero Zero
Global zero (0)
dimensionedScalar sin(const dimensionedScalar &ds)
static scalar setSkewThreshold(const scalar)
Set the skewness warning threshold as percentage.
bool checkDuplicateFaces(const label, const Map< label > &, label &nBaffleFaces, labelHashSet *) const
Check if all points on face are shared with another face.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
bool checkFaceAngles(const pointField &points, const vectorField &faceAreas, const bool report, const scalar maxDeg, labelHashSet *setPtr) const
Check face angles.
Unit conversion functions.
bool checkCellVolumes(const scalarField &vols, const bool report, const bool detailedReport, labelHashSet *setPtr) const
Check for negative cell volumes.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Type gSum(const FieldField< Field, Type > &f)
label checkTopology(const polyMesh &mesh, const bool allTopology, const bool allGeometry, autoPtr< surfaceWriter > &surfWriter, const autoPtr< writer< scalar >> &setWriter)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
virtual bool checkTopology(const bool report=false) const
Check mesh topology for correctness.
#define forAll(list, i)
Loop across all elements in list.
bool checkConcaveCells(const vectorField &fAreas, const pointField &fCentres, const bool report, labelHashSet *setPtr) const
Check for concave cells by the planes of faces.
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
static scalar setClosedThreshold(const scalar)
Set the closedness ratio warning threshold.
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
virtual bool checkGeometry(const bool report=false) const
Check mesh geometry (& implicitly topology) for correctness.
static scalar setNonOrthThreshold(const scalar)
Set the non-orthogonality warning threshold in degrees.
messageStream Info
Information stream (stdout output on master, null elsewhere)
#define DebugInFunction
Report an information message using Foam::Info.
constexpr scalar radToDeg(const scalar rad) noexcept
Conversion from radians to degrees.
label checkGeometry(const polyMesh &mesh, const bool allGeometry, autoPtr< surfaceWriter > &surfWriter, const autoPtr< writer< scalar >> &setWriter)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
reduce(hasMovingMesh, orOp< bool >())
A list that is sorted upon construction or when explicitly requested with the sort() method.
static scalar aspectThreshold_
Aspect ratio warning threshold.
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
errorManipArg< error, int > exit(error &err, const int errNo=1)
virtual bool checkPoints(const bool report=false, labelHashSet *setPtr=nullptr) const
Check for unused points.
static scalar closedThreshold_
Static data to control mesh checking.
static scalar nonOrthThreshold_
Non-orthogonality warning threshold in deg.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
forAllConstIters(mixture.phases(), phase)
static scalar planarCosAngle_
Threshold where faces are considered coplanar.
static scalar setAspectThreshold(const scalar)
Set the aspect ratio warning threshold.
static scalar skewThreshold_
Skewness warning threshold.
bool checkClosedCells(const vectorField &faceAreas, const scalarField &cellVolumes, const bool report, labelHashSet *setPtr, labelHashSet *aspectSetPtr, const Vector< label > &meshD) const
Check cells for closedness.
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
bool checkClosedBoundary(const vectorField &areas, const bool report, const bitSet &internalOrCoupledFaces) const
Check boundary for closedness.
dimensionedScalar acos(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label size() const noexcept
Number of entries.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
const dimensionedScalar c
Speed of light in a vacuum.
A face is a list of labels corresponding to mesh vertices.
bool checkFaceOrthogonality(const vectorField &fAreas, const vectorField &cellCtrs, const bool report, labelHashSet *setPtr) const
Check for non-orthogonality.
Smooth ATC in cells having a point to a set of patches supplied by type.
Various functions to operate on Lists.
virtual bool checkFaceVertices(const bool report=false, labelHashSet *setPtr=nullptr) const
Check uniqueness of face vertices.
A cell is defined as a list of faces with extra functionality.
Class to handle errors and exceptions in a simple, consistent stream-based manner.
dimensionedScalar asin(const dimensionedScalar &ds)
virtual bool checkCellsZipUp(const bool report=false, labelHashSet *setPtr=nullptr) const
Check cell zip-up.
bool checkCommonOrder(const label, const Map< label > &, labelHashSet *) const
Check that shared points are in consecutive order.
dimensionedScalar cos(const dimensionedScalar &ds)
const vectorField & faceAreas() const
virtual bool checkUpperTriangular(const bool report=false, labelHashSet *setPtr=nullptr) const
Check face ordering.