44 <<
"Checking whether the boundary is closed" <<
endl;
53 scalar sumMagClosedBoundary = 0;
57 sumClosed += areas[facei];
58 sumMagClosedBoundary +=
mag(areas[facei]);
64 vector openness = sumClosed/(sumMagClosedBoundary + vSmall);
70 Info<<
" ***Boundary openness " << openness
71 <<
" possible hole in boundary description."
81 Info<<
" Boundary openness " << openness <<
" OK."
101 <<
"Checking whether cells are closed" <<
endl;
105 const scalarField& cellVolumes = this->cellVolumes();
110 label nErrorClosed = 0;
114 const cell& curCell =
c[cI];
116 if (
min(curCell) < 0 ||
max(curCell) > nFaces())
127 if (nErrorClosed > 0)
131 Info<<
" ***Cells with invalid face labels found, number of cells "
132 << nErrorClosed <<
endl;
152 scalar maxOpennessCell =
max(openness);
154 scalar maxAspectRatio =
max(aspectRatio);
173 aspectSetPtr->
insert(celli);
191 Info<<
" ***Open cells found, max cell openness: "
192 << maxOpennessCell <<
", number of open cells " << nOpen
203 Info<<
" ***High aspect ratio cells found, Max aspect ratio: "
205 <<
", number of cells " << nAspect
214 Info<<
" Max cell openness = " << maxOpennessCell <<
" OK." <<
nl
215 <<
" Max aspect ratio = " << maxAspectRatio <<
" OK."
237 scalar minArea = great;
238 scalar maxArea = -great;
240 forAll(magFaceAreas, facei)
242 if (magFaceAreas[facei] < vSmall)
250 minArea =
min(minArea, magFaceAreas[facei]);
251 maxArea =
max(maxArea, magFaceAreas[facei]);
257 if (minArea < vSmall)
261 Info<<
" ***Zero or negative face area detected. "
262 "Minimum area: " << minArea <<
endl;
271 Info<<
" Minimum face area = " << minArea
272 <<
". Maximum face area = " << maxArea
273 <<
". Face area magnitudes OK." <<
endl;
293 scalar minVolume = great;
294 scalar maxVolume = -great;
296 label nNegVolCells = 0;
300 if (vols[celli] < vSmall)
310 minVolume =
min(minVolume, vols[celli]);
311 maxVolume =
max(maxVolume, vols[celli]);
318 if (minVolume < vSmall)
322 Info<<
" ***Zero or negative cell volume detected. "
323 <<
"Minimum negative volume: " << minVolume
324 <<
", Number of negative volume cells: " << nNegVolCells
334 Info<<
" Min volume = " << minVolume
335 <<
". Max volume = " << maxVolume
336 <<
". Total volume = " <<
gSum(vols)
337 <<
". Cell volumes OK." <<
endl;
348 const scalar minPyrVol,
373 label nErrorPyrs = 0;
377 if (ownPyrVol[facei] < minPyrVol)
387 if (isInternalFace(facei))
389 if (neiPyrVol[facei] < minPyrVol)
406 Info<<
" ***Error in face pyramids: "
407 << nErrorPyrs <<
" faces are incorrectly oriented."
417 Info<<
" Face pyramids OK." <<
endl;
437 if (maxDeg < -small || maxDeg > 180+small)
440 <<
"maxDeg should be [0..180] but is now " << maxDeg
458 scalar maxEdgeSin =
max(faceAngles);
464 if (faceAngles[facei] > small)
480 scalar maxConcaveDegr =
485 Info<<
" *There are " << nConcave
486 <<
" faces with concave angles between consecutive"
487 <<
" edges. Max concave angle = " << maxConcaveDegr
488 <<
" degrees." <<
endl;
497 Info<<
" All angles in faces OK." <<
endl;
508 const scalar warnFlatness,
517 if (warnFlatness < 0 || warnFlatness > 1)
520 <<
"warnFlatness should be [0..1] but is now " << warnFlatness
525 const vectorField& faceCentres = this->faceCentres();
540 scalar minFlatness = great;
541 scalar sumFlatness = 0;
545 forAll(faceFlatness, facei)
547 if (fcs[facei].size() > 3 && magAreas[facei] > vSmall)
549 sumFlatness += faceFlatness[facei];
552 minFlatness =
min(minFlatness, faceFlatness[facei]);
554 if (faceFlatness[facei] < warnFlatness)
577 Info<<
" Face flatness (1 = flat, 0 = butterfly) : min = "
578 << minFlatness <<
" average = " << sumFlatness / nSummed
588 Info<<
" *There are " << nWarped
589 <<
" faces with ratio between projected and actual area < "
590 << warnFlatness <<
endl;
592 Info<<
" Minimum ratio (minimum flatness, maximum warpage) = "
593 << minFlatness <<
endl;
602 Info<<
" All face flatness OK." <<
endl;
626 label nConcaveCells = 0;
630 const cell& cFaces =
c[celli];
632 bool concave =
false;
641 label fI = cFaces[i];
643 const point& fC = fCentres[fI];
647 fN /=
max(
mag(fN), vSmall);
651 if (fOwner[fI] != celli)
663 label fJ = cFaces[j];
665 const point& pt = fCentres[fJ];
675 pC /=
max(
mag(pC), vSmall);
699 if (nConcaveCells > 0)
703 Info<<
" ***Concave cells (using face planes) found,"
704 <<
" number of cells: " << nConcaveCells <<
endl;
713 Info<<
" Concave cell check OK." <<
endl;
740 label internal = nInternalFaces();
745 label nMultipleCells =
false;
749 for (
label facei = 0; facei <
internal; facei++)
751 if (own[facei] >= nei[facei])
775 label facei = curFaces[i];
777 if (facei >= nInternalFaces())
784 label nbrCelli = nei[facei];
786 if (nbrCelli == celli)
788 nbrCelli = own[facei];
791 if (celli < nbrCelli)
811 label prevCell = nbr[0];
814 bool hasMultipleFaces =
false;
818 label thisCell = nbr[i];
826 if (thisCell == prevCell)
828 hasMultipleFaces =
true;
836 else if (thisFace < prevFace)
850 if (hasMultipleFaces)
859 if ((debug || report) && nMultipleCells > 0)
861 Info<<
" <<Found " << nMultipleCells
862 <<
" neighbouring cells with multiple in between faces." <<
endl;
869 Info<<
" ***Faces not in upper triangular order." <<
endl;
878 Info<<
" Upper triangular ordering OK." <<
endl;
897 label nOpenCells = 0;
912 edgeList curFaceEdges =
f[curFaces[facei]].edges();
914 forAll(curFaceEdges, faceEdgeI)
916 const edge& curEdge = curFaceEdges[faceEdgeI];
918 forAll(cellEdges, cellEdgeI)
920 if (cellEdges[cellEdgeI] == curEdge)
922 edgeUsage[cellEdgeI]++;
930 label nSingleEdges = 0;
934 if (edgeUsage[edgeI] == 1)
936 singleEdges[nSingleEdges] = cellEdges[edgeI];
939 else if (edgeUsage[edgeI] != 2)
948 if (nSingleEdges > 0)
965 Info<<
" ***Open cells found, number of cells: " << nOpenCells
966 <<
". This problem may be fixable using the zipUpMesh utility."
976 Info<<
" Topological cell zip-up check OK." <<
endl;
998 label nErrorFaces = 0;
1002 const face& curFace =
f[fI];
1019 bool inserted = facePoints.
insert(curFace[fp]);
1035 if (nErrorFaces > 0)
1037 if (debug || report)
1039 Info<<
" Faces with invalid vertex labels found, "
1040 <<
" number of faces: " << nErrorFaces <<
endl;
1047 if (debug || report)
1049 Info<<
" Face vertices OK." <<
endl;
1068 label nFaceErrors = 0;
1069 label nCellErrors = 0;
1075 if (pf[pointi].empty())
1089 const labelList& pc = pointCells(pointi);
1105 if (nFaceErrors > 0 || nCellErrors > 0)
1107 if (debug || report)
1109 Info<<
" ***Unused points found in the mesh, "
1110 "number unused by faces: " << nFaceErrors
1111 <<
" number unused by cells: " << nCellErrors
1119 if (debug || report)
1133 label& nBaffleFaces,
1141 label nbFacei = iter.key();
1142 label nCommon = iter();
1144 const face& curFace = faces()[facei];
1145 const face& nbFace = faces()[nbFacei];
1147 if (nCommon == nbFace.
size() || nCommon == curFace.
size())
1149 if (nbFace.
size() != curFace.
size())
1181 label nbFacei = iter.key();
1182 label nCommon = iter();
1184 const face& curFace = faces()[facei];
1185 const face& nbFace = faces()[nbFacei];
1190 && nCommon != nbFace.
size()
1191 && nCommon != curFace.
size()
1223 if (nbFace[nbPlus1] == curFace[fpPlus1])
1228 else if (nbFace[nbPlus1] == curFace[fpMin1])
1233 else if (nbFace[nbMin1] == curFace[fpMin1])
1253 if (curFp >= curFace.
size())
1259 curFp = curFace.
size()-1;
1264 if (curNb >= nbFace.
size())
1270 curNb = nbFace.
size()-1;
1272 }
while (curFace[curFp] == nbFace[curNb]);
1281 for (
label commonI = 0; commonI < nCommon; commonI++)
1285 if (curFp >= curFace.
size())
1291 curFp = curFace.
size()-1;
1296 if (curNb >= nbFace.
size())
1302 curNb = nbFace.
size()-1;
1305 if (curFace[curFp] != nbFace[curNb])
1344 label nBaffleFaces = 0;
1345 label nErrorDuplicate = 0;
1346 label nErrorOrder = 0;
1349 for (
label facei = 0; facei < nFaces(); facei++)
1351 const face& curFace = faces()[facei];
1355 nCommonPoints.
clear();
1359 label pointi = curFace[fp];
1365 label nbFacei = nbs[nbI];
1367 if (facei < nbFacei)
1373 if (fnd == nCommonPoints.
end())
1376 nCommonPoints.
insert(nbFacei, 1);
1389 if (checkDuplicateFaces(facei, nCommonPoints, nBaffleFaces, setPtr))
1395 if (checkCommonOrder(facei, nCommonPoints, setPtr))
1407 Info<<
" Number of identical duplicate faces (baffle faces): "
1408 << nBaffleFaces <<
endl;
1411 if (nErrorDuplicate > 0 || nErrorOrder > 0)
1414 if (nErrorDuplicate > 0)
1416 Info<<
" <<Number of duplicate (not baffle) faces found: "
1418 <<
". This might indicate a problem." <<
endl;
1421 if (nErrorOrder > 0)
1423 Info<<
" <<Number of faces with non-consecutive shared points: "
1424 << nErrorOrder <<
". This might indicate a problem." <<
endl;
1431 if (debug || report)
1433 Info<<
" Face-face connectivity OK." <<
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.
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 checkClosedCells(const bool report=false, labelHashSet *setPtr=nullptr, labelHashSet *highAspectSetPtr=nullptr, const Vector< label > &solutionD=Vector< label >::one) const
Check cells for closedness.
bool checkFacePyramids(const bool report=false, const scalar minPyrVol=-small, labelHashSet *setPtr=nullptr) const
Check face pyramid volume.
bool checkCellsZipUp(const bool report=false, labelHashSet *setPtr=nullptr) const
Check cell zip-up.
bool checkFaceAreas(const bool report=false, labelHashSet *setPtr=nullptr) const
Check for negative face areas.
bool checkCommonOrder(const label, const Map< label > &, labelHashSet *) const
Check that shared points are in consecutive order.
bool checkFaceAngles(const bool report=false, const scalar maxSin=10, labelHashSet *setPtr=nullptr) const
Check face angles.
bool checkUpperTriangular(const bool report=false, labelHashSet *setPtr=nullptr) const
Check face ordering.
bool checkConcaveCells(const bool report=false, labelHashSet *setPtr=nullptr) const
Check for concave cells by the planes of faces.
bool checkClosedBoundary(const bool report=false) const
Check boundary for closedness.
bool checkCellVolumes(const bool report=false, labelHashSet *setPtr=nullptr) const
Check for negative cell volumes.
bool checkDuplicateFaces(const label, const Map< label > &, label &nBaffleFaces, labelHashSet *) const
Check if all points on face are shared with another face.
bool checkFaceVertices(const bool report=false, labelHashSet *setPtr=nullptr) const
Check uniqueness of face vertices.
label nInternalFaces() const
bool checkPoints(const bool report=false, labelHashSet *setPtr=nullptr) const
Check for unused points.
const vectorField & faceAreas() const
bool checkFaceFlatness(const bool report, const scalar warnFlatness, labelHashSet *setPtr) const
Check face warpage: decompose face and check ratio between.
bool checkFaceFaces(const bool report=false, labelHashSet *setPtr=nullptr) const
Check face-face connectivity.
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.
scalar closedThreshold
Data to control mesh checking.
scalar aspectThreshold
Aspect ratio warning threshold.
scalar planarCosAngle
Threshold where faces are considered coplanar.
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.
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
dimensionedScalar sin(const dimensionedScalar &ds)
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)
dimensionSet cmptMag(const dimensionSet &)
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
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
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.