42 void Foam::polyMeshGeometry::updateFaceCentresAndAreas
52 label facei = changedFaces[i];
61 faceCentres_[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
62 faceAreas_[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
70 point fCentre = p[f[0]];
80 const point& nextPoint = p[f[(
pi + 1) % nPoints]];
82 vector c = p[f[
pi]] + nextPoint + fCentre;
83 vector n = (nextPoint - p[f[
pi]])^(fCentre - p[f[
pi]]);
91 faceCentres_[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL);
92 faceAreas_[facei] = 0.5*sumN;
98 void Foam::polyMeshGeometry::updateCellCentresAndVols
105 UIndirectList<vector>(cellCentres_, changedCells) =
vector::zero;
106 UIndirectList<scalar>(cellVolumes_, changedCells) = 0.0;
108 const labelList& own = mesh_.faceOwner();
109 const labelList& nei = mesh_.faceNeighbour();
114 UIndirectList<vector>(cEst, changedCells) =
vector::zero;
116 UIndirectList<scalar>(nCellFaces, changedCells) = 0.0;
120 label faceI = changedFaces[i];
121 cEst[own[faceI]] += faceCentres_[faceI];
122 nCellFaces[own[faceI]] += 1;
124 if (mesh_.isInternalFace(faceI))
126 cEst[nei[faceI]] += faceCentres_[faceI];
127 nCellFaces[nei[faceI]] += 1;
133 label cellI = changedCells[i];
134 cEst[cellI] /= nCellFaces[cellI];
139 label faceI = changedFaces[i];
144 faceAreas_[faceI] & (faceCentres_[faceI] - cEst[own[faceI]]),
149 vector pc = (3.0/4.0)*faceCentres_[faceI] + (1.0/4.0)*cEst[own[faceI]];
152 cellCentres_[own[faceI]] += pyr3Vol*
pc;
155 cellVolumes_[own[faceI]] += pyr3Vol;
157 if (mesh_.isInternalFace(faceI))
162 faceAreas_[faceI] & (cEst[nei[faceI]] - faceCentres_[faceI]),
168 (3.0/4.0)*faceCentres_[faceI]
169 + (1.0/4.0)*cEst[nei[faceI]];
172 cellCentres_[nei[faceI]] += pyr3Vol*
pc;
175 cellVolumes_[nei[faceI]] += pyr3Vol;
181 label cellI = changedCells[i];
183 cellCentres_[cellI] /= cellVolumes_[cellI] + VSMALL;
184 cellVolumes_[cellI] *= (1.0/3.0);
202 label faceI = changedFaces[i];
204 affectedCells.insert(own[faceI]);
208 affectedCells.insert(nei[faceI]);
211 return affectedCells.toc();
215 Foam::scalar Foam::polyMeshGeometry::checkNonOrtho
219 const scalar severeNonorthogonalityThreshold,
224 label& severeNonOrth,
229 scalar dDotS = (d & s)/(
mag(d)*
mag(s) + VSMALL);
231 if (dDotS < severeNonorthogonalityThreshold)
245 Pout<<
"Severe non-orthogonality for face " << faceI
246 <<
" between cells " << mesh.
faceOwner()[faceI]
262 "polyMeshGeometry::checkFaceDotProduct"
263 "(const bool, const scalar, const labelList&"
265 ) <<
"Severe non-orthogonality detected for face "
267 <<
" between cells " << mesh.
faceOwner()[faceI]
286 Foam::scalar Foam::polyMeshGeometry::calcSkewness
293 scalar dOwn =
mag(fc - ownCc);
294 scalar dNei =
mag(fc - neiCc);
296 point faceIntersection =
297 ownCc*dNei/(dOwn+dNei+VSMALL)
298 + neiCc*dOwn/(dOwn+dNei+VSMALL);
301 mag(fc - faceIntersection)
328 faceAreas_ = mesh_.faceAreas();
329 faceCentres_ = mesh_.faceCentres();
330 cellCentres_ = mesh_.cellCentres();
331 cellVolumes_ = mesh_.cellVolumes();
343 updateFaceCentresAndAreas(p, changedFaces);
345 updateCellCentresAndVols(affectedCells(mesh_, changedFaces), changedFaces);
352 const scalar orthWarn,
369 const scalar severeNonorthogonalityThreshold =
383 scalar minDDotS = GREAT;
388 label severeNonOrth = 0;
390 label errorNonOrth = 0;
394 label faceI = checkFaces[i];
396 const point& ownCc = cellCentres[own[faceI]];
400 scalar dDotS = checkNonOrtho
404 severeNonorthogonalityThreshold,
407 cellCentres[nei[faceI]] - ownCc,
414 if (dDotS < minDDotS)
426 if (patches[patchI].coupled())
428 scalar dDotS = checkNonOrtho
432 severeNonorthogonalityThreshold,
442 if (dDotS < minDDotS)
455 label face0 = baffles[i].first();
456 label face1 = baffles[i].second();
458 const point& ownCc = cellCentres[own[face0]];
460 scalar dDotS = checkNonOrtho
464 severeNonorthogonalityThreshold,
467 cellCentres[own[face1]] - ownCc,
474 if (dDotS < minDDotS)
492 if (report && minDDotS < severeNonorthogonalityThreshold)
494 Info<<
"Number of non-orthogonality errors: " << errorNonOrth
495 <<
". Number of severely non-orthogonal faces: "
496 << severeNonOrth <<
"." <<
endl;
504 Info<<
"Mesh non-orthogonality Max: "
512 if (errorNonOrth > 0)
518 "polyMeshGeometry::checkFaceDotProduct"
519 "(const bool, const scalar, const labelList&, labelHashSet*)"
520 ) <<
"Error in non-orthogonality detected" <<
endl;
529 Info<<
"Non-orthogonality check OK.\n" <<
endl;
540 const scalar minPyrVol,
555 label nErrorPyrs = 0;
559 label faceI = checkFaces[i];
565 cellCentres[own[faceI]]
568 if (pyrVol > -minPyrVol)
572 Pout<<
"bool polyMeshGeometry::checkFacePyramids("
573 <<
"const bool, const scalar, const pointField&"
574 <<
", const labelList&, labelHashSet*): "
575 <<
"face " << faceI <<
" points the wrong way. " <<
endl
576 <<
"Pyramid volume: " << -pyrVol
577 <<
" Face " << f[faceI] <<
" area: " << f[faceI].mag(p)
578 <<
" Owner cell: " << own[faceI] <<
endl
579 <<
"Owner cell vertex labels: "
580 << mesh.
cells()[own[faceI]].labels(f)
599 if (pyrVol < minPyrVol)
603 Pout<<
"bool polyMeshGeometry::checkFacePyramids("
604 <<
"const bool, const scalar, const pointField&"
605 <<
", const labelList&, labelHashSet*): "
606 <<
"face " << faceI <<
" points the wrong way. " <<
endl
607 <<
"Pyramid volume: " << -pyrVol
608 <<
" Face " << f[faceI] <<
" area: " << f[faceI].mag(p)
609 <<
" Neighbour cell: " << nei[faceI] <<
endl
610 <<
"Neighbour cell vertex labels: "
611 << mesh.
cells()[nei[faceI]].labels(f)
627 label face0 = baffles[i].first();
628 label face1 = baffles[i].second();
630 const point& ownCc = cellCentres[own[face0]];
639 if (pyrVolOwn > -minPyrVol)
643 Pout<<
"bool polyMeshGeometry::checkFacePyramids("
644 <<
"const bool, const scalar, const pointField&"
645 <<
", const labelList&, labelHashSet*): "
646 <<
"face " << face0 <<
" points the wrong way. " <<
endl
647 <<
"Pyramid volume: " << -pyrVolOwn
648 <<
" Face " << f[face0] <<
" area: " << f[face0].mag(p)
649 <<
" Owner cell: " << own[face0] <<
endl
650 <<
"Owner cell vertex labels: "
651 << mesh.
cells()[own[face0]].labels(f)
668 if (pyrVolNbr < minPyrVol)
672 Pout<<
"bool polyMeshGeometry::checkFacePyramids("
673 <<
"const bool, const scalar, const pointField&"
674 <<
", const labelList&, labelHashSet*): "
675 <<
"face " << face0 <<
" points the wrong way. " <<
endl
676 <<
"Pyramid volume: " << -pyrVolNbr
677 <<
" Face " << f[face0] <<
" area: " << f[face0].mag(p)
678 <<
" Neighbour cell: " << own[face1] <<
endl
679 <<
"Neighbour cell vertex labels: "
680 << mesh.
cells()[own[face1]].labels(f)
701 "polyMeshGeometry::checkFacePyramids("
702 "const bool, const scalar, const pointField&"
703 ", const labelList&, labelHashSet*)"
704 ) <<
"Error in face pyramids: faces pointing the wrong way!"
714 Info<<
"Face pyramids OK.\n" <<
endl;
725 const scalar internalSkew,
726 const scalar boundarySkew,
759 label faceI = checkFaces[i];
763 scalar skewness = calcSkewness
765 cellCentres[own[faceI]],
766 cellCentres[nei[faceI]],
773 if (skewness > internalSkew)
777 Pout<<
"Severe skewness for face " << faceI
778 <<
" skewness = " << skewness <<
endl;
789 maxSkew =
max(maxSkew, skewness);
791 else if (patches[patches.
whichPatch(faceI)].coupled())
793 scalar skewness = calcSkewness
795 cellCentres[own[faceI]],
803 if (skewness > internalSkew)
807 Pout<<
"Severe skewness for coupled face " << faceI
808 <<
" skewness = " << skewness <<
endl;
819 maxSkew =
max(maxSkew, skewness);
826 vector faceNormal = faceAreas[faceI];
827 faceNormal /=
mag(faceNormal) + ROOTVSMALL;
829 vector dOwn = faceCentres[faceI] - cellCentres[own[faceI]];
831 vector dWall = faceNormal*(faceNormal & dOwn);
833 point faceIntersection = cellCentres[own[faceI]] + dWall;
836 mag(faceCentres[faceI] - faceIntersection)
837 /(2*
mag(dWall) + ROOTVSMALL);
842 if (skewness > boundarySkew)
846 Pout<<
"Severe skewness for boundary face " << faceI
847 <<
" skewness = " << skewness <<
endl;
858 maxSkew =
max(maxSkew, skewness);
864 label face0 = baffles[i].first();
865 label face1 = baffles[i].second();
867 const point& ownCc = cellCentres[own[face0]];
869 scalar skewness = calcSkewness
872 cellCentres[own[face1]],
879 if (skewness > internalSkew)
883 Pout<<
"Severe skewness for face " << face0
884 <<
" skewness = " << skewness <<
endl;
895 maxSkew =
max(maxSkew, skewness);
908 "polyMeshGeometry::checkFaceSkewness"
909 "(const bool, const scalar, const labelList&, labelHashSet*)"
910 ) <<
"Large face skewness detected. Max skewness = "
912 <<
" percent.\nThis may impair the quality of the result." <<
nl
913 << nWarnSkew <<
" highly skew faces detected."
923 Info<<
"Max skewness = " << 100*maxSkew
924 <<
" percent. Face skewness OK.\n" <<
endl;
935 const scalar warnWeight,
961 scalar minWeight = GREAT;
963 label nWarnWeight = 0;
967 label faceI = checkFaces[i];
969 const point& fc = faceCentres[faceI];
970 const vector& fa = faceAreas[faceI];
972 scalar dOwn =
mag(fa & (fc-cellCentres[own[faceI]]));
976 scalar dNei =
mag(fa & (cellCentres[nei[faceI]]-fc));
977 scalar weight =
min(dNei,dOwn)/(dNei+dOwn+VSMALL);
979 if (weight < warnWeight)
983 Pout<<
"Small weighting factor for face " << faceI
984 <<
" weight = " << weight <<
endl;
995 minWeight =
min(minWeight, weight);
1001 if (patches[patchI].coupled())
1004 scalar weight =
min(dNei,dOwn)/(dNei+dOwn+VSMALL);
1006 if (weight < warnWeight)
1010 Pout<<
"Small weighting factor for face " << faceI
1011 <<
" weight = " << weight <<
endl;
1022 minWeight =
min(minWeight, weight);
1029 label face0 = baffles[i].first();
1030 label face1 = baffles[i].second();
1032 const point& ownCc = cellCentres[own[face0]];
1033 const point& fc = faceCentres[face0];
1034 const vector& fa = faceAreas[face0];
1036 scalar dOwn =
mag(fa & (fc-ownCc));
1037 scalar dNei =
mag(fa & (cellCentres[own[face1]]-fc));
1038 scalar weight =
min(dNei,dOwn)/(dNei+dOwn+VSMALL);
1040 if (weight < warnWeight)
1044 Pout<<
"Small weighting factor for face " << face0
1045 <<
" weight = " << weight <<
endl;
1056 minWeight =
min(minWeight, weight);
1062 if (minWeight < warnWeight)
1068 "polyMeshGeometry::checkFaceWeights"
1069 "(const bool, const scalar, const labelList&, labelHashSet*)"
1070 ) <<
"Small interpolation weight detected. Min weight = "
1071 << minWeight <<
'.' <<
nl
1072 << nWarnWeight <<
" faces with small weights detected."
1082 Info<<
"Min weight = " << minWeight
1083 <<
". Weights OK.\n" <<
endl;
1094 const scalar warnRatio,
1118 scalar minRatio = GREAT;
1120 label nWarnRatio = 0;
1124 label faceI = checkFaces[i];
1126 scalar ownVol =
mag(cellVolumes[own[faceI]]);
1128 scalar neiVol = -GREAT;
1132 neiVol =
mag(cellVolumes[nei[faceI]]);
1138 if (patches[patchI].coupled())
1146 scalar ratio =
min(ownVol, neiVol) / (
max(ownVol, neiVol) + VSMALL);
1148 if (ratio < warnRatio)
1152 Pout<<
"Small ratio for face " << faceI
1153 <<
" ratio = " << ratio <<
endl;
1164 minRatio =
min(minRatio, ratio);
1170 label face0 = baffles[i].first();
1171 label face1 = baffles[i].second();
1173 scalar ownVol =
mag(cellVolumes[own[face0]]);
1175 scalar neiVol =
mag(cellVolumes[own[face1]]);
1179 scalar ratio =
min(ownVol, neiVol) / (
max(ownVol, neiVol) + VSMALL);
1181 if (ratio < warnRatio)
1185 Pout<<
"Small ratio for face " << face0
1186 <<
" ratio = " << ratio <<
endl;
1197 minRatio =
min(minRatio, ratio);
1204 if (minRatio < warnRatio)
1210 "polyMeshGeometry::checkVolRatio"
1211 "(const bool, const scalar, const labelList&, labelHashSet*)"
1212 ) <<
"Small volume ratio detected. Min ratio = "
1213 << minRatio <<
'.' <<
nl
1214 << nWarnRatio <<
" faces with small ratios detected."
1224 Info<<
"Min ratio = " << minRatio
1225 <<
". Ratios OK.\n" <<
endl;
1240 const scalar maxDeg,
1248 if (maxDeg < -SMALL || maxDeg > 180+SMALL)
1252 "polyMeshGeometry::checkFaceAngles"
1253 "(const bool, const scalar, const pointField&, const labelList&"
1255 ) <<
"maxDeg should be [0..180] but is now " << maxDeg
1263 scalar maxEdgeSin = 0.0;
1267 label errorFaceI = -1;
1271 label faceI = checkFaces[i];
1273 const face& f = fcs[faceI];
1275 vector faceNormal = faceAreas[faceI];
1276 faceNormal /=
mag(faceNormal) + VSMALL;
1280 scalar magEPrev =
mag(ePrev);
1281 ePrev /= magEPrev + VSMALL;
1289 vector e10(p[f[fp1]] - p[f[fp0]]);
1290 scalar magE10 =
mag(e10);
1291 e10 /= magE10 + VSMALL;
1293 if (magEPrev > SMALL && magE10 > SMALL)
1295 vector edgeNormal = ePrev ^ e10;
1296 scalar magEdgeNormal =
mag(edgeNormal);
1298 if (magEdgeNormal < maxSin)
1305 edgeNormal /= magEdgeNormal;
1307 if ((edgeNormal & faceNormal) < SMALL)
1309 if (faceI != errorFaceI)
1321 maxEdgeSin =
max(maxEdgeSin, magEdgeNormal);
1336 if (maxEdgeSin > SMALL)
1338 scalar maxConcaveDegr =
1342 Info<<
"There are " << nConcave
1343 <<
" faces with concave angles between consecutive"
1344 <<
" edges. Max concave angle = " << maxConcaveDegr
1345 <<
" degrees.\n" <<
endl;
1349 Info<<
"All angles in faces are convex or less than " << maxDeg
1350 <<
" degrees concave.\n" <<
endl;
1360 "polyMeshGeometry::checkFaceAngles"
1361 "(const bool, const scalar, const pointField&"
1362 ", const labelList&, labelHashSet*)"
1363 ) << nConcave <<
" face points with severe concave angle (> "
1364 << maxDeg <<
" deg) found.\n"
1382 const scalar minTwist,
1392 if (minTwist < -1-SMALL || minTwist > 1+SMALL)
1396 "polyMeshGeometry::checkFaceTwist"
1397 "(const bool, const scalar, const polyMesh&, const pointField&"
1398 ", const pointField&, const pointField&, const pointField&"
1399 ", const labelList&, labelHashSet*)"
1400 ) <<
"minTwist should be [-1..1] but is now " << minTwist
1469 label faceI = checkFaces[i];
1471 const face& f = fcs[faceI];
1479 nf = cellCentres[nei[faceI]] - cellCentres[own[faceI]];
1480 nf /=
mag(nf) + VSMALL;
1482 else if (patches[patches.
whichPatch(faceI)].coupled())
1486 - cellCentres[own[faceI]];
1487 nf /=
mag(nf) + VSMALL;
1491 nf = faceCentres[faceI] - cellCentres[own[faceI]];
1492 nf /=
mag(nf) + VSMALL;
1497 const point& fc = faceCentres[faceI];
1511 scalar magTri =
mag(triArea);
1513 if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
1535 Info<<
"There are " << nWarped
1536 <<
" faces with cosine of the angle"
1537 <<
" between triangle normal and face normal less than "
1538 << minTwist <<
nl <<
endl;
1542 Info<<
"All faces are flat in that the cosine of the angle"
1543 <<
" between triangle normal and face normal less than "
1544 << minTwist <<
nl <<
endl;
1554 "polyMeshGeometry::checkFaceTwist"
1555 "(const bool, const scalar, const polyMesh&, const pointField&"
1556 ", const pointField&, const pointField&, const pointField&"
1557 ", const labelList&, labelHashSet*)"
1558 ) << nWarped <<
" faces with severe warpage "
1559 <<
"(cosine of the angle between triangle normal and face normal"
1560 <<
" < " << minTwist <<
") found.\n"
1577 const scalar minTwist,
1586 if (minTwist < -1-SMALL || minTwist > 1+SMALL)
1590 "polyMeshGeometry::checkTriangleTwist"
1591 "(const bool, const scalar, const polyMesh&, const pointField&"
1592 ", const labelList&, labelHashSet*)"
1593 ) <<
"minTwist should be [-1..1] but is now " << minTwist
1603 label faceI = checkFaces[i];
1605 const face& f = fcs[faceI];
1609 const point& fc = faceCentres[faceI];
1624 scalar magTri =
mag(prevN);
1626 if (magTri > VSMALL)
1651 scalar magTri =
mag(triN);
1653 if (magTri > VSMALL)
1657 if ((prevN & triN) < minTwist)
1671 else if (minTwist > 0)
1683 while (fp != startFp);
1695 Info<<
"There are " << nWarped
1696 <<
" faces with cosine of the angle"
1697 <<
" between consecutive triangle normals less than "
1698 << minTwist <<
nl <<
endl;
1702 Info<<
"All faces are flat in that the cosine of the angle"
1703 <<
" between consecutive triangle normals is less than "
1704 << minTwist <<
nl <<
endl;
1714 "polyMeshGeometry::checkTriangleTwist"
1715 "(const bool, const scalar, const polyMesh&"
1716 ", const pointField&, const labelList&, labelHashSet*)"
1717 ) << nWarped <<
" faces with severe warpage "
1718 <<
"(cosine of the angle between consecutive triangle normals"
1719 <<
" < " << minTwist <<
") found.\n"
1735 const scalar minArea,
1742 label nZeroArea = 0;
1746 label faceI = checkFaces[i];
1748 if (
mag(faceAreas[faceI]) < minArea)
1765 Info<<
"There are " << nZeroArea
1766 <<
" faces with area < " << minArea <<
'.' <<
nl <<
endl;
1770 Info<<
"All faces have area > " << minArea <<
'.' <<
nl <<
endl;
1780 "polyMeshGeometry::checkFaceArea"
1781 "(const bool, const scalar, const polyMesh&"
1782 ", const pointField&, const labelList&, labelHashSet*)"
1783 ) << nZeroArea <<
" faces with area < " << minArea
1800 const scalar warnDet,
1810 scalar minDet = GREAT;
1811 scalar sumDet = 0.0;
1817 const cell& cFaces = cells[affectedCells[i]];
1820 scalar magAreaSum = 0;
1824 label faceI = cFaces[cFaceI];
1826 scalar magArea =
mag(faceAreas[faceI]);
1828 magAreaSum += magArea;
1829 areaSum += faceAreas[faceI]*(faceAreas[faceI]/(magArea+VSMALL));
1832 scalar scaledDet =
det(areaSum/(magAreaSum+VSMALL))/0.037037037037037;
1834 minDet =
min(minDet, scaledDet);
1835 sumDet += scaledDet;
1838 if (scaledDet < warnDet)
1845 label faceI = cFaces[cFaceI];
1862 Info<<
"Cell determinant (1 = uniform cube) : average = "
1863 << sumDet / nSumDet <<
" min = " << minDet <<
endl;
1868 Info<<
"There are " << nWarnDet
1869 <<
" cells with determinant < " << warnDet <<
'.' <<
nl
1874 Info<<
"All faces have determinant > " << warnDet <<
'.' <<
nl
1885 "polyMeshGeometry::checkCellDeterminant"
1886 "(const bool, const scalar, const polyMesh&"
1887 ", const pointField&, const labelList&, const labelList&"
1889 ) << nWarnDet <<
" cells with determinant < " << warnDet
1906 const scalar orthWarn,
1912 return checkFaceDotProduct
1929 const scalar minPyrVol,
1936 return checkFacePyramids
1953 const scalar internalSkew,
1954 const scalar boundarySkew,
1960 return checkFaceSkewness
1979 const scalar warnWeight,
1985 return checkFaceWeights
2003 const scalar warnRatio,
2009 return checkVolRatio
2025 const scalar maxDeg,
2031 return checkFaceAngles
2047 const scalar minTwist,
2053 return checkFaceTwist
2071 const scalar minTwist,
2077 return checkTriangleTwist
2094 const scalar minArea,
2099 return checkFaceArea
2114 const scalar warnDet,
2120 return checkCellDeterminant