55 <<
"Checking whether the boundary is closed" <<
endl;
62 scalar sumMagClosedBoundary = 0;
66 if (!internalOrCoupledFaces.
size() || !internalOrCoupledFaces[facei])
68 sumClosed += areas[facei];
69 sumMagClosedBoundary +=
mag(areas[facei]);
76 vector openness = sumClosed/(sumMagClosedBoundary + vSmall);
82 Info<<
" ***Boundary openness " << openness
83 <<
" possible hole in boundary description."
93 Info<<
" Boundary openness " << openness <<
" OK."
115 <<
"Checking whether cells are closed" <<
endl;
121 label nErrorClosed = 0;
125 const cell& curCell =
c[cI];
127 if (
min(curCell) < 0 ||
max(curCell) > nFaces())
138 if (nErrorClosed > 0)
142 Info<<
" ***Cells with invalid face labels found, number of cells "
143 << nErrorClosed <<
endl;
163 scalar maxOpennessCell =
max(openness);
165 scalar maxAspectRatio =
max(aspectRatio);
184 aspectSetPtr->
insert(celli);
202 Info<<
" ***Open cells found, max cell openness: "
203 << maxOpennessCell <<
", number of open cells " << nOpen
214 Info<<
" ***High aspect ratio cells found, Max aspect ratio: "
216 <<
", number of cells " << nAspect
225 Info<<
" Max cell openness = " << maxOpennessCell <<
" OK." <<
nl
226 <<
" Max aspect ratio = " << maxAspectRatio <<
" OK."
238 const bool detailedReport,
249 scalar minArea = great;
250 scalar maxArea = -great;
252 forAll(magFaceAreas, facei)
254 if (magFaceAreas[facei] < vSmall)
262 if (isInternalFace(facei))
264 Pout<<
"Zero or negative face area detected for "
265 <<
"internal face "<< facei <<
" between cells "
266 << faceOwner()[facei] <<
" and "
267 << faceNeighbour()[facei]
268 <<
". Face area magnitude = " << magFaceAreas[facei]
273 Pout<<
"Zero or negative face area detected for "
274 <<
"boundary face " << facei <<
" next to cell "
275 << faceOwner()[facei] <<
". Face area magnitude = "
276 << magFaceAreas[facei] <<
endl;
281 minArea =
min(minArea, magFaceAreas[facei]);
282 maxArea =
max(maxArea, magFaceAreas[facei]);
288 if (minArea < vSmall)
292 Info<<
" ***Zero or negative face area detected. "
293 "Minimum area: " << minArea <<
endl;
302 Info<<
" Minimum face area = " << minArea
303 <<
". Maximum face area = " << maxArea
304 <<
". Face area magnitudes OK." <<
endl;
316 const bool detailedReport,
325 scalar minVolume = great;
326 scalar maxVolume = -great;
328 label nNegVolCells = 0;
332 if (vols[celli] < vSmall)
340 Pout<<
"Zero or negative cell volume detected for cell "
341 << celli <<
". Volume = " << vols[celli] <<
endl;
347 minVolume =
min(minVolume, vols[celli]);
348 maxVolume =
max(maxVolume, vols[celli]);
355 if (minVolume < vSmall)
359 Info<<
" ***Zero or negative cell volume detected. "
360 <<
"Minimum negative volume: " << minVolume
361 <<
", Number of negative volume cells: " << nNegVolCells
371 Info<<
" Min volume = " << minVolume
372 <<
". Max volume = " << maxVolume
373 <<
". Total volume = " <<
gSum(vols)
374 <<
". Cell volumes OK." <<
endl;
405 const scalar severeNonorthogonalityThreshold =
408 scalar minDDotS =
min(ortho);
410 scalar sumDDotS =
sum(ortho);
412 label severeNonOrth = 0;
414 label errorNonOrth = 0;
419 if (ortho[facei] < severeNonorthogonalityThreshold)
421 if (ortho[facei] > small)
456 Info<<
" Mesh non-orthogonality Max: "
463 if (severeNonOrth > 0)
465 Info<<
" *Number of severely non-orthogonal faces: "
466 << severeNonOrth <<
"." <<
endl;
470 if (errorNonOrth > 0)
474 Info<<
" ***Number of non-orthogonality errors: "
475 << errorNonOrth <<
"." <<
endl;
484 Info<<
" Non-orthogonality check OK." <<
endl;
497 const bool detailedReport,
498 const scalar minPyrVol,
524 label nErrorPyrs = 0;
528 if (ownPyrVol[facei] < minPyrVol)
536 Pout<<
"Negative pyramid volume: " << ownPyrVol[facei]
537 <<
" for face " << facei <<
" " <<
f[facei]
538 <<
" and owner cell: " << own[facei] <<
endl
539 <<
"Owner cell vertex labels: "
540 <<
cells()[own[facei]].labels(faces())
547 if (isInternalFace(facei))
549 if (neiPyrVol[facei] < minPyrVol)
557 Pout<<
"Negative pyramid volume: " << neiPyrVol[facei]
558 <<
" for face " << facei <<
" " <<
f[facei]
559 <<
" and neighbour cell: " << nei[facei] <<
nl
560 <<
"Neighbour cell vertex labels: "
561 <<
cells()[nei[facei]].labels(faces())
575 Info<<
" ***Error in face pyramids: "
576 << nErrorPyrs <<
" faces are incorrectly oriented."
586 Info<<
" Face pyramids OK." <<
endl;
622 scalar maxSkew =
max(skewness);
647 Info<<
" ***Max skewness = " << maxSkew
648 <<
", " << nWarnSkew <<
" highly skew faces detected"
649 " which may impair the quality of the results"
659 Info<<
" Max skewness = " << maxSkew <<
" OK." <<
endl;
681 if (maxDeg < -small || maxDeg > 180+small)
684 <<
"maxDeg should be [0..180] but is now " << maxDeg
700 scalar maxEdgeSin =
max(faceAngles);
706 if (faceAngles[facei] > small)
722 scalar maxConcaveDegr =
727 Info<<
" *There are " << nConcave
728 <<
" faces with concave angles between consecutive"
729 <<
" edges. Max concave angle = " << maxConcaveDegr
730 <<
" degrees." <<
endl;
739 Info<<
" All angles in faces OK." <<
endl;
753 const scalar warnFlatness,
762 if (warnFlatness < 0 || warnFlatness > 1)
765 <<
"warnFlatness should be [0..1] but is now " << warnFlatness
782 scalar minFlatness = great;
783 scalar sumFlatness = 0;
787 forAll(faceFlatness, facei)
789 if (fcs[facei].size() > 3 && magAreas[facei] > vSmall)
791 sumFlatness += faceFlatness[facei];
794 minFlatness =
min(minFlatness, faceFlatness[facei]);
796 if (faceFlatness[facei] < warnFlatness)
819 Info<<
" Face flatness (1 = flat, 0 = butterfly) : min = "
820 << minFlatness <<
" average = " << sumFlatness / nSummed
830 Info<<
" *There are " << nWarped
831 <<
" faces with ratio between projected and actual area < "
832 << warnFlatness <<
endl;
834 Info<<
" Minimum ratio (minimum flatness, maximum warpage) = "
835 << minFlatness <<
endl;
844 Info<<
" All face flatness OK." <<
endl;
868 label nConcaveCells = 0;
872 const cell& cFaces =
c[celli];
874 bool concave =
false;
883 label fI = cFaces[i];
885 const point& fC = fCentres[fI];
889 fN /=
max(
mag(fN), vSmall);
893 if (fOwner[fI] != celli)
905 label fJ = cFaces[j];
907 const point& pt = fCentres[fJ];
917 pC /=
max(
mag(pC), vSmall);
941 if (nConcaveCells > 0)
945 Info<<
" ***Concave cells (using face planes) found,"
946 <<
" number of cells: " << nConcaveCells <<
endl;
955 Info<<
" Concave cell check OK." <<
endl;
982 label internal = nInternalFaces();
987 label nMultipleCells =
false;
991 for (
label facei = 0; facei <
internal; facei++)
993 if (own[facei] >= nei[facei])
1017 label facei = curFaces[i];
1019 if (facei >= nInternalFaces())
1026 label nbrCelli = nei[facei];
1028 if (nbrCelli == celli)
1030 nbrCelli = own[facei];
1033 if (celli < nbrCelli)
1053 label prevCell = nbr[0];
1056 bool hasMultipleFaces =
false;
1060 label thisCell = nbr[i];
1068 if (thisCell == prevCell)
1070 hasMultipleFaces =
true;
1074 setPtr->
insert(prevFace);
1075 setPtr->
insert(thisFace);
1078 else if (thisFace < prevFace)
1084 setPtr->
insert(thisFace);
1088 prevCell = thisCell;
1089 prevFace = thisFace;
1092 if (hasMultipleFaces)
1101 if ((debug || report) && nMultipleCells > 0)
1103 Info<<
" <<Found " << nMultipleCells
1104 <<
" neighbouring cells with multiple in between faces." <<
endl;
1109 if (debug || report)
1111 Info<<
" ***Faces not in upper triangular order." <<
endl;
1118 if (debug || report)
1120 Info<<
" Upper triangular ordering OK." <<
endl;
1139 label nOpenCells = 0;
1148 const edgeList cellEdges =
c[celli].edges(
f);
1154 edgeList curFaceEdges =
f[curFaces[facei]].edges();
1156 forAll(curFaceEdges, faceEdgeI)
1158 const edge& curEdge = curFaceEdges[faceEdgeI];
1160 forAll(cellEdges, cellEdgeI)
1162 if (cellEdges[cellEdgeI] == curEdge)
1164 edgeUsage[cellEdgeI]++;
1172 label nSingleEdges = 0;
1176 if (edgeUsage[edgeI] == 1)
1178 singleEdges[nSingleEdges] = cellEdges[edgeI];
1181 else if (edgeUsage[edgeI] != 2)
1190 if (nSingleEdges > 0)
1205 if (debug || report)
1207 Info<<
" ***Open cells found, number of cells: " << nOpenCells
1208 <<
". This problem may be fixable using the zipUpMesh utility."
1216 if (debug || report)
1218 Info<<
" Topological cell zip-up check OK." <<
endl;
1240 label nErrorFaces = 0;
1244 const face& curFace =
f[fI];
1261 bool inserted = facePoints.
insert(curFace[fp]);
1277 if (nErrorFaces > 0)
1279 if (debug || report)
1281 Info<<
" Faces with invalid vertex labels found, "
1282 <<
" number of faces: " << nErrorFaces <<
endl;
1289 if (debug || report)
1291 Info<<
" Face vertices OK." <<
endl;
1310 label nFaceErrors = 0;
1311 label nCellErrors = 0;
1317 if (pf[pointi].empty())
1331 const labelList& pc = pointCells(pointi);
1347 if (nFaceErrors > 0 || nCellErrors > 0)
1349 if (debug || report)
1351 Info<<
" ***Unused points found in the mesh, "
1352 "number unused by faces: " << nFaceErrors
1353 <<
" number unused by cells: " << nCellErrors
1361 if (debug || report)
1375 label& nBaffleFaces,
1383 label nbFacei = iter.key();
1384 label nCommon = iter();
1386 const face& curFace = faces()[facei];
1387 const face& nbFace = faces()[nbFacei];
1389 if (nCommon == nbFace.
size() || nCommon == curFace.
size())
1391 if (nbFace.
size() != curFace.
size())
1423 label nbFacei = iter.key();
1424 label nCommon = iter();
1426 const face& curFace = faces()[facei];
1427 const face& nbFace = faces()[nbFacei];
1432 && nCommon != nbFace.
size()
1433 && nCommon != curFace.
size()
1465 if (nbFace[nbPlus1] == curFace[fpPlus1])
1470 else if (nbFace[nbPlus1] == curFace[fpMin1])
1475 else if (nbFace[nbMin1] == curFace[fpMin1])
1495 if (curFp >= curFace.
size())
1501 curFp = curFace.
size()-1;
1506 if (curNb >= nbFace.
size())
1512 curNb = nbFace.
size()-1;
1514 }
while (curFace[curFp] == nbFace[curNb]);
1523 for (
label commonI = 0; commonI < nCommon; commonI++)
1527 if (curFp >= curFace.
size())
1533 curFp = curFace.
size()-1;
1538 if (curNb >= nbFace.
size())
1544 curNb = nbFace.
size()-1;
1547 if (curFace[curFp] != nbFace[curNb])
1586 label nBaffleFaces = 0;
1587 label nErrorDuplicate = 0;
1588 label nErrorOrder = 0;
1591 for (
label facei = 0; facei < nFaces(); facei++)
1593 const face& curFace = faces()[facei];
1597 nCommonPoints.
clear();
1601 label pointi = curFace[fp];
1607 label nbFacei = nbs[nbI];
1609 if (facei < nbFacei)
1615 if (fnd == nCommonPoints.
end())
1618 nCommonPoints.
insert(nbFacei, 1);
1631 if (checkDuplicateFaces(facei, nCommonPoints, nBaffleFaces, setPtr))
1637 if (checkCommonOrder(facei, nCommonPoints, setPtr))
1649 Info<<
" Number of identical duplicate faces (baffle faces): "
1650 << nBaffleFaces <<
endl;
1653 if (nErrorDuplicate > 0 || nErrorOrder > 0)
1656 if (nErrorDuplicate > 0)
1658 Info<<
" <<Number of duplicate (not baffle) faces found: "
1660 <<
". This might indicate a problem." <<
endl;
1663 if (nErrorOrder > 0)
1665 Info<<
" <<Number of faces with non-consecutive shared points: "
1666 << nErrorOrder <<
". This might indicate a problem." <<
endl;
1673 if (debug || report)
1675 Info<<
" Face-face connectivity OK." <<
endl;
1687 return checkClosedBoundary(faceAreas(), report,
PackedBoolList(0));
1699 return checkClosedCells
1717 return checkFaceAreas
1733 return checkCellVolumes
1749 return checkFaceOrthogonality
1762 const scalar minPyrVol,
1799 const scalar maxDeg,
1817 const scalar warnFlatness,
1839 return checkConcaveCells
1851 label noFailedChecks = 0;
1853 if (checkPoints(report)) noFailedChecks++;
1854 if (checkUpperTriangular(report)) noFailedChecks++;
1855 if (checkCellsZipUp(report)) noFailedChecks++;
1856 if (checkFaceVertices(report)) noFailedChecks++;
1857 if (checkFaceFaces(report)) noFailedChecks++;
1859 if (noFailedChecks == 0)
1861 if (debug || report)
1863 Info<<
" Mesh topology OK." <<
endl;
1870 if (debug || report)
1872 Info<<
" Failed " << noFailedChecks
1873 <<
" mesh topology checks." <<
endl;
1883 label noFailedChecks = 0;
1885 if (checkClosedBoundary(report)) noFailedChecks++;
1886 if (checkClosedCells(report)) noFailedChecks++;
1887 if (checkFaceAreas(report)) noFailedChecks++;
1888 if (checkCellVolumes(report)) noFailedChecks++;
1889 if (checkFaceOrthogonality(report)) noFailedChecks++;
1893 if (noFailedChecks == 0)
1895 if (debug || report)
1897 Info<<
" Mesh geometry OK." <<
endl;
1903 if (debug || report)
1905 Info<<
" Failed " << noFailedChecks
1906 <<
" mesh geometry checks." <<
endl;
1923 if (noFailedChecks == 0)
1925 if (debug || report)
1934 if (debug || report)
1936 Info<<
" Failed " << noFailedChecks
1937 <<
" mesh checks." <<
endl;
Various functions to operate on Lists.
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
bool insert(const Key &key)
Insert a new entry.
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
void clear()
Clear all entries from table.
void size(const label)
Override size to be inconsistent with allocated storage.
label size() const
Number of entries.
A list that is sorted upon construction or when explicitly requested with the sort() method.
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
void sort()
(stable) sort the list (if changed after construction time)
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
bool empty() const
Return true if the UList is empty (ie, size() is zero)
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.
bool checkFaceOrthogonality(const vectorField &fAreas, const vectorField &cellCtrs, const bool report, labelHashSet *setPtr) const
Check for non-orthogonality.
virtual bool checkCellsZipUp(const bool report=false, labelHashSet *setPtr=nullptr) const
Check cell zip-up.
virtual bool checkMesh(const bool report=false) const
Check mesh for correctness. Returns false for no error.
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.
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.
bool checkClosedBoundary(const vectorField &, const bool, const PackedBoolList &) const
Check boundary for closedness.
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.
label nInternalFaces() const
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.
bool checkConcaveCells(const vectorField &fAreas, const pointField &fCentres, const bool report, labelHashSet *setPtr) const
Check for concave cells by the planes of faces.
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.
A class for managing temporary objects.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
#define InfoInFunction
Report an information message using Foam::Info.
const dimensionedScalar c
Speed of light in a vacuum.
bool checkFacePyramids(const bool report, const scalar minPyrVol, const polyMesh &, const vectorField &cellCentres, const pointField &p, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *)
Check face pyramid volumes.
bool checkFaceSkewness(const bool report, const scalar internalSkew, const scalar boundarySkew, const polyMesh &mesh, const pointField &points, const vectorField &cellCentres, const vectorField &faceCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
Check face skewness.
bool checkFaceAngles(const bool report, const scalar maxDeg, const polyMesh &mesh, const vectorField &faceAreas, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Check convexity of angles in a face. See primitiveMesh for explanation.
bool checkFaceFlatness(const bool report, const scalar minFlatness, const polyMesh &, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Check for face areas v.s. sum of face-triangle (from face-centre.
scalar closedThreshold
Data to control mesh checking.
scalar skewThreshold
Skewness warning threshold.
scalar aspectThreshold
Aspect ratio warning threshold.
scalar planarCosAngle
Threshold where faces are considered coplanar.
scalar nonOrthThreshold
Non-orthogonality warning threshold in deg.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensionedScalar asin(const dimensionedScalar &ds)
Type gSum(const FieldField< Field, Type > &f)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
dimensionedScalar sin(const dimensionedScalar &ds)
label checkGeometry(const polyMesh &mesh, const bool allGeometry, const autoPtr< surfaceWriter > &, const autoPtr< setWriter > &)
Check the geometry.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
label checkTopology(const polyMesh &mesh, const bool allTopology, const bool allGeometry, const autoPtr< surfaceWriter > &surfWriter, const autoPtr< Foam::setWriter > &setWriter)
dimensionSet cmptMag(const dimensionSet &)
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
prefixOSstream Pout(cout, "Pout")
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
static const label labelMax
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar acos(const dimensionedScalar &ds)
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Unit conversion functions.