55 <<
"Checking whether the boundary is closed" <<
endl;
62 scalar sumMagClosedBoundary = 0;
64 for (
label facei = nInternalFaces(); facei < areas.
size(); facei++)
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;
152 primitiveMeshTools::cellClosedness
163 scalar maxOpennessCell =
max(openness);
165 scalar maxAspectRatio =
max(aspectRatio);
170 if (openness[celli] > closedThreshold_)
180 if (aspectRatio[celli] > aspectThreshold_)
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,
514 primitiveMeshTools::facePyramidVolume
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);
629 if (skewness[facei] > skewThreshold_)
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);
919 if ((pC & fN) > -planarCosAngle_)
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];
1054 label prevFace = curFaces[nbr.indices()[0]];
1056 bool hasMultipleFaces =
false;
1058 for (
label i = 1; i < nbr.size(); i++)
1060 label thisCell = nbr[i];
1061 label thisFace = curFaces[nbr.indices()[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;
1762 const scalar minPyrVol,
1799 const scalar maxDeg,
1817 const scalar warnFlatness,
1851 label noFailedChecks = 0;
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;
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;
bool checkFaceAreas(const vectorField &faceAreas, const bool report, const bool detailedReport, labelHashSet *setPtr) const
Check for negative face areas.
static scalar aspectThreshold_
Aspect ratio warning threshold.
dimensionedScalar acos(const dimensionedScalar &ds)
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
#define forAll(list, i)
Loop across all elements in list.
bool checkFaceAngles(const pointField &points, const vectorField &faceAreas, const bool report, const scalar maxDeg, labelHashSet *setPtr) const
Check face angles.
bool empty() const
Return true if the UList is empty (ie, size() is zero)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
errorManipArg< error, int > exit(error &err, const int errNo=1)
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
A face is a list of labels corresponding to mesh vertices.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
A list that is sorted upon construction or when explicitly requested with the sort() method...
Unit conversion functions.
void size(const label)
Override size to be inconsistent with allocated storage.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
static scalar setClosedThreshold(const scalar)
Set the closedness ratio warning threshold.
Ostream & endl(Ostream &os)
Add newline and flush stream.
virtual const pointField & points() const =0
Return mesh points.
bool insert(const Key &key)
Insert a new entry.
dimensionedScalar asin(const dimensionedScalar &ds)
static scalar setAspectThreshold(const scalar)
Set the aspect ratio warning threshold.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Various functions to operate on Lists.
static scalar setSkewThreshold(const scalar)
Set the skewness warning threshold as percentage.
Class to handle errors and exceptions in a simple, consistent stream-based manner.
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
bool checkFaceSkewness(const pointField &points, const vectorField &fCtrs, const vectorField &fAreas, const vectorField &cellCtrs, const bool report, labelHashSet *setPtr) const
Check face skewness.
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
bool checkClosedCells(const vectorField &faceAreas, const scalarField &cellVolumes, const bool report, labelHashSet *setPtr, labelHashSet *aspectSetPtr, const Vector< label > &meshD) const
Check cells for closedness.
bool checkFacePyramids(const pointField &points, const vectorField &ctrs, const bool report, const bool detailedReport, const scalar minPyrVol, labelHashSet *setPtr) const
Check face pyramid volume.
Type gSum(const FieldField< Field, Type > &f)
dimensionedScalar cos(const dimensionedScalar &ds)
bool checkFaceOrthogonality(const vectorField &fAreas, const vectorField &cellCtrs, const bool report, labelHashSet *setPtr) const
Check for non-orthogonality.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
static scalar closedThreshold_
Static data to control mesh checking.
void clear()
Clear all entries from table.
virtual bool checkMesh(const bool report=false) const
Check mesh for correctness. Returns false for no error.
static const label labelMax
static scalar setNonOrthThreshold(const scalar)
Set the non-orthogonality warning threshold in degrees.
const vectorField & cellCentres() const
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
bool checkCellVolumes(const scalarField &vols, const bool report, const bool detailedReport, labelHashSet *setPtr) const
Check for negative cell volumes.
static scalar planarCosAngle_
Threshold where faces are considered coplanar.
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
dimensionedScalar sin(const dimensionedScalar &ds)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
virtual bool checkFaceVertices(const bool report=false, labelHashSet *setPtr=nullptr) const
Check uniqueness of face vertices.
const vectorField & faceCentres() const
virtual bool checkTopology(const bool report=false) const
Check mesh topology for correctness.
bool checkConcaveCells(const vectorField &fAreas, const pointField &fCentres, const bool report, labelHashSet *setPtr) const
Check for concave cells by the planes of faces.
A cell is defined as a list of faces with extra functionality.
prefixOSstream Pout(cout, "Pout")
bool checkCommonOrder(const label, const Map< label > &, labelHashSet *) const
Check that shared points are in consecutive order.
const vectorField & faceAreas() const
const dimensionedScalar c
Speed of light in a vacuum.
bool checkDuplicateFaces(const label, const Map< label > &, label &nBaffleFaces, labelHashSet *) const
Check if all points on face are shared with another face.
virtual bool checkCellsZipUp(const bool report=false, labelHashSet *setPtr=nullptr) const
Check cell zip-up.
virtual bool checkUpperTriangular(const bool report=false, labelHashSet *setPtr=nullptr) const
Check face ordering.
virtual bool checkPoints(const bool report=false, labelHashSet *setPtr=nullptr) const
Check for unused points.
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual bool checkGeometry(const bool report=false) const
Check mesh geometry (& implicitly topology) for correctness.
bool checkFaceFlatness(const pointField &points, const vectorField &faceCentres, const vectorField &faceAreas, const bool report, const scalar warnFlatness, labelHashSet *setPtr) const
Check face warpage.
label size() const
Number of entries.
virtual bool checkFaceFaces(const bool report=false, labelHashSet *setPtr=nullptr) const
Check face-face connectivity.
static scalar skewThreshold_
Skewness warning threshold.
static scalar nonOrthThreshold_
Non-orthogonality warning threshold in deg.
A class for managing temporary objects.
bool checkClosedBoundary(const vectorField &, const bool, const PackedBoolList &) const
Check boundary for closedness.
const scalarField & cellVolumes() const
#define InfoInFunction
Report an information message using Foam::Info.