Go to the documentation of this file.
52 const bitSet& internalOrCoupledFaces
58 <<
"Checking whether the boundary is closed" <<
endl;
65 solveScalar sumMagClosedBoundary = 0;
67 for (
label facei = nInternalFaces(); facei < areas.size(); facei++)
69 if (!internalOrCoupledFaces.
size() || !internalOrCoupledFaces[facei])
72 sumMagClosedBoundary +=
mag(areas[facei]);
85 Info<<
" ***Boundary openness " << openness
86 <<
" possible hole in boundary description."
95 Info<<
" Boundary openness " << openness <<
" OK."
116 <<
"Checking whether cells are closed" <<
endl;
122 label nErrorClosed = 0;
126 const cell& curCell =
c[cI];
128 if (
min(curCell) < 0 ||
max(curCell) > nFaces())
139 if (nErrorClosed > 0)
143 Info<<
" ***Cells with invalid face labels found, number of cells "
144 << nErrorClosed <<
endl;
153 primitiveMeshTools::cellClosedness
164 scalar maxOpennessCell =
max(openness);
166 scalar maxAspectRatio =
max(aspectRatio);
171 if (openness[celli] > closedThreshold_)
181 if (aspectRatio[celli] > aspectThreshold_)
185 aspectSetPtr->
insert(celli);
203 Info<<
" ***Open cells found, max cell openness: "
204 << maxOpennessCell <<
", number of open cells " << nOpen
215 Info<<
" ***High aspect ratio cells found, Max aspect ratio: "
217 <<
", number of cells " << nAspect
226 Info<<
" Max cell openness = " << maxOpennessCell <<
" OK." <<
nl
227 <<
" Max aspect ratio = " << maxAspectRatio <<
" OK."
239 const bool detailedReport,
250 scalar minArea = GREAT;
251 scalar maxArea = -GREAT;
253 forAll(magFaceAreas, facei)
255 if (magFaceAreas[facei] < VSMALL)
263 if (isInternalFace(facei))
265 Pout<<
"Zero or negative face area detected for "
266 <<
"internal face "<< facei <<
" between cells "
267 << faceOwner()[facei] <<
" and "
268 << faceNeighbour()[facei]
269 <<
". Face area magnitude = " << magFaceAreas[facei]
274 Pout<<
"Zero or negative face area detected for "
275 <<
"boundary face " << facei <<
" next to cell "
276 << faceOwner()[facei] <<
". Face area magnitude = "
277 << magFaceAreas[facei] <<
endl;
282 minArea =
min(minArea, magFaceAreas[facei]);
283 maxArea =
max(maxArea, magFaceAreas[facei]);
289 if (minArea < VSMALL)
293 Info<<
" ***Zero or negative face area detected. "
294 "Minimum area: " << minArea <<
endl;
302 Info<<
" Minimum face area = " << minArea
303 <<
". Maximum face area = " << maxArea
304 <<
". Face area magnitudes OK." <<
endl;
315 const bool detailedReport,
324 scalar minVolume = GREAT;
325 scalar maxVolume = -GREAT;
327 label nNegVolCells = 0;
331 if (vols[celli] < VSMALL)
339 Pout<<
"Zero or negative cell volume detected for cell "
340 << celli <<
". Volume = " << vols[celli] <<
endl;
346 minVolume =
min(minVolume, vols[celli]);
347 maxVolume =
max(maxVolume, vols[celli]);
354 if (minVolume < VSMALL)
358 Info<<
" ***Zero or negative cell volume detected. "
359 <<
"Minimum negative volume: " << minVolume
360 <<
", Number of negative volume cells: " << nNegVolCells
369 Info<<
" Min volume = " << minVolume
370 <<
". Max volume = " << maxVolume
371 <<
". Total volume = " <<
gSum(vols)
372 <<
". Cell volumes OK." <<
endl;
402 const scalar severeNonorthogonalityThreshold =
405 scalar minDDotS =
min(ortho);
407 scalar sumDDotS =
sum(ortho);
409 label severeNonOrth = 0;
411 label errorNonOrth = 0;
416 if (ortho[facei] < severeNonorthogonalityThreshold)
418 if (ortho[facei] > SMALL)
446 label neiSize = ortho.size();
453 Info<<
" Mesh non-orthogonality Max: "
460 if (severeNonOrth > 0)
462 Info<<
" *Number of severely non-orthogonal faces: "
463 << severeNonOrth <<
"." <<
endl;
467 if (errorNonOrth > 0)
471 Info<<
" ***Number of non-orthogonality errors: "
472 << errorNonOrth <<
"." <<
endl;
480 Info<<
" Non-orthogonality check OK." <<
endl;
492 const bool detailedReport,
493 const scalar minPyrVol,
509 primitiveMeshTools::facePyramidVolume
519 label nErrorPyrs = 0;
523 if (ownPyrVol[facei] < minPyrVol)
531 Pout<<
"Negative pyramid volume: " << ownPyrVol[facei]
532 <<
" for face " << facei <<
" " <<
f[facei]
533 <<
" and owner cell: " << own[facei] <<
endl
534 <<
"Owner cell vertex labels: "
535 <<
cells()[own[facei]].labels(faces())
542 if (isInternalFace(facei))
544 if (neiPyrVol[facei] < minPyrVol)
552 Pout<<
"Negative pyramid volume: " << neiPyrVol[facei]
553 <<
" for face " << facei <<
" " <<
f[facei]
554 <<
" and neighbour cell: " << nei[facei] <<
nl
555 <<
"Neighbour cell vertex labels: "
556 <<
cells()[nei[facei]].labels(faces())
570 Info<<
" ***Error in face pyramids: "
571 << nErrorPyrs <<
" faces are incorrectly oriented."
580 Info<<
" Face pyramids OK." <<
endl;
615 scalar maxSkew =
max(skewness);
622 if (skewness[facei] > skewThreshold_)
640 Info<<
" ***Max skewness = " << maxSkew
641 <<
", " << nWarnSkew <<
" highly skew faces detected"
642 " which may impair the quality of the results"
651 Info<<
" Max skewness = " << maxSkew <<
" OK." <<
endl;
672 if (maxDeg < -SMALL || maxDeg > 180+SMALL)
675 <<
"maxDeg should be [0..180] but is now " << maxDeg
691 scalar maxEdgeSin =
max(faceAngles);
697 if (faceAngles[facei] > SMALL)
713 scalar maxConcaveDegr =
718 Info<<
" *There are " << nConcave
719 <<
" faces with concave angles between consecutive"
720 <<
" edges. Max concave angle = " << maxConcaveDegr
721 <<
" degrees." <<
endl;
729 Info<<
" All angles in faces OK." <<
endl;
742 const scalar warnFlatness,
751 if (warnFlatness < 0 || warnFlatness > 1)
754 <<
"warnFlatness should be [0..1] but is now " << warnFlatness
771 scalar minFlatness = GREAT;
772 scalar sumFlatness = 0;
776 forAll(faceFlatness, facei)
778 if (fcs[facei].size() > 3 && magAreas[facei] > VSMALL)
780 sumFlatness += faceFlatness[facei];
783 minFlatness =
min(minFlatness, faceFlatness[facei]);
785 if (faceFlatness[facei] < warnFlatness)
808 Info<<
" Face flatness (1 = flat, 0 = butterfly) : min = "
809 << minFlatness <<
" average = " << sumFlatness / nSummed
819 Info<<
" *There are " << nWarped
820 <<
" faces with ratio between projected and actual area < "
821 << warnFlatness <<
endl;
823 Info<<
" Minimum ratio (minimum flatness, maximum warpage) = "
824 << minFlatness <<
endl;
832 Info<<
" All face flatness OK." <<
endl;
855 label nConcaveCells = 0;
859 const cell& cFaces =
c[celli];
861 bool concave =
false;
870 label fI = cFaces[i];
872 const point& fC = fCentres[fI];
876 fN /=
max(
mag(fN), VSMALL);
880 if (fOwner[fI] != celli)
892 label fJ = cFaces[j];
894 const point& pt = fCentres[fJ];
904 pC /=
max(
mag(pC), VSMALL);
906 if ((pC & fN) > -planarCosAngle_)
928 if (nConcaveCells > 0)
932 Info<<
" ***Concave cells (using face planes) found,"
933 <<
" number of cells: " << nConcaveCells <<
endl;
941 Info<<
" Concave cell check OK." <<
endl;
965 label internal = nInternalFaces();
970 label nMultipleCells =
false;
974 for (
label facei = 0; facei <
internal; facei++)
976 if (own[facei] >= nei[facei])
1000 label facei = curFaces[i];
1002 if (facei >= nInternalFaces())
1009 label nbrCelli = nei[facei];
1011 if (nbrCelli == celli)
1013 nbrCelli = own[facei];
1016 if (celli < nbrCelli)
1036 label prevCell = nbr[0];
1037 label prevFace = curFaces[nbr.indices()[0]];
1039 bool hasMultipleFaces =
false;
1041 for (
label i = 1; i < nbr.size(); i++)
1043 label thisCell = nbr[i];
1044 label thisFace = curFaces[nbr.indices()[i]];
1051 if (thisCell == prevCell)
1053 hasMultipleFaces =
true;
1057 setPtr->
insert(prevFace);
1058 setPtr->
insert(thisFace);
1061 else if (thisFace < prevFace)
1067 setPtr->
insert(thisFace);
1071 prevCell = thisCell;
1072 prevFace = thisFace;
1075 if (hasMultipleFaces)
1084 if ((
debug || report) && nMultipleCells > 0)
1086 Info<<
" <<Found " << nMultipleCells
1087 <<
" neighbouring cells with multiple inbetween faces." <<
endl;
1092 if (
debug || report)
1094 Info<<
" ***Faces not in upper triangular order." <<
endl;
1100 if (
debug || report)
1102 Info<<
" Upper triangular ordering OK." <<
endl;
1120 label nOpenCells = 0;
1129 const edgeList cellEdges =
c[celli].edges(
f);
1135 edgeList curFaceEdges =
f[curFaces[facei]].edges();
1137 forAll(curFaceEdges, faceEdgeI)
1139 const edge& curEdge = curFaceEdges[faceEdgeI];
1141 forAll(cellEdges, cellEdgeI)
1143 if (cellEdges[cellEdgeI] == curEdge)
1145 edgeUsage[cellEdgeI]++;
1152 edgeList singleEdges(cellEdges.size());
1153 label nSingleEdges = 0;
1157 if (edgeUsage[edgeI] == 1)
1159 singleEdges[nSingleEdges] = cellEdges[edgeI];
1162 else if (edgeUsage[edgeI] != 2)
1171 if (nSingleEdges > 0)
1186 if (
debug || report)
1188 Info<<
" ***Open cells found, number of cells: " << nOpenCells
1189 <<
". This problem may be fixable using the zipUpMesh utility."
1196 if (
debug || report)
1198 Info<<
" Topological cell zip-up check OK." <<
endl;
1219 label nErrorFaces = 0;
1223 const face& curFace =
f[fI];
1240 bool inserted = facePoints.insert(curFace[fp]);
1256 if (nErrorFaces > 0)
1258 if (
debug || report)
1260 Info<<
" Faces with invalid vertex labels found, "
1261 <<
" number of faces: " << nErrorFaces <<
endl;
1267 if (
debug || report)
1269 Info<<
" Face vertices OK." <<
endl;
1287 label nFaceErrors = 0;
1288 label nCellErrors = 0;
1294 if (pf[pointi].empty())
1324 if (nFaceErrors > 0 || nCellErrors > 0)
1326 if (
debug || report)
1328 Info<<
" ***Unused points found in the mesh, "
1329 "number unused by faces: " << nFaceErrors
1330 <<
" number unused by cells: " << nCellErrors
1337 if (
debug || report)
1350 label& nBaffleFaces,
1358 label nbFacei = iter.key();
1359 label nCommon = iter();
1361 const face& curFace = faces()[facei];
1362 const face& nbFace = faces()[nbFacei];
1364 if (nCommon == nbFace.size() || nCommon == curFace.size())
1366 if (nbFace.size() != curFace.size())
1398 label nbFacei = iter.key();
1399 label nCommon = iter();
1401 const face& curFace = faces()[facei];
1402 const face& nbFace = faces()[nbFacei];
1407 && nCommon != nbFace.size()
1408 && nCommon != curFace.size()
1414 label nb = nbFace.find(curFace[fp]);
1428 label fpPlus1 = curFace.fcIndex(fp);
1429 label fpMin1 = curFace.rcIndex(fp);
1432 label nbPlus1 = nbFace.fcIndex(nb);
1433 label nbMin1 = nbFace.rcIndex(nb);
1440 if (nbFace[nbPlus1] == curFace[fpPlus1])
1445 else if (nbFace[nbPlus1] == curFace[fpMin1])
1450 else if (nbFace[nbMin1] == curFace[fpMin1])
1470 if (curFp >= curFace.size())
1476 curFp = curFace.size()-1;
1481 if (curNb >= nbFace.size())
1487 curNb = nbFace.size()-1;
1489 }
while (curFace[curFp] == nbFace[curNb]);
1498 for (
label commonI = 0; commonI < nCommon; commonI++)
1502 if (curFp >= curFace.size())
1508 curFp = curFace.size()-1;
1513 if (curNb >= nbFace.size())
1519 curNb = nbFace.size()-1;
1522 if (curFace[curFp] != nbFace[curNb])
1561 label nBaffleFaces = 0;
1562 label nErrorDuplicate = 0;
1563 label nErrorOrder = 0;
1566 for (
label facei = 0; facei < nFaces(); facei++)
1568 const face& curFace = faces()[facei];
1572 nCommonPoints.clear();
1576 label pointi = curFace[fp];
1582 label nbFacei = nbs[nbI];
1584 if (facei < nbFacei)
1587 ++(nCommonPoints(nbFacei, 0));
1595 if (checkDuplicateFaces(facei, nCommonPoints, nBaffleFaces, setPtr))
1601 if (checkCommonOrder(facei, nCommonPoints, setPtr))
1613 Info<<
" Number of identical duplicate faces (baffle faces): "
1614 << nBaffleFaces <<
endl;
1617 if (nErrorDuplicate > 0 || nErrorOrder > 0)
1620 if (nErrorDuplicate > 0)
1622 Info<<
" <<Number of duplicate (not baffle) faces found: "
1624 <<
". This might indicate a problem." <<
endl;
1627 if (nErrorOrder > 0)
1629 Info<<
" <<Number of faces with non-consecutive shared points: "
1630 << nErrorOrder <<
". This might indicate a problem." <<
endl;
1636 if (
debug || report)
1638 Info<<
" Face-face connectivity OK." <<
endl;
1661 return checkClosedCells
1679 return checkFaceAreas
1695 return checkCellVolumes
1711 return checkFaceOrthogonality
1724 const scalar minPyrVol,
1728 return checkFacePyramids
1746 return checkFaceSkewness
1761 const scalar maxDeg,
1765 return checkFaceAngles
1779 const scalar warnFlatness,
1783 return checkFaceFlatness
1801 return checkConcaveCells
1813 label nFailedChecks = 0;
1815 if (checkPoints(report)) ++nFailedChecks;
1816 if (checkUpperTriangular(report)) ++nFailedChecks;
1817 if (checkCellsZipUp(report)) ++nFailedChecks;
1818 if (checkFaceVertices(report)) ++nFailedChecks;
1819 if (checkFaceFaces(report)) ++nFailedChecks;
1823 if (
debug || report)
1825 Info<<
" Failed " << nFailedChecks
1826 <<
" mesh topology checks." <<
endl;
1832 if (
debug || report)
1834 Info<<
" Mesh topology OK." <<
endl;
1843 label nFailedChecks = 0;
1845 if (checkClosedBoundary(report)) ++nFailedChecks;
1846 if (checkClosedCells(report)) ++nFailedChecks;
1847 if (checkFaceAreas(report)) ++nFailedChecks;
1848 if (checkCellVolumes(report)) ++nFailedChecks;
1849 if (checkFaceOrthogonality(report)) ++nFailedChecks;
1850 if (checkFacePyramids(report)) ++nFailedChecks;
1851 if (checkFaceSkewness(report)) ++nFailedChecks;
1855 if (
debug || report)
1857 Info<<
" Failed " << nFailedChecks
1858 <<
" mesh geometry checks." <<
endl;
1864 if (
debug || report)
1866 Info<<
" Mesh geometry OK." <<
endl;
1884 if (
debug || report)
1886 Info<<
" Failed " << nFailedChecks
1887 <<
" mesh checks." <<
endl;
1893 if (
debug || report)
1904 scalar prev = closedThreshold_;
1905 closedThreshold_ =
val;
1913 scalar prev = aspectThreshold_;
1914 aspectThreshold_ =
val;
1922 scalar prev = nonOrthThreshold_;
1923 nonOrthThreshold_ =
val;
1931 scalar prev = skewThreshold_;
1932 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.
#define InfoInFunction
Report an information message using Foam::Info.
label ListType::const_reference val
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.
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
An Ostream wrapper for parallel output to std::cout.
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.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
static scalar setNonOrthThreshold(const scalar)
Set the non-orthogonality warning threshold in degrees.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
messageStream Info
Information stream (uses stdout - output is on the master only)
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.
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.