44 void Foam::primitiveMeshGeometry::updateFaceCentresAndAreas
54 label facei = changedFaces[i];
63 faceCentres_[facei] = (1.0/3.0)*(
p[
f[0]] +
p[
f[1]] +
p[
f[2]]);
64 faceAreas_[facei] = 0.5*((
p[
f[1]] -
p[
f[0]])^(
p[
f[2]] -
p[
f[0]]));
93 faceCentres_[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL);
94 faceAreas_[facei] = 0.5*sumN;
100 void Foam::primitiveMeshGeometry::updateCellCentresAndVols
107 UIndirectList<vector>(cellCentres_, changedCells) =
Zero;
108 UIndirectList<scalar>(cellVolumes_, changedCells) = 0.0;
110 const labelList& own = mesh_.faceOwner();
111 const labelList& nei = mesh_.faceNeighbour();
116 UIndirectList<vector>(cEst, changedCells) =
Zero;
118 UIndirectList<scalar>(nCellFaces, changedCells) = 0.0;
122 label facei = changedFaces[i];
123 cEst[own[facei]] += faceCentres_[facei];
124 nCellFaces[own[facei]] += 1;
126 if (mesh_.isInternalFace(facei))
128 cEst[nei[facei]] += faceCentres_[facei];
129 nCellFaces[nei[facei]] += 1;
135 label celli = changedCells[i];
136 cEst[celli] /= nCellFaces[celli];
141 label facei = changedFaces[i];
146 faceAreas_[facei] & (faceCentres_[facei] - cEst[own[facei]]),
151 vector pc = (3.0/4.0)*faceCentres_[facei] + (1.0/4.0)*cEst[own[facei]];
154 cellCentres_[own[facei]] += pyr3Vol*pc;
157 cellVolumes_[own[facei]] += pyr3Vol;
159 if (mesh_.isInternalFace(facei))
164 faceAreas_[facei] & (cEst[nei[facei]] - faceCentres_[facei]),
170 (3.0/4.0)*faceCentres_[facei]
171 + (1.0/4.0)*cEst[nei[facei]];
174 cellCentres_[nei[facei]] += pyr3Vol*pc;
177 cellVolumes_[nei[facei]] += pyr3Vol;
183 label celli = changedCells[i];
185 cellCentres_[celli] /= cellVolumes_[celli];
186 cellVolumes_[celli] *= (1.0/3.0);
196 const labelList& own = mesh_.faceOwner();
197 const labelList& nei = mesh_.faceNeighbour();
203 label facei = changedFaces[i];
205 affectedCells.insert(own[facei]);
207 if (mesh_.isInternalFace(facei))
209 affectedCells.insert(nei[facei]);
212 return affectedCells.toc();
238 faceAreas_ = mesh_.faceAreas();
239 faceCentres_ = mesh_.faceCentres();
240 cellCentres_ = mesh_.cellCentres();
241 cellVolumes_ = mesh_.cellVolumes();
252 updateFaceCentresAndAreas(
p, changedFaces);
254 updateCellCentresAndVols(affectedCells(changedFaces), changedFaces);
261 const scalar orthWarn,
275 const scalar severeNonorthogonalityThreshold =
::cos(
degToRad(orthWarn));
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]
323 <<
"Severe non-orthogonality detected for face "
325 <<
" between cells " << own[facei] <<
" and "
340 if (dDotS < minDDotS)
354 label neiSize = nei.size();
360 if (report && minDDotS < severeNonorthogonalityThreshold)
362 Info<<
"Number of non-orthogonality errors: " << errorNonOrth
363 <<
". Number of severely non-orthogonal faces: "
364 << severeNonOrth <<
"." <<
endl;
372 Info<<
"Mesh non-orthogonality Max: "
379 if (errorNonOrth > 0)
384 <<
"Error in non-orthogonality detected" <<
endl;
392 Info<<
"Non-orthogonality check OK.\n" <<
endl;
402 const scalar minPyrVol,
416 label nErrorPyrs = 0;
420 label facei = checkFaces[i];
426 cellCentres[own[facei]]
429 if (pyrVol > -minPyrVol)
433 Pout<<
"bool primitiveMeshGeometry::checkFacePyramids("
434 <<
"const bool, const scalar, const pointField&"
435 <<
", const labelList&, labelHashSet*): "
436 <<
"face " << facei <<
" points the wrong way. " <<
endl
437 <<
"Pyramid volume: " << -pyrVol
438 <<
" Face " <<
f[facei] <<
" area: " <<
f[facei].mag(
p)
439 <<
" Owner cell: " << own[facei] <<
endl
440 <<
"Owner cell vertex labels: "
460 if (pyrVol < minPyrVol)
464 Pout<<
"bool primitiveMeshGeometry::checkFacePyramids("
465 <<
"const bool, const scalar, const pointField&"
466 <<
", const labelList&, labelHashSet*): "
467 <<
"face " << facei <<
" points the wrong way. " <<
endl
468 <<
"Pyramid volume: " << -pyrVol
469 <<
" Face " <<
f[facei] <<
" area: " <<
f[facei].mag(
p)
470 <<
" Neighbour cell: " << nei[facei] <<
endl
471 <<
"Neighbour cell vertex labels: "
493 <<
"Error in face pyramids: faces pointing the wrong way!"
502 Info<<
"Face pyramids OK.\n" <<
endl;
512 const scalar internalSkew,
513 const scalar boundarySkew,
534 label facei = checkFaces[i];
538 scalar dOwn =
mag(faceCentres[facei] - cellCentres[own[facei]]);
539 scalar dNei =
mag(faceCentres[facei] - cellCentres[nei[facei]]);
541 point faceIntersection =
542 cellCentres[own[facei]]*dNei/(dOwn+dNei)
543 + cellCentres[nei[facei]]*dOwn/(dOwn+dNei);
546 mag(faceCentres[facei] - faceIntersection)
548 mag(cellCentres[nei[facei]]-cellCentres[own[facei]])
555 if (skewness > internalSkew)
559 Pout<<
"Severe skewness for face " << facei
560 <<
" skewness = " << skewness <<
endl;
571 if (skewness > maxSkew)
583 vector dOwn = faceCentres[facei] - cellCentres[own[facei]];
585 vector dWall = faceNormal*(faceNormal & dOwn);
587 point faceIntersection = cellCentres[own[facei]] + dWall;
590 mag(faceCentres[facei] - faceIntersection)
591 /(2*
mag(dWall) + VSMALL);
596 if (skewness > boundarySkew)
600 Pout<<
"Severe skewness for boundary face " << facei
601 <<
" skewness = " << skewness <<
endl;
612 if (skewness > maxSkew)
628 <<
" percent.\nThis may impair the quality of the result." <<
nl
629 << nWarnSkew <<
" highly skew faces detected."
638 Info<<
"Max skewness = " << 100*maxSkew
639 <<
" percent. Face skewness OK.\n" <<
endl;
649 const scalar warnWeight,
663 scalar minWeight = GREAT;
665 label nWarnWeight = 0;
669 label facei = checkFaces[i];
673 const point& fc = faceCentres[facei];
675 scalar dOwn =
mag(faceAreas[facei] & (fc-cellCentres[own[facei]]));
676 scalar dNei =
mag(faceAreas[facei] & (cellCentres[nei[facei]]-fc));
678 scalar weight =
min(dNei,dOwn)/(dNei+dOwn);
680 if (weight < warnWeight)
684 Pout<<
"Small weighting factor for face " << facei
685 <<
" weight = " << weight <<
endl;
696 minWeight =
min(minWeight, weight);
703 if (minWeight < warnWeight)
708 << minWeight <<
'.' <<
nl
709 << nWarnWeight <<
" faces with small weights detected."
718 Info<<
"Min weight = " << minWeight
719 <<
" percent. Weights OK.\n" <<
endl;
741 if (maxDeg < -SMALL || maxDeg > 180+SMALL)
744 <<
"maxDeg should be [0..180] but is now " << maxDeg
752 scalar maxEdgeSin = 0.0;
756 label errorFacei = -1;
760 label facei = checkFaces[i];
762 const face&
f = fcs[facei];
768 scalar magEPrev =
mag(ePrev);
769 ePrev /= magEPrev + VSMALL;
774 label fp1 =
f.fcIndex(fp0);
778 scalar magE10 =
mag(e10);
779 e10 /= magE10 + VSMALL;
781 if (magEPrev > SMALL && magE10 > SMALL)
783 vector edgeNormal = ePrev ^ e10;
784 scalar magEdgeNormal =
mag(edgeNormal);
786 if (magEdgeNormal < maxSin)
793 edgeNormal /= magEdgeNormal;
795 if ((edgeNormal & faceNormal) < SMALL)
797 if (facei != errorFacei)
809 maxEdgeSin =
max(maxEdgeSin, magEdgeNormal);
824 if (maxEdgeSin > SMALL)
826 scalar maxConcaveDegr =
829 Info<<
"There are " << nConcave
830 <<
" faces with concave angles between consecutive"
831 <<
" edges. Max concave angle = " << maxConcaveDegr
832 <<
" degrees.\n" <<
endl;
836 Info<<
"All angles in faces are convex or less than " << maxDeg
837 <<
" degrees concave.\n" <<
endl;
846 << nConcave <<
" face points with severe concave angle (> "
847 << maxDeg <<
" deg) found.\n"
991 const scalar minTwist,
1000 if (minTwist < -1-SMALL || minTwist > 1+SMALL)
1003 <<
"minTwist should be [-1..1] but is now " << minTwist
1015 for (
const label facei : checkFaces)
1017 const face&
f = fcs[facei];
1019 const scalar magArea =
mag(faceAreas[facei]);
1021 if (
f.size() > 3 && magArea > VSMALL)
1023 const vector nf = faceAreas[facei] / magArea;
1025 const point& fc = faceCentres[facei];
1034 p[
f.nextLabel(fpI)],
1039 const scalar magTri =
mag(triArea);
1041 if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
1061 Info<<
"There are " << nWarped
1062 <<
" faces with cosine of the angle"
1063 <<
" between triangle normal and face normal less than "
1064 << minTwist <<
nl <<
endl;
1068 Info<<
"All faces are flat in that the cosine of the angle"
1069 <<
" between triangle normal and face normal less than "
1070 << minTwist <<
nl <<
endl;
1079 << nWarped <<
" faces with severe warpage "
1080 <<
"(cosine of the angle between triangle normal and "
1081 <<
"face normal < " << minTwist <<
") found.\n"
1095 const scalar minArea,
1102 label nZeroArea = 0;
1106 label facei = checkFaces[i];
1108 if (
mag(faceAreas[facei]) < minArea)
1125 Info<<
"There are " << nZeroArea
1126 <<
" faces with area < " << minArea <<
'.' <<
nl <<
endl;
1130 Info<<
"All faces have area > " << minArea <<
'.' <<
nl <<
endl;
1139 << nZeroArea <<
" faces with area < " << minArea
1154 const scalar warnDet,
1164 scalar minDet = GREAT;
1165 scalar sumDet = 0.0;
1171 const cell& cFaces =
cells[affectedCells[i]];
1174 scalar magAreaSum = 0;
1178 label facei = cFaces[cFacei];
1180 scalar magArea =
mag(faceAreas[facei]);
1182 magAreaSum += magArea;
1183 areaSum += faceAreas[facei]*(faceAreas[facei]/magArea);
1186 scalar scaledDet =
det(areaSum/magAreaSum)/0.037037037037037;
1188 minDet =
min(minDet, scaledDet);
1189 sumDet += scaledDet;
1192 if (scaledDet < warnDet)
1199 label facei = cFaces[cFacei];
1216 Info<<
"Cell determinant (1 = uniform cube) : average = "
1217 << sumDet / nSumDet <<
" min = " << minDet <<
endl;
1222 Info<<
"There are " << nWarnDet
1223 <<
" cells with determinant < " << warnDet <<
'.' <<
nl
1228 Info<<
"All faces have determinant > " << warnDet <<
'.' <<
nl
1238 << nWarnDet <<
" cells with determinant < " << warnDet
1253 const scalar orthWarn,
1258 return checkFaceDotProduct
1274 const scalar minPyrVol,
1280 return checkFacePyramids
1296 const scalar internalSkew,
1297 const scalar boundarySkew,
1302 return checkFaceSkewness
1320 const scalar warnWeight,
1325 return checkFaceWeights
1342 const scalar maxDeg,
1348 return checkFaceAngles
1387 const scalar minTwist,
1393 return checkFaceTwist
1410 const scalar minArea,
1415 return checkFaceArea
1430 const scalar warnDet,
1436 return checkCellDeterminant