35 Foam::scalar Foam::primitiveMesh::closedThreshold_ = 1.0e-6;
36 Foam::scalar Foam::primitiveMesh::aspectThreshold_ = 1000;
37 Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70;
38 Foam::scalar Foam::primitiveMesh::skewThreshold_ = 4;
47 Info<<
"bool primitiveMesh::checkClosedBoundary("
48 <<
"const bool) const: "
49 <<
"checking whether the boundary is closed" <<
endl;
56 scalar sumMagClosedBoundary = 0;
62 sumClosed += areas[faceI];
63 sumMagClosedBoundary +=
mag(areas[faceI]);
69 vector openness = sumClosed/(sumMagClosedBoundary + VSMALL);
75 Info<<
" ***Boundary openness " << openness
76 <<
" possible hole in boundary description."
86 Info<<
" Boundary openness " << openness <<
" OK."
105 Info<<
"bool primitiveMesh::checkClosedCells("
106 <<
"const bool, labelHashSet*, labelHashSet*"
107 <<
", const Vector<label>&) const: "
108 <<
"checking whether cells are closed" <<
endl;
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;
155 sumClosed[own[faceI]] += areas[faceI];
156 sumMagClosed[own[faceI]] +=
cmptMag(areas[faceI]);
162 sumClosed[nei[faceI]] -= areas[faceI];
163 sumMagClosed[nei[faceI]] +=
cmptMag(areas[faceI]);
167 scalar maxOpennessCell = 0;
170 scalar maxAspectRatio = 0;
187 scalar maxOpenness = 0;
194 mag(sumClosed[cellI][cmpt])
195 /(sumMagClosed[cellI][cmpt] + VSMALL)
199 maxOpennessCell =
max(maxOpennessCell, maxOpenness);
201 if (maxOpenness > closedThreshold_)
213 scalar minCmpt = VGREAT;
214 scalar maxCmpt = -VGREAT;
219 minCmpt =
min(minCmpt, sumMagClosed[cellI][dir]);
220 maxCmpt =
max(maxCmpt, sumMagClosed[cellI][dir]);
224 scalar aspectRatio = maxCmpt/(minCmpt + VSMALL);
230 1.0/6.0*
cmptSum(sumMagClosed[cellI])/
pow(vols[cellI], 2.0/3.0)
234 maxAspectRatio =
max(maxAspectRatio, aspectRatio);
236 if (aspectRatio > aspectThreshold_)
240 aspectSetPtr->
insert(cellI);
258 Info<<
" ***Open cells found, max cell openness: "
259 << maxOpennessCell <<
", number of open cells " << nOpen
270 Info<<
" ***High aspect ratio cells found, Max aspect ratio: "
272 <<
", number of cells " << nAspect
281 Info<<
" Max cell openness = " << maxOpennessCell <<
" OK." <<
nl
282 <<
" Max aspect ratio = " << maxAspectRatio <<
" OK."
298 Info<<
"bool primitiveMesh::checkFaceAreas("
299 <<
"const bool, labelHashSet*) const: "
300 <<
"checking face area magnitudes" <<
endl;
305 scalar minArea = GREAT;
306 scalar maxArea = -GREAT;
308 forAll (magFaceAreas, faceI)
310 if (magFaceAreas[faceI] < VSMALL)
318 minArea =
min(minArea, magFaceAreas[faceI]);
319 maxArea =
max(maxArea, magFaceAreas[faceI]);
325 if (minArea < VSMALL)
329 Info<<
" ***Zero or negative face area detected. "
330 "Minimum area: " << minArea <<
endl;
339 Info<<
" Minumum face area = " << minArea
340 <<
". Maximum face area = " << maxArea
341 <<
". Face area magnitudes OK." <<
endl;
357 Info<<
"bool primitiveMesh::checkCellVolumes("
358 <<
"const bool, labelHashSet*) const: "
359 <<
"checking cell volumes" <<
endl;
364 scalar minVolume = GREAT;
365 scalar maxVolume = -GREAT;
367 label nNegVolCells = 0;
371 if (vols[cellI] < VSMALL)
381 minVolume =
min(minVolume, vols[cellI]);
382 maxVolume =
max(maxVolume, vols[cellI]);
389 if (minVolume < VSMALL)
393 Info<<
" ***Zero or negative cell volume detected. "
394 <<
"Minimum negative volume: " << minVolume
395 <<
", Number of negative volume cells: " << nNegVolCells
405 Info<<
" Min volume = " << minVolume
406 <<
". Max volume = " << maxVolume
407 <<
". Total volume = " <<
gSum(vols)
408 <<
". Cell volumes OK." <<
endl;
424 Info<<
"bool primitiveMesh::checkFaceOrthogonality("
425 <<
"const bool, labelHashSet*) const: "
426 <<
"checking mesh non-orthogonality" <<
endl;
437 const scalar severeNonorthogonalityThreshold =
440 scalar minDDotS = GREAT;
444 label severeNonOrth = 0;
446 label errorNonOrth = 0;
450 vector d = centres[nei[faceI]] - centres[own[faceI]];
451 const vector& s = areas[faceI];
453 scalar dDotS = (d & s)/(
mag(d)*
mag(s) + VSMALL);
455 if (dDotS < severeNonorthogonalityThreshold)
477 if (dDotS < minDDotS)
492 label neiSize = nei.
size();
499 Info<<
" Mesh non-orthogonality Max: "
507 if (severeNonOrth > 0)
509 Info<<
" *Number of severely non-orthogonal faces: "
510 << severeNonOrth <<
"." <<
endl;
514 if (errorNonOrth > 0)
518 Info<<
" ***Number of non-orthogonality errors: "
519 << errorNonOrth <<
"." <<
endl;
528 Info<<
" Non-orthogonality check OK." <<
endl;
539 const scalar minPyrVol,
545 Info<<
"bool primitiveMesh::checkFacePyramids("
546 <<
"const bool, const scalar, labelHashSet*) const: "
547 <<
"checking face orientation" <<
endl;
560 label nErrorPyrs = 0;
567 if (pyrVol > -minPyrVol)
577 if (isInternalFace(faceI))
583 if (pyrVol < minPyrVol)
601 Info<<
" ***Error in face pyramids: "
602 << nErrorPyrs <<
" faces are incorrectly oriented."
612 Info<<
" Face pyramids OK." <<
endl;
628 Info<<
"bool primitiveMesh::checkFaceSkewnesss("
629 <<
"const bool, labelHashSet*) const: "
630 <<
"checking face skewness" <<
endl;
650 vector Cpf = faceCtrs[faceI] - cellCtrs[own[faceI]];
651 vector d = cellCtrs[nei[faceI]] - cellCtrs[own[faceI]];
655 Cpf - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + SMALL))*
d;
661 scalar fd = 0.2*
mag(d) + VSMALL;
662 const face&
f = fcs[faceI];
665 fd =
max(fd,
mag(svHat & (p[f[
pi]] - faceCtrs[faceI])));
669 scalar skewness =
mag(sv)/fd;
673 if (skewness > skewThreshold_)
683 if(skewness > maxSkew)
693 for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
695 vector Cpf = faceCtrs[faceI] - cellCtrs[own[faceI]];
698 normal /=
mag(normal) + VSMALL;
699 vector d = normal*(normal & Cpf);
703 vector sv = Cpf - ((fAreas[faceI]&Cpf)/((fAreas[faceI]&d)+VSMALL))*
d;
709 scalar fd = 0.4*
mag(d) + VSMALL;
710 const face&
f = fcs[faceI];
713 fd =
max(fd,
mag(svHat & (p[f[
pi]] - faceCtrs[faceI])));
717 scalar skewness =
mag(sv)/fd;
721 if (skewness > skewThreshold_)
731 if(skewness > maxSkew)
745 Info<<
" ***Max skewness = " << maxSkew
746 <<
", " << nWarnSkew <<
" highly skew faces detected"
747 " which may impair the quality of the results"
757 Info<<
" Max skewness = " << maxSkew <<
" OK." <<
endl;
773 Info<<
"bool primitiveMesh::checkPoints"
774 <<
"(const bool, labelHashSet*) const: "
775 <<
"checking points" <<
endl;
778 label nFaceErrors = 0;
779 label nCellErrors = 0;
785 if (pf[pointI].empty())
815 if (nFaceErrors > 0 || nCellErrors > 0)
819 Info<<
" ***Unused points found in the mesh, "
820 "number unused by faces: " << nFaceErrors
821 <<
" number unused by cells: " << nCellErrors
852 Info<<
"bool primitiveMesh::checkFaceAngles"
853 <<
"(const bool, const scalar, labelHashSet*) const: "
854 <<
"checking face angles" <<
endl;
857 if (maxDeg < -SMALL || maxDeg > 180+SMALL)
861 "primitiveMesh::checkFaceAngles"
862 "(const bool, const scalar, labelHashSet*)"
863 ) <<
"maxDeg should be [0..180] but is now " << maxDeg
872 faceNormals /=
mag(faceNormals) + VSMALL;
874 scalar maxEdgeSin = 0.0;
878 label errorFaceI = -1;
882 const face&
f = fcs[faceI];
886 scalar magEPrev =
mag(ePrev);
887 ePrev /= magEPrev + VSMALL;
895 vector e10(p[f[fp1]] - p[f[fp0]]);
896 scalar magE10 =
mag(e10);
897 e10 /= magE10 + VSMALL;
899 if (magEPrev > SMALL && magE10 > SMALL)
901 vector edgeNormal = ePrev ^ e10;
902 scalar magEdgeNormal =
mag(edgeNormal);
904 if (magEdgeNormal < maxSin)
911 edgeNormal /= magEdgeNormal;
913 if ((edgeNormal & faceNormals[faceI]) < SMALL)
915 if (faceI != errorFaceI)
927 maxEdgeSin =
max(maxEdgeSin, magEdgeNormal);
942 scalar maxConcaveDegr =
948 Info<<
" *There are " << nConcave
949 <<
" faces with concave angles between consecutive"
950 <<
" edges. Max concave angle = " << maxConcaveDegr
951 <<
" degrees." <<
endl;
960 Info<<
" All angles in faces OK." <<
endl;
974 const scalar warnFlatness,
980 Info<<
"bool primitiveMesh::checkFaceFlatness"
981 <<
"(const bool, const scalar, labelHashSet*) const: "
982 <<
"checking face flatness" <<
endl;
985 if (warnFlatness < 0 || warnFlatness > 1)
989 "primitiveMesh::checkFaceFlatness"
990 "(const bool, const scalar, labelHashSet*)"
991 ) <<
"warnFlatness should be [0..1] but is now " << warnFlatness
1006 scalar minFlatness = GREAT;
1007 scalar sumFlatness = 0;
1012 const face&
f = fcs[faceI];
1014 if (f.
size() > 3 && magAreas[faceI] > VSMALL)
1016 const point& fc = fctrs[faceI];
1025 const point& thisPoint = p[f[fp]];
1026 const point& nextPoint = p[f.nextLabel(fp)];
1029 vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
1033 scalar flatness = magAreas[faceI] / (sumA+VSMALL);
1035 sumFlatness += flatness;
1038 minFlatness =
min(minFlatness, flatness);
1040 if (flatness < warnFlatness)
1059 if (debug || report)
1063 Info<<
" Face flatness (1 = flat, 0 = butterfly) : average = "
1064 << sumFlatness / nSummed <<
" min = " << minFlatness <<
endl;
1071 if (debug || report)
1073 Info<<
" *There are " << nWarped
1074 <<
" faces with ratio between projected and actual area < "
1075 << warnFlatness <<
endl;
1077 Info<<
" Minimum ratio (minimum flatness, maximum warpage) = "
1078 << minFlatness <<
endl;
1085 if (debug || report)
1087 Info<<
" All face flatness OK." <<
endl;
1110 Info<<
"bool primitiveMesh::checkEdgeAlignment("
1111 <<
"const bool, const Vector<label>&, labelHashSet*) const: "
1112 <<
"checking edge alignment" <<
endl;
1118 if (directions[cmpt] == 1)
1122 else if (directions[cmpt] != 0)
1126 "primitiveMesh::checkEdgeAlignment"
1127 "(const bool, const Vector<label>&, labelHashSet*)"
1128 ) <<
"directions should contain 0 or 1 but is now " << directions
1133 if (nDirs == vector::nComponents)
1147 const face&
f = fcs[faceI];
1156 scalar magD =
mag(d);
1158 if (magD > ROOTVSMALL)
1163 label nEmptyDirs = 0;
1164 label nNonEmptyDirs = 0;
1167 if (
mag(d[cmpt]) > 1
e-6)
1169 if (directions[cmpt] == 0)
1180 if (nEmptyDirs == 0)
1184 else if (nEmptyDirs == 1)
1187 if (nNonEmptyDirs > 0)
1192 else if (nEmptyDirs > 1)
1204 if (nErrorEdges > 0)
1206 if (debug || report)
1208 Info<<
" ***Number of edges not aligned with or perpendicular to "
1209 <<
"non-empty directions: " << nErrorEdges <<
endl;
1217 setPtr->
insert(iter.key()[0]);
1218 setPtr->
insert(iter.key()[1]);
1226 if (debug || report)
1228 Info<<
" All edges aligned with or perpendicular to "
1229 <<
"non-empty directions." <<
endl;
1244 Info<<
"bool primitiveMesh::checkUpperTriangular("
1245 <<
"const bool, labelHashSet*) const: "
1246 <<
"checking face ordering" <<
endl;
1255 label
internal = nInternalFaces();
1260 label nMultipleCells =
false;
1264 for (label faceI = 0; faceI <
internal; faceI++)
1266 if (own[faceI] >= nei[faceI])
1290 label faceI = curFaces[i];
1292 if (faceI >= nInternalFaces())
1299 label nbrCellI = nei[faceI];
1301 if (nbrCellI == cellI)
1303 nbrCellI = own[faceI];
1306 if (cellI < nbrCellI)
1326 label prevCell = nbr[0];
1327 label prevFace = curFaces[nbr.indices()[0]];
1329 bool hasMultipleFaces =
false;
1331 for (label i = 1; i < nbr.size(); i++)
1333 label thisCell = nbr[i];
1334 label thisFace = curFaces[nbr.indices()[i]];
1336 if (thisCell == labelMax)
1341 if (thisCell == prevCell)
1343 hasMultipleFaces =
true;
1347 setPtr->
insert(prevFace);
1348 setPtr->
insert(thisFace);
1351 else if (thisFace < prevFace)
1357 setPtr->
insert(thisFace);
1361 prevCell = thisCell;
1362 prevFace = thisFace;
1365 if (hasMultipleFaces)
1374 if ((debug || report) && nMultipleCells > 0)
1376 Info<<
" <<Found " << nMultipleCells
1377 <<
" neighbouring cells with multiple inbetween faces." <<
endl;
1382 if (debug || report)
1384 Info<<
" ***Faces not in upper triangular order." <<
endl;
1391 if (debug || report)
1393 Info<<
" Upper triangular ordering OK." <<
endl;
1409 Info<<
"bool primitiveMesh::checkCellsZipUp("
1410 <<
"const bool, labelHashSet*) const: "
1411 <<
"checking topological cell openness" <<
endl;
1414 label nOpenCells = 0;
1423 const edgeList cellEdges = c[cellI].edges(f);
1429 edgeList curFaceEdges = f[curFaces[faceI]].edges();
1431 forAll (curFaceEdges, faceEdgeI)
1433 const edge& curEdge = curFaceEdges[faceEdgeI];
1435 forAll (cellEdges, cellEdgeI)
1437 if (cellEdges[cellEdgeI] == curEdge)
1439 edgeUsage[cellEdgeI]++;
1447 label nSingleEdges = 0;
1449 forAll (edgeUsage, edgeI)
1451 if (edgeUsage[edgeI] == 1)
1453 singleEdges[nSingleEdges] = cellEdges[edgeI];
1456 else if (edgeUsage[edgeI] != 2)
1465 if (nSingleEdges > 0)
1480 if (debug || report)
1482 Info<<
" ***Open cells found, number of cells: " << nOpenCells
1483 <<
". This problem may be fixable using the zipUpMesh utility."
1491 if (debug || report)
1493 Info<<
" Topological cell zip-up check OK." <<
endl;
1510 Info<<
"bool primitiveMesh::checkFaceVertices("
1511 <<
"const bool, labelHashSet*) const: "
1512 <<
"checking face vertices" <<
endl;
1518 label nErrorFaces = 0;
1522 const face& curFace = f[fI];
1539 bool inserted = facePoints.insert(curFace[fp]);
1555 if (nErrorFaces > 0)
1557 if (debug || report)
1559 Info<<
" Faces with invalid vertex labels found, "
1560 <<
" number of faces: " << nErrorFaces <<
endl;
1567 if (debug || report)
1569 Info<<
" Face vertices OK." <<
endl;
1578 bool Foam::primitiveMesh::checkDuplicateFaces
1582 label& nBaffleFaces,
1590 label nbFaceI = iter.key();
1591 label nCommon = iter();
1593 const face& curFace = faces()[faceI];
1594 const face& nbFace = faces()[nbFaceI];
1596 if (nCommon == nbFace.
size() || nCommon == curFace.
size())
1598 if (nbFace.
size() != curFace.
size())
1620 bool Foam::primitiveMesh::checkCommonOrder
1623 const Map<label>& nCommonPoints,
1631 label nbFaceI = iter.key();
1632 label nCommon = iter();
1634 const face& curFace = faces()[faceI];
1635 const face& nbFace = faces()[nbFaceI];
1640 && nCommon != nbFace.size()
1641 && nCommon != curFace.size()
1647 label nb =
findIndex(nbFace, curFace[fp]);
1661 label fpPlus1 = curFace.fcIndex(fp);
1662 label fpMin1 = curFace.rcIndex(fp);
1665 label nbPlus1 = nbFace.fcIndex(nb);
1666 label nbMin1 = nbFace.rcIndex(nb);
1670 label curInc = labelMax;
1671 label nbInc = labelMax;
1673 if (nbFace[nbPlus1] == curFace[fpPlus1])
1678 else if (nbFace[nbPlus1] == curFace[fpMin1])
1683 else if (nbFace[nbMin1] == curFace[fpMin1])
1703 if (curFp >= curFace.size())
1709 curFp = curFace.size()-1;
1714 if (curNb >= nbFace.size())
1720 curNb = nbFace.size()-1;
1722 }
while (curFace[curFp] == nbFace[curNb]);
1731 for (label commonI = 0; commonI < nCommon; commonI++)
1735 if (curFp >= curFace.size())
1741 curFp = curFace.size()-1;
1746 if (curNb >= nbFace.size())
1752 curNb = nbFace.size()-1;
1755 if (curFace[curFp] != nbFace[curNb])
1759 setPtr->insert(faceI);
1760 setPtr->insert(nbFaceI);
1791 Info<<
"bool primitiveMesh::checkFaceFaces(const bool, labelHashSet*)"
1792 <<
" const: " <<
"checking face-face connectivity" <<
endl;
1797 label nBaffleFaces = 0;
1798 label nErrorDuplicate = 0;
1799 label nErrorOrder = 0;
1802 for (label faceI = 0; faceI < nFaces(); faceI++)
1804 const face& curFace = faces()[faceI];
1808 nCommonPoints.
clear();
1812 label pointI = curFace[fp];
1818 label nbFaceI = nbs[nbI];
1820 if (faceI < nbFaceI)
1826 if (fnd == nCommonPoints.
end())
1829 nCommonPoints.
insert(nbFaceI, 1);
1842 if (checkDuplicateFaces(faceI, nCommonPoints, nBaffleFaces, setPtr))
1848 if (checkCommonOrder(faceI, nCommonPoints, setPtr))
1860 Info<<
" Number of identical duplicate faces (baffle faces): "
1861 << nBaffleFaces <<
endl;
1864 if (nErrorDuplicate > 0 || nErrorOrder > 0)
1866 if (nErrorDuplicate > 0)
1868 Info<<
" ***Number of duplicate (not baffle) faces found: "
1869 << nErrorDuplicate <<
endl;
1872 if (nErrorOrder > 0)
1874 Info<<
" ***Number of faces with non-consecutive shared points: "
1875 << nErrorOrder <<
endl;
1882 if (debug || report)
1884 Info<<
" Face-face connectivity OK." <<
endl;
1902 Info<<
"bool primitiveMesh::checkCellDeterminant(const bool"
1903 <<
", labelHashSet*) const: "
1904 <<
"checking for under-determined cells" <<
endl;
1912 if (meshD[dir] == 1)
1925 label nErrorCells = 0;
1927 scalar minDet = GREAT;
1934 sumDet = c.
size()*minDet;
1946 label nInternalFaces = 0;
1950 if (isInternalFace(curFaces[i]))
1952 avgArea +=
mag(faceAreas()[curFaces[i]]);
1958 if (nInternalFaces == 0)
1969 avgArea /= nInternalFaces;
1975 if (isInternalFace(curFaces[i]))
1977 areaTensor +=
sqr(faceAreas()[curFaces[i]]/avgArea);
1987 areaTensor.
xx() = 1;
1991 areaTensor.
yy() = 1;
1995 areaTensor.
zz() = 1;
1999 scalar determinant =
mag(
det(areaTensor));
2001 minDet =
min(determinant, minDet);
2002 sumDet += determinant;
2005 if (determinant < 1
e-3)
2023 if (debug || report)
2027 Info<<
" Cell determinant (wellposedness) : minimum: " << minDet
2028 <<
" average: " << sumDet/nSummed
2033 if (nErrorCells > 0)
2035 if (debug || report)
2037 Info<<
" ***Cells with small determinant found, number of cells: "
2038 << nErrorCells <<
endl;
2045 if (debug || report)
2047 Info<<
" Cell determinant check OK." <<
endl;
2061 label noFailedChecks = 0;
2063 if (checkPoints(report)) noFailedChecks++;
2064 if (checkUpperTriangular(report)) noFailedChecks++;
2065 if (checkCellsZipUp(report)) noFailedChecks++;
2066 if (checkFaceVertices(report)) noFailedChecks++;
2067 if (checkFaceFaces(report)) noFailedChecks++;
2070 if (noFailedChecks == 0)
2072 if (debug || report)
2074 Info<<
" Mesh topology OK." <<
endl;
2081 if (debug || report)
2083 Info<<
" Failed " << noFailedChecks
2084 <<
" mesh topology checks." <<
endl;
2094 label noFailedChecks = 0;
2096 if (checkClosedBoundary(report)) noFailedChecks++;
2097 if (checkClosedCells(report)) noFailedChecks++;
2098 if (checkFaceAreas(report)) noFailedChecks++;
2099 if (checkCellVolumes(report)) noFailedChecks++;
2100 if (checkFaceOrthogonality(report)) noFailedChecks++;
2101 if (checkFacePyramids(report)) noFailedChecks++;
2102 if (checkFaceSkewness(report)) noFailedChecks++;
2104 if (noFailedChecks == 0)
2106 if (debug || report)
2108 Info<<
" Mesh geometry OK." <<
endl;
2114 if (debug || report)
2116 Info<<
" Failed " << noFailedChecks
2117 <<
" mesh geometry checks." <<
endl;
2129 Info<<
"bool primitiveMesh::checkMesh(const bool report) const: "
2130 <<
"checking primitiveMesh" <<
endl;
2135 if (noFailedChecks == 0)
2137 if (debug || report)
2146 if (debug || report)
2148 Info<<
" Failed " << noFailedChecks
2149 <<
" mesh checks." <<
endl;
2159 scalar prev = closedThreshold_;
2160 closedThreshold_ = val;
2168 scalar prev = aspectThreshold_;
2169 aspectThreshold_ = val;
2177 scalar prev = nonOrthThreshold_;
2178 nonOrthThreshold_ = val;
2186 scalar prev = skewThreshold_;
2187 skewThreshold_ = val;