44void Foam::primitiveMeshGeometry::updateFaceCentresAndAreas
47 const labelList& changedFaces
52 for (label facei : changedFaces)
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]]));
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;
98void Foam::primitiveMeshGeometry::updateCellCentresAndVols
105 UIndirectList<vector>(cellCentres_, changedCells) =
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) =
Zero;
116 UIndirectList<scalar>(nCellFaces, changedCells) = 0.0;
118 for (label facei : changedFaces)
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;
130 for (label celli : changedCells)
132 cEst[celli] /= nCellFaces[celli];
135 for (label facei : changedFaces)
140 faceAreas_[facei] & (faceCentres_[facei] - cEst[own[facei]]),
145 vector pc = (3.0/4.0)*faceCentres_[facei] + (1.0/4.0)*cEst[own[facei]];
148 cellCentres_[own[facei]] += pyr3Vol*pc;
151 cellVolumes_[own[facei]] += pyr3Vol;
153 if (mesh_.isInternalFace(facei))
158 faceAreas_[facei] & (cEst[nei[facei]] - faceCentres_[facei]),
164 (3.0/4.0)*faceCentres_[facei]
165 + (1.0/4.0)*cEst[nei[facei]];
168 cellCentres_[nei[facei]] += pyr3Vol*pc;
171 cellVolumes_[nei[facei]] += pyr3Vol;
177 label celli = changedCells[i];
179 cellCentres_[celli] /= cellVolumes_[celli];
180 cellVolumes_[celli] *= (1.0/3.0);
190 const labelList& own = mesh_.faceOwner();
191 const labelList& nei = mesh_.faceNeighbour();
195 for (label facei : changedFaces)
197 affectedCells.
insert(own[facei]);
199 if (mesh_.isInternalFace(facei))
201 affectedCells.
insert(nei[facei]);
204 return affectedCells.
toc();
225 faceAreas_ = mesh_.faceAreas();
226 faceCentres_ = mesh_.faceCentres();
227 cellCentres_ = mesh_.cellCentres();
228 cellVolumes_ = mesh_.cellVolumes();
239 updateFaceCentresAndAreas(
p, changedFaces);
241 updateCellCentresAndVols(affectedCells(changedFaces), changedFaces);
248 const scalar orthWarn,
262 const scalar severeNonorthogonalityThreshold =
::cos(
degToRad(orthWarn));
264 scalar minDDotS = GREAT;
268 label severeNonOrth = 0;
270 label errorNonOrth = 0;
274 label facei = checkFaces[i];
278 vector d = cellCentres[nei[facei]] - cellCentres[own[facei]];
279 const vector&
s = faceAreas[facei];
281 scalar dDotS = (d &
s)/(
mag(d)*
mag(
s) + VSMALL);
283 if (dDotS < severeNonorthogonalityThreshold)
290 Pout<<
"Severe non-orthogonality for face " << facei
291 <<
" between cells " << own[facei]
292 <<
" and " << nei[facei]
310 <<
"Severe non-orthogonality detected for face "
312 <<
" between cells " << own[facei] <<
" and "
327 if (dDotS < minDDotS)
341 label neiSize = nei.
size();
347 if (report && minDDotS < severeNonorthogonalityThreshold)
349 Info<<
"Number of non-orthogonality errors: " << errorNonOrth
350 <<
". Number of severely non-orthogonal faces: "
351 << severeNonOrth <<
"." <<
endl;
359 Info<<
"Mesh non-orthogonality Max: "
366 if (errorNonOrth > 0)
371 <<
"Error in non-orthogonality detected" <<
endl;
379 Info<<
"Non-orthogonality check OK.\n" <<
endl;
389 const scalar minPyrVol,
403 label nErrorPyrs = 0;
405 for (label facei : checkFaces)
411 cellCentres[own[facei]]
414 if (pyrVol > -minPyrVol)
418 Pout<<
"bool primitiveMeshGeometry::checkFacePyramids("
419 <<
"const bool, const scalar, const pointField&"
420 <<
", const labelList&, labelHashSet*): "
421 <<
"face " << facei <<
" points the wrong way. " <<
endl
422 <<
"Pyramid volume: " << -pyrVol
423 <<
" Face " <<
f[facei] <<
" area: " <<
f[facei].mag(
p)
424 <<
" Owner cell: " << own[facei] <<
endl
425 <<
"Owner cell vertex labels: "
445 if (pyrVol < minPyrVol)
449 Pout<<
"bool primitiveMeshGeometry::checkFacePyramids("
450 <<
"const bool, const scalar, const pointField&"
451 <<
", const labelList&, labelHashSet*): "
452 <<
"face " << facei <<
" points the wrong way. " <<
endl
453 <<
"Pyramid volume: " << -pyrVol
454 <<
" Face " <<
f[facei] <<
" area: " <<
f[facei].mag(
p)
455 <<
" Neighbour cell: " << nei[facei] <<
endl
456 <<
"Neighbour cell vertex labels: "
478 <<
"Error in face pyramids: faces pointing the wrong way!"
487 Info<<
"Face pyramids OK.\n" <<
endl;
497 const scalar internalSkew,
498 const scalar boundarySkew,
517 for (label facei : checkFaces)
521 scalar dOwn =
mag(faceCentres[facei] - cellCentres[own[facei]]);
522 scalar dNei =
mag(faceCentres[facei] - cellCentres[nei[facei]]);
524 point faceIntersection =
525 cellCentres[own[facei]]*dNei/(dOwn+dNei)
526 + cellCentres[nei[facei]]*dOwn/(dOwn+dNei);
529 mag(faceCentres[facei] - faceIntersection)
531 mag(cellCentres[nei[facei]]-cellCentres[own[facei]])
538 if (skewness > internalSkew)
542 Pout<<
"Severe skewness for face " << facei
543 <<
" skewness = " << skewness <<
endl;
554 if (skewness > maxSkew)
566 vector dOwn = faceCentres[facei] - cellCentres[own[facei]];
568 vector dWall = faceNormal*(faceNormal & dOwn);
570 point faceIntersection = cellCentres[own[facei]] + dWall;
573 mag(faceCentres[facei] - faceIntersection)
574 /(2*
mag(dWall) + VSMALL);
579 if (skewness > boundarySkew)
583 Pout<<
"Severe skewness for boundary face " << facei
584 <<
" skewness = " << skewness <<
endl;
595 if (skewness > maxSkew)
611 <<
" percent.\nThis may impair the quality of the result." <<
nl
612 << nWarnSkew <<
" highly skew faces detected."
621 Info<<
"Max skewness = " << 100*maxSkew
622 <<
" percent. Face skewness OK.\n" <<
endl;
632 const scalar warnWeight,
646 scalar minWeight = GREAT;
648 label nWarnWeight = 0;
650 for (label facei : checkFaces)
654 const point& fc = faceCentres[facei];
656 scalar dOwn =
mag(faceAreas[facei] & (fc-cellCentres[own[facei]]));
657 scalar dNei =
mag(faceAreas[facei] & (cellCentres[nei[facei]]-fc));
659 scalar weight =
min(dNei,dOwn)/(dNei+dOwn);
661 if (weight < warnWeight)
665 Pout<<
"Small weighting factor for face " << facei
666 <<
" weight = " << weight <<
endl;
677 minWeight =
min(minWeight, weight);
684 if (minWeight < warnWeight)
689 << minWeight <<
'.' <<
nl
690 << nWarnWeight <<
" faces with small weights detected."
699 Info<<
"Min weight = " << minWeight
700 <<
" percent. Weights OK.\n" <<
endl;
722 if (maxDeg < -SMALL || maxDeg > 180+SMALL)
725 <<
"maxDeg should be [0..180] but is now " << maxDeg
733 scalar maxEdgeSin = 0.0;
737 label errorFacei = -1;
739 for (label facei : checkFaces)
741 const face&
f = fcs[facei];
747 scalar magEPrev =
mag(ePrev);
748 ePrev /= magEPrev + VSMALL;
753 vector e10(
p[
f.nextLabel(fp0)] -
p[
f.thisLabel(fp0)]);
754 scalar magE10 =
mag(e10);
755 e10 /= magE10 + VSMALL;
757 if (magEPrev > SMALL && magE10 > SMALL)
759 vector edgeNormal = ePrev ^ e10;
760 scalar magEdgeNormal =
mag(edgeNormal);
762 if (magEdgeNormal < maxSin)
769 edgeNormal /= magEdgeNormal;
771 if ((edgeNormal & faceNormal) < SMALL)
773 if (facei != errorFacei)
785 maxEdgeSin =
max(maxEdgeSin, magEdgeNormal);
800 if (maxEdgeSin > SMALL)
802 scalar maxConcaveDegr =
805 Info<<
"There are " << nConcave
806 <<
" faces with concave angles between consecutive"
807 <<
" edges. Max concave angle = " << maxConcaveDegr
808 <<
" degrees.\n" <<
endl;
812 Info<<
"All angles in faces are convex or less than " << maxDeg
813 <<
" degrees concave.\n" <<
endl;
822 << nConcave <<
" face points with severe concave angle (> "
823 << maxDeg <<
" deg) found.\n"
967 const scalar minTwist,
976 if (minTwist < -1-SMALL || minTwist > 1+SMALL)
979 <<
"minTwist should be [-1..1] but is now " << minTwist
991 for (
const label facei : checkFaces)
993 const face&
f = fcs[facei];
995 const scalar magArea =
mag(faceAreas[facei]);
997 if (
f.
size() > 3 && magArea > VSMALL)
999 const vector nf = faceAreas[facei] / magArea;
1001 const point& fc = faceCentres[facei];
1010 p[
f.nextLabel(fpI)],
1015 const scalar magTri =
mag(triArea);
1017 if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
1037 Info<<
"There are " << nWarped
1038 <<
" faces with cosine of the angle"
1039 <<
" between triangle normal and face normal less than "
1040 << minTwist <<
nl <<
endl;
1044 Info<<
"All faces are flat in that the cosine of the angle"
1045 <<
" between triangle normal and face normal less than "
1046 << minTwist <<
nl <<
endl;
1055 << nWarped <<
" faces with severe warpage "
1056 <<
"(cosine of the angle between triangle normal and "
1057 <<
"face normal < " << minTwist <<
") found.\n"
1071 const scalar minArea,
1078 label nZeroArea = 0;
1080 for (label facei : checkFaces)
1082 if (
mag(faceAreas[facei]) < minArea)
1099 Info<<
"There are " << nZeroArea
1100 <<
" faces with area < " << minArea <<
'.' <<
nl <<
endl;
1104 Info<<
"All faces have area > " << minArea <<
'.' <<
nl <<
endl;
1113 << nZeroArea <<
" faces with area < " << minArea
1128 const scalar warnDet,
1138 scalar minDet = GREAT;
1139 scalar sumDet = 0.0;
1143 for (label celli : affectedCells)
1148 scalar magAreaSum = 0;
1152 label facei = cFaces[cFacei];
1154 scalar magArea =
mag(faceAreas[facei]);
1156 magAreaSum += magArea;
1157 areaSum += faceAreas[facei]*(faceAreas[facei]/magArea);
1160 scalar scaledDet =
det(areaSum/magAreaSum)/0.037037037037037;
1162 minDet =
min(minDet, scaledDet);
1163 sumDet += scaledDet;
1166 if (scaledDet < warnDet)
1173 label facei = cFaces[cFacei];
1190 Info<<
"Cell determinant (1 = uniform cube) : average = "
1191 << sumDet / nSumDet <<
" min = " << minDet <<
endl;
1196 Info<<
"There are " << nWarnDet
1197 <<
" cells with determinant < " << warnDet <<
'.' <<
nl
1202 Info<<
"All faces have determinant > " << warnDet <<
'.' <<
nl
1212 << nWarnDet <<
" cells with determinant < " << warnDet
1227 const scalar orthWarn,
1232 return checkFaceDotProduct
1248 const scalar minPyrVol,
1254 return checkFacePyramids
1270 const scalar internalSkew,
1271 const scalar boundarySkew,
1276 return checkFaceSkewness
1294 const scalar warnWeight,
1299 return checkFaceWeights
1316 const scalar maxDeg,
1322 return checkFaceAngles
1361 const scalar minTwist,
1367 return checkFaceTwist
1384 const scalar minArea,
1389 return checkFaceArea
1404 const scalar warnDet,
1410 return checkCellDeterminant
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
List< Key > toc() const
The table of contents (the keys) in unsorted order.
T & first()
Return the first element of the list.
void size(const label n)
Older name for setAddressableSize.
T & last()
Return the last element of the list.
A cell is defined as a list of faces with extra functionality.
A face is a list of labels corresponding to mesh vertices.
virtual const faceList & faces() const
Return raw faces.
virtual const labelList & faceOwner() const
Return face owner.
virtual const labelList & faceNeighbour() const
Return face neighbour.
Updateable mesh geometry + checking routines.
static bool checkFaceTwist(const bool report, const scalar minTwist, const primitiveMesh &, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
static bool checkFaceAngles(const bool report, const scalar maxDeg, const primitiveMesh &mesh, const vectorField &faceAreas, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
static bool checkFaceSkewness(const bool report, const scalar internalSkew, const scalar boundarySkew, const primitiveMesh &mesh, const vectorField &cellCentres, const vectorField &faceCentres, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
void correct()
Take over properties from mesh.
static bool checkFaceWeights(const bool report, const scalar warnWeight, const primitiveMesh &mesh, const vectorField &cellCentres, const vectorField &faceCentres, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
labelList affectedCells(const labelList &changedFaces) const
Helper function: get affected cells from faces.
static bool checkCellDeterminant(const bool report, const scalar minDet, const primitiveMesh &, const vectorField &faceAreas, const labelList &checkFaces, const labelList &affectedCells, labelHashSet *setPtr)
static bool checkFaceArea(const bool report, const scalar minArea, const primitiveMesh &, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
static bool checkFaceDotProduct(const bool report, const scalar orthWarn, const primitiveMesh &, const vectorField &cellCentres, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
static bool checkFacePyramids(const bool report, const scalar minPyrVol, const primitiveMesh &, const vectorField &cellCentres, const pointField &p, const labelList &checkFaces, labelHashSet *)
Cell-face mesh analysis engine.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
virtual const faceList & faces() const =0
Return faces.
const cellList & cells() const
A triangle primitive used to calculate face normals and swept volumes.
vector areaNormal() const
The area normal - with magnitude equal to area of triangle.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
#define WarningInFunction
Report a warning using Foam::Warning.
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
constexpr scalar pi(M_PI)
const dimensionedScalar c
Speed of light in a vacuum.
pyramid< point, const point &, const face & > pyramidPointFaceRef
A pyramid using referred points and faces.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedScalar asin(const dimensionedScalar &ds)
List< label > labelList
A List of labels.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
dimensionedScalar sin(const dimensionedScalar &ds)
messageStream Info
Information stream (stdout output on master, null elsewhere)
vector point
Point is a vector.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
constexpr scalar radToDeg() noexcept
Multiplication factor for radians to degrees conversion.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
Field< vector > vectorField
Specialisation of Field<T> for vector.
static constexpr const zero Zero
Global zero (0)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< face > faceList
A List of faces.
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar acos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.
Unit conversion functions.