41 void Foam::primitiveMeshGeometry::updateFaceCentresAndAreas
51 label facei = changedFaces[i];
60 faceCentres_[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
61 faceAreas_[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
69 point fCentre = p[f[0]];
79 const point& nextPoint = p[f[(
pi + 1) % nPoints]];
81 vector c = p[f[
pi]] + nextPoint + fCentre;
82 vector n = (nextPoint - p[f[
pi]])^(fCentre - p[f[
pi]]);
90 faceCentres_[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL);
91 faceAreas_[facei] = 0.5*sumN;
97 void Foam::primitiveMeshGeometry::updateCellCentresAndVols
104 UIndirectList<vector>(cellCentres_, changedCells) =
vector::zero;
105 UIndirectList<scalar>(cellVolumes_, changedCells) = 0.0;
107 const labelList& own = mesh_.faceOwner();
108 const labelList& nei = mesh_.faceNeighbour();
113 UIndirectList<vector>(cEst, changedCells) =
vector::zero;
115 UIndirectList<scalar>(nCellFaces, changedCells) = 0.0;
119 label faceI = changedFaces[i];
120 cEst[own[faceI]] += faceCentres_[faceI];
121 nCellFaces[own[faceI]] += 1;
123 if (mesh_.isInternalFace(faceI))
125 cEst[nei[faceI]] += faceCentres_[faceI];
126 nCellFaces[nei[faceI]] += 1;
132 label cellI = changedCells[i];
133 cEst[cellI] /= nCellFaces[cellI];
138 label faceI = changedFaces[i];
143 faceAreas_[faceI] & (faceCentres_[faceI] - cEst[own[faceI]]),
148 vector pc = (3.0/4.0)*faceCentres_[faceI] + (1.0/4.0)*cEst[own[faceI]];
151 cellCentres_[own[faceI]] += pyr3Vol*
pc;
154 cellVolumes_[own[faceI]] += pyr3Vol;
156 if (mesh_.isInternalFace(faceI))
161 faceAreas_[faceI] & (cEst[nei[faceI]] - faceCentres_[faceI]),
167 (3.0/4.0)*faceCentres_[faceI]
168 + (1.0/4.0)*cEst[nei[faceI]];
171 cellCentres_[nei[faceI]] += pyr3Vol*
pc;
174 cellVolumes_[nei[faceI]] += pyr3Vol;
180 label cellI = changedCells[i];
182 cellCentres_[cellI] /= cellVolumes_[cellI];
183 cellVolumes_[cellI] *= (1.0/3.0);
193 const labelList& own = mesh_.faceOwner();
194 const labelList& nei = mesh_.faceNeighbour();
200 label faceI = changedFaces[i];
202 affectedCells.insert(own[faceI]);
204 if (mesh_.isInternalFace(faceI))
206 affectedCells.insert(nei[faceI]);
209 return affectedCells.toc();
236 faceAreas_ = mesh_.faceAreas();
237 faceCentres_ = mesh_.faceCentres();
238 cellCentres_ = mesh_.cellCentres();
239 cellVolumes_ = mesh_.cellVolumes();
251 updateFaceCentresAndAreas(p, changedFaces);
253 updateCellCentresAndVols(affectedCells(changedFaces), changedFaces);
260 const scalar orthWarn,
274 const scalar severeNonorthogonalityThreshold =
277 scalar minDDotS = GREAT;
281 label severeNonOrth = 0;
283 label errorNonOrth = 0;
287 label faceI = checkFaces[i];
291 vector d = cellCentres[nei[faceI]] - cellCentres[own[faceI]];
292 const vector& s = faceAreas[faceI];
294 scalar dDotS = (d & s)/(
mag(d)*
mag(s) + VSMALL);
296 if (dDotS < severeNonorthogonalityThreshold)
303 Pout<<
"Severe non-orthogonality for face " << faceI
304 <<
" between cells " << own[faceI]
305 <<
" and " << nei[faceI]
325 "primitiveMeshGeometry::checkFaceDotProduct"
326 "(const bool, const scalar, const labelList&"
328 ) <<
"Severe non-orthogonality detected for face "
330 <<
" between cells " << own[faceI] <<
" and "
346 if (dDotS < minDDotS)
360 label neiSize = nei.
size();
366 if (report && minDDotS < severeNonorthogonalityThreshold)
368 Info<<
"Number of non-orthogonality errors: " << errorNonOrth
369 <<
". Number of severely non-orthogonal faces: "
370 << severeNonOrth <<
"." <<
endl;
378 Info<<
"Mesh non-orthogonality Max: "
386 if (errorNonOrth > 0)
392 "primitiveMeshGeometry::checkFaceDotProduct"
393 "(const bool, const scalar, const labelList&, labelHashSet*)"
394 ) <<
"Error in non-orthogonality detected" <<
endl;
403 Info<<
"Non-orthogonality check OK.\n" <<
endl;
414 const scalar minPyrVol,
428 label nErrorPyrs = 0;
432 label faceI = checkFaces[i];
438 cellCentres[own[faceI]]
441 if (pyrVol > -minPyrVol)
445 Pout<<
"bool primitiveMeshGeometry::checkFacePyramids("
446 <<
"const bool, const scalar, const pointField&"
447 <<
", const labelList&, labelHashSet*): "
448 <<
"face " << faceI <<
" points the wrong way. " <<
endl
449 <<
"Pyramid volume: " << -pyrVol
450 <<
" Face " << f[faceI] <<
" area: " << f[faceI].mag(p)
451 <<
" Owner cell: " << own[faceI] <<
endl
452 <<
"Owner cell vertex labels: "
453 << mesh.
cells()[own[faceI]].labels(f)
472 if (pyrVol < minPyrVol)
476 Pout<<
"bool primitiveMeshGeometry::checkFacePyramids("
477 <<
"const bool, const scalar, const pointField&"
478 <<
", const labelList&, labelHashSet*): "
479 <<
"face " << faceI <<
" points the wrong way. " <<
endl
480 <<
"Pyramid volume: " << -pyrVol
481 <<
" Face " << f[faceI] <<
" area: " << f[faceI].mag(p)
482 <<
" Neighbour cell: " << nei[faceI] <<
endl
483 <<
"Neighbour cell vertex labels: "
484 << mesh.
cells()[nei[faceI]].labels(f)
506 "primitiveMeshGeometry::checkFacePyramids("
507 "const bool, const scalar, const pointField&"
508 ", const labelList&, labelHashSet*)"
509 ) <<
"Error in face pyramids: faces pointing the wrong way!"
519 Info<<
"Face pyramids OK.\n" <<
endl;
530 const scalar internalSkew,
531 const scalar boundarySkew,
552 label faceI = checkFaces[i];
556 scalar dOwn =
mag(faceCentres[faceI] - cellCentres[own[faceI]]);
557 scalar dNei =
mag(faceCentres[faceI] - cellCentres[nei[faceI]]);
559 point faceIntersection =
560 cellCentres[own[faceI]]*dNei/(dOwn+dNei)
561 + cellCentres[nei[faceI]]*dOwn/(dOwn+dNei);
564 mag(faceCentres[faceI] - faceIntersection)
566 mag(cellCentres[nei[faceI]]-cellCentres[own[faceI]])
573 if (skewness > internalSkew)
577 Pout<<
"Severe skewness for face " << faceI
578 <<
" skewness = " << skewness <<
endl;
589 if (skewness > maxSkew)
599 vector faceNormal = faceAreas[faceI];
600 faceNormal /=
mag(faceNormal) + VSMALL;
602 vector dOwn = faceCentres[faceI] - cellCentres[own[faceI]];
604 vector dWall = faceNormal*(faceNormal & dOwn);
606 point faceIntersection = cellCentres[own[faceI]] + dWall;
609 mag(faceCentres[faceI] - faceIntersection)
610 /(2*
mag(dWall) + VSMALL);
615 if (skewness > boundarySkew)
619 Pout<<
"Severe skewness for boundary face " << faceI
620 <<
" skewness = " << skewness <<
endl;
631 if (skewness > maxSkew)
647 "primitiveMeshGeometry::checkFaceSkewness"
648 "(const bool, const scalar, const labelList&, labelHashSet*)"
649 ) <<
"Large face skewness detected. Max skewness = "
651 <<
" percent.\nThis may impair the quality of the result." <<
nl
652 << nWarnSkew <<
" highly skew faces detected."
662 Info<<
"Max skewness = " << 100*maxSkew
663 <<
" percent. Face skewness OK.\n" <<
endl;
674 const scalar warnWeight,
688 scalar minWeight = GREAT;
690 label nWarnWeight = 0;
694 label faceI = checkFaces[i];
698 const point& fc = faceCentres[faceI];
700 scalar dOwn =
mag(faceAreas[faceI] & (fc-cellCentres[own[faceI]]));
701 scalar dNei =
mag(faceAreas[faceI] & (cellCentres[nei[faceI]]-fc));
703 scalar weight =
min(dNei,dOwn)/(dNei+dOwn);
705 if (weight < warnWeight)
709 Pout<<
"Small weighting factor for face " << faceI
710 <<
" weight = " << weight <<
endl;
721 minWeight =
min(minWeight, weight);
728 if (minWeight < warnWeight)
734 "primitiveMeshGeometry::checkFaceWeights"
735 "(const bool, const scalar, const labelList&, labelHashSet*)"
736 ) <<
"Small interpolation weight detected. Min weight = "
737 << minWeight <<
'.' <<
nl
738 << nWarnWeight <<
" faces with small weights detected."
748 Info<<
"Min weight = " << minWeight
749 <<
" percent. Weights OK.\n" <<
endl;
772 if (maxDeg < -SMALL || maxDeg > 180+SMALL)
776 "primitiveMeshGeometry::checkFaceAngles"
777 "(const bool, const scalar, const pointField&, const labelList&"
779 ) <<
"maxDeg should be [0..180] but is now " << maxDeg
787 scalar maxEdgeSin = 0.0;
791 label errorFaceI = -1;
795 label faceI = checkFaces[i];
797 const face& f = fcs[faceI];
799 vector faceNormal = faceAreas[faceI];
800 faceNormal /=
mag(faceNormal) + VSMALL;
804 scalar magEPrev =
mag(ePrev);
805 ePrev /= magEPrev + VSMALL;
813 vector e10(p[f[fp1]] - p[f[fp0]]);
814 scalar magE10 =
mag(e10);
815 e10 /= magE10 + VSMALL;
817 if (magEPrev > SMALL && magE10 > SMALL)
819 vector edgeNormal = ePrev ^ e10;
820 scalar magEdgeNormal =
mag(edgeNormal);
822 if (magEdgeNormal < maxSin)
829 edgeNormal /= magEdgeNormal;
831 if ((edgeNormal & faceNormal) < SMALL)
833 if (faceI != errorFaceI)
845 maxEdgeSin =
max(maxEdgeSin, magEdgeNormal);
860 if (maxEdgeSin > SMALL)
862 scalar maxConcaveDegr =
866 Info<<
"There are " << nConcave
867 <<
" faces with concave angles between consecutive"
868 <<
" edges. Max concave angle = " << maxConcaveDegr
869 <<
" degrees.\n" <<
endl;
873 Info<<
"All angles in faces are convex or less than " << maxDeg
874 <<
" degrees concave.\n" <<
endl;
884 "primitiveMeshGeometry::checkFaceAngles"
885 "(const bool, const scalar, const pointField&"
886 ", const labelList&, labelHashSet*)"
887 ) << nConcave <<
" face points with severe concave angle (> "
888 << maxDeg <<
" deg) found.\n"
1044 const scalar minTwist,
1053 if (minTwist < -1-SMALL || minTwist > 1+SMALL)
1057 "primitiveMeshGeometry::checkFaceTwist"
1058 "(const bool, const scalar, const primitiveMesh&, const pointField&"
1059 ", const labelList&, labelHashSet*)"
1060 ) <<
"minTwist should be [-1..1] but is now " << minTwist
1074 label faceI = checkFaces[i];
1076 const face& f = fcs[faceI];
1078 scalar magArea =
mag(faceAreas[faceI]);
1080 if (f.
size() > 3 && magArea > VSMALL)
1082 const vector nf = faceAreas[faceI] / magArea;
1084 const point& fc = faceCentres[faceI];
1098 scalar magTri =
mag(triArea);
1100 if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
1120 Info<<
"There are " << nWarped
1121 <<
" faces with cosine of the angle"
1122 <<
" between triangle normal and face normal less than "
1123 << minTwist <<
nl <<
endl;
1127 Info<<
"All faces are flat in that the cosine of the angle"
1128 <<
" between triangle normal and face normal less than "
1129 << minTwist <<
nl <<
endl;
1139 "primitiveMeshGeometry::checkFaceTwist"
1140 "(const bool, const scalar, const primitiveMesh&"
1141 ", const pointField&, const labelList&, labelHashSet*)"
1142 ) << nWarped <<
" faces with severe warpage "
1143 <<
"(cosine of the angle between triangle normal and face normal"
1144 <<
" < " << minTwist <<
") found.\n"
1160 const scalar minArea,
1167 label nZeroArea = 0;
1171 label faceI = checkFaces[i];
1173 if (
mag(faceAreas[faceI]) < minArea)
1190 Info<<
"There are " << nZeroArea
1191 <<
" faces with area < " << minArea <<
'.' <<
nl <<
endl;
1195 Info<<
"All faces have area > " << minArea <<
'.' <<
nl <<
endl;
1205 "primitiveMeshGeometry::checkFaceArea"
1206 "(const bool, const scalar, const primitiveMesh&"
1207 ", const pointField&, const labelList&, labelHashSet*)"
1208 ) << nZeroArea <<
" faces with area < " << minArea
1225 const scalar warnDet,
1235 scalar minDet = GREAT;
1236 scalar sumDet = 0.0;
1242 const cell& cFaces = cells[affectedCells[i]];
1245 scalar magAreaSum = 0;
1249 label faceI = cFaces[cFaceI];
1251 scalar magArea =
mag(faceAreas[faceI]);
1253 magAreaSum += magArea;
1254 areaSum += faceAreas[faceI]*(faceAreas[faceI]/magArea);
1257 scalar scaledDet =
det(areaSum/magAreaSum)/0.037037037037037;
1259 minDet =
min(minDet, scaledDet);
1260 sumDet += scaledDet;
1263 if (scaledDet < warnDet)
1270 label faceI = cFaces[cFaceI];
1287 Info<<
"Cell determinant (1 = uniform cube) : average = "
1288 << sumDet / nSumDet <<
" min = " << minDet <<
endl;
1293 Info<<
"There are " << nWarnDet
1294 <<
" cells with determinant < " << warnDet <<
'.' <<
nl
1299 Info<<
"All faces have determinant > " << warnDet <<
'.' <<
nl
1310 "primitiveMeshGeometry::checkCellDeterminant"
1311 "(const bool, const scalar, const primitiveMesh&"
1312 ", const pointField&, const labelList&, const labelList&"
1314 ) << nWarnDet <<
" cells with determinant < " << warnDet
1331 const scalar orthWarn,
1336 return checkFaceDotProduct
1352 const scalar minPyrVol,
1358 return checkFacePyramids
1374 const scalar internalSkew,
1375 const scalar boundarySkew,
1380 return checkFaceSkewness
1398 const scalar warnWeight,
1403 return checkFaceWeights
1420 const scalar maxDeg,
1426 return checkFaceAngles
1465 const scalar minTwist,
1471 return checkFaceTwist
1488 const scalar minArea,
1493 return checkFaceArea
1508 const scalar warnDet,
1514 return checkCellDeterminant