44 void Foam::primitiveMeshGeometry::updateFaceCentresAndAreas
52 for (label facei : changedFaces)
94 faceCentres_[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL);
95 faceAreas_[facei] = 0.5*sumN;
101 void Foam::primitiveMeshGeometry::updateCellCentresAndVols
108 UIndirectList<vector>(cellCentres_, changedCells) =
Zero;
109 UIndirectList<scalar>(cellVolumes_, changedCells) = 0.0;
111 const labelList& own = mesh_.faceOwner();
112 const labelList& nei = mesh_.faceNeighbour();
117 UIndirectList<vector>(cEst, changedCells) =
Zero;
119 UIndirectList<scalar>(nCellFaces, changedCells) = 0.0;
121 for (label facei : changedFaces)
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;
133 for (label celli : changedCells)
135 cEst[celli] /= nCellFaces[celli];
138 for (label facei : changedFaces)
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();
198 for (label facei : changedFaces)
200 affectedCells.insert(own[facei]);
202 if (mesh_.isInternalFace(facei))
204 affectedCells.insert(nei[facei]);
207 return affectedCells.toc();
228 faceAreas_ = mesh_.faceAreas();
229 faceCentres_ = mesh_.faceCentres();
230 cellCentres_ = mesh_.cellCentres();
231 cellVolumes_ = mesh_.cellVolumes();
242 updateFaceCentresAndAreas(
p, changedFaces);
244 updateCellCentresAndVols(affectedCells(changedFaces), changedFaces);
251 const scalar orthWarn,
265 const scalar severeNonorthogonalityThreshold =
::cos(
degToRad(orthWarn));
267 scalar minDDotS = GREAT;
271 label severeNonOrth = 0;
273 label errorNonOrth = 0;
277 label facei = checkFaces[i];
281 vector d = cellCentres[nei[facei]] - cellCentres[own[facei]];
282 const vector&
s = faceAreas[facei];
284 scalar dDotS = (d &
s)/(
mag(d)*
mag(
s) + VSMALL);
286 if (dDotS < severeNonorthogonalityThreshold)
293 Pout<<
"Severe non-orthogonality for face " << facei
294 <<
" between cells " << own[facei]
295 <<
" and " << nei[facei]
313 <<
"Severe non-orthogonality detected for face " 315 <<
" between cells " << own[facei] <<
" and " 330 if (dDotS < minDDotS)
339 reduce(minDDotS, minOp<scalar>());
340 reduce(sumDDotS, sumOp<scalar>());
341 reduce(severeNonOrth, sumOp<label>());
342 reduce(errorNonOrth, sumOp<label>());
349 if (report && minDDotS < severeNonorthogonalityThreshold)
351 Info<<
"Number of non-orthogonality errors: " << errorNonOrth
352 <<
". Number of severely non-orthogonal faces: " 353 << severeNonOrth <<
"." <<
endl;
361 Info<<
"Mesh non-orthogonality Max: " 368 if (errorNonOrth > 0)
373 <<
"Error in non-orthogonality detected" <<
endl;
381 Info<<
"Non-orthogonality check OK.\n" <<
endl;
391 const scalar minPyrVol,
392 const primitiveMesh&
mesh,
405 label nErrorPyrs = 0;
407 for (label facei : checkFaces)
413 cellCentres[own[facei]]
416 if (pyrVol > -minPyrVol)
420 Pout<<
"bool primitiveMeshGeometry::checkFacePyramids(" 421 <<
"const bool, const scalar, const pointField&" 422 <<
", const labelList&, labelHashSet*): " 423 <<
"face " << facei <<
" points the wrong way. " <<
endl 424 <<
"Pyramid volume: " << -pyrVol
425 <<
" Face " <<
f[facei] <<
" area: " <<
f[facei].mag(
p)
426 <<
" Owner cell: " << own[facei] <<
endl 427 <<
"Owner cell vertex labels: " 435 setPtr->insert(facei);
447 if (pyrVol < minPyrVol)
451 Pout<<
"bool primitiveMeshGeometry::checkFacePyramids(" 452 <<
"const bool, const scalar, const pointField&" 453 <<
", const labelList&, labelHashSet*): " 454 <<
"face " << facei <<
" points the wrong way. " <<
endl 455 <<
"Pyramid volume: " << -pyrVol
456 <<
" Face " <<
f[facei] <<
" area: " <<
f[facei].mag(
p)
457 <<
" Neighbour cell: " << nei[facei] <<
endl 458 <<
"Neighbour cell vertex labels: " 465 setPtr->insert(facei);
473 reduce(nErrorPyrs, sumOp<label>());
480 <<
"Error in face pyramids: faces pointing the wrong way!" 489 Info<<
"Face pyramids OK.\n" <<
endl;
499 const scalar internalSkew,
500 const scalar boundarySkew,
501 const primitiveMesh&
mesh,
519 for (label facei : checkFaces)
523 scalar dOwn =
mag(faceCentres[facei] - cellCentres[own[facei]]);
524 scalar dNei =
mag(faceCentres[facei] - cellCentres[nei[facei]]);
526 point faceIntersection =
527 cellCentres[own[facei]]*dNei/(dOwn+dNei)
528 + cellCentres[nei[facei]]*dOwn/(dOwn+dNei);
531 mag(faceCentres[facei] - faceIntersection)
533 mag(cellCentres[nei[facei]]-cellCentres[own[facei]])
540 if (skewness > internalSkew)
544 Pout<<
"Severe skewness for face " << facei
545 <<
" skewness = " << skewness <<
endl;
550 setPtr->insert(facei);
556 if (skewness > maxSkew)
568 vector dOwn = faceCentres[facei] - cellCentres[own[facei]];
570 vector dWall = faceNormal*(faceNormal & dOwn);
572 point faceIntersection = cellCentres[own[facei]] + dWall;
575 mag(faceCentres[facei] - faceIntersection)
576 /(2*
mag(dWall) + VSMALL);
581 if (skewness > boundarySkew)
585 Pout<<
"Severe skewness for boundary face " << facei
586 <<
" skewness = " << skewness <<
endl;
591 setPtr->insert(facei);
597 if (skewness > maxSkew)
604 reduce(maxSkew, maxOp<scalar>());
605 reduce(nWarnSkew, sumOp<label>());
613 <<
" percent.\nThis may impair the quality of the result." <<
nl 614 << nWarnSkew <<
" highly skew faces detected." 623 Info<<
"Max skewness = " << 100*maxSkew
624 <<
" percent. Face skewness OK.\n" <<
endl;
634 const scalar warnWeight,
635 const primitiveMesh&
mesh,
648 scalar minWeight = GREAT;
650 label nWarnWeight = 0;
652 for (label facei : checkFaces)
656 const point& fc = faceCentres[facei];
658 scalar dOwn =
mag(faceAreas[facei] & (fc-cellCentres[own[facei]]));
659 scalar dNei =
mag(faceAreas[facei] & (cellCentres[nei[facei]]-fc));
661 scalar weight =
min(dNei,dOwn)/(dNei+dOwn);
663 if (weight < warnWeight)
667 Pout<<
"Small weighting factor for face " << facei
668 <<
" weight = " << weight <<
endl;
673 setPtr->insert(facei);
679 minWeight =
min(minWeight, weight);
683 reduce(minWeight, minOp<scalar>());
684 reduce(nWarnWeight, sumOp<label>());
686 if (minWeight < warnWeight)
691 << minWeight <<
'.' <<
nl 692 << nWarnWeight <<
" faces with small weights detected." 701 Info<<
"Min weight = " << minWeight
702 <<
" percent. Weights OK.\n" <<
endl;
724 if (maxDeg < -SMALL || maxDeg > 180+SMALL)
727 <<
"maxDeg should be [0..180] but is now " << maxDeg
735 scalar maxEdgeSin = 0.0;
739 label errorFacei = -1;
741 for (label facei : checkFaces)
743 const face&
f = fcs[facei];
749 scalar magEPrev =
mag(ePrev);
750 ePrev /= magEPrev + VSMALL;
755 vector e10(
p[
f.nextLabel(fp0)] -
p[
f.thisLabel(fp0)]);
756 scalar magE10 =
mag(e10);
757 e10 /= magE10 + VSMALL;
759 if (magEPrev > SMALL && magE10 > SMALL)
761 vector edgeNormal = ePrev ^ e10;
762 scalar magEdgeNormal =
mag(edgeNormal);
764 if (magEdgeNormal < maxSin)
771 edgeNormal /= magEdgeNormal;
773 if ((edgeNormal & faceNormal) < SMALL)
775 if (facei != errorFacei)
787 maxEdgeSin =
max(maxEdgeSin, magEdgeNormal);
802 if (maxEdgeSin > SMALL)
804 scalar maxConcaveDegr =
807 Info<<
"There are " << nConcave
808 <<
" faces with concave angles between consecutive" 809 <<
" edges. Max concave angle = " << maxConcaveDegr
810 <<
" degrees.\n" <<
endl;
814 Info<<
"All angles in faces are convex or less than " << maxDeg
815 <<
" degrees concave.\n" <<
endl;
824 << nConcave <<
" face points with severe concave angle (> " 825 << maxDeg <<
" deg) found.\n" 969 const scalar minTwist,
978 if (minTwist < -1-SMALL || minTwist > 1+SMALL)
981 <<
"minTwist should be [-1..1] but is now " << minTwist
993 for (
const label facei : checkFaces)
995 const face&
f = fcs[facei];
997 const scalar magArea =
mag(faceAreas[facei]);
999 if (
f.
size() > 3 && magArea > VSMALL)
1001 const vector nf = faceAreas[facei] / magArea;
1003 const point& fc = faceCentres[facei];
1012 p[
f.nextLabel(fpI)],
1017 const scalar magTri =
mag(triArea);
1019 if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
1039 Info<<
"There are " << nWarped
1040 <<
" faces with cosine of the angle" 1041 <<
" between triangle normal and face normal less than " 1042 << minTwist <<
nl <<
endl;
1046 Info<<
"All faces are flat in that the cosine of the angle" 1047 <<
" between triangle normal and face normal less than " 1048 << minTwist <<
nl <<
endl;
1057 << nWarped <<
" faces with severe warpage " 1058 <<
"(cosine of the angle between triangle normal and " 1059 <<
"face normal < " << minTwist <<
") found.\n" 1073 const scalar minArea,
1074 const primitiveMesh&
mesh,
1080 label nZeroArea = 0;
1082 for (label facei : checkFaces)
1084 if (
mag(faceAreas[facei]) < minArea)
1088 setPtr->insert(facei);
1095 reduce(nZeroArea, sumOp<label>());
1101 Info<<
"There are " << nZeroArea
1102 <<
" faces with area < " << minArea <<
'.' <<
nl <<
endl;
1106 Info<<
"All faces have area > " << minArea <<
'.' <<
nl <<
endl;
1115 << nZeroArea <<
" faces with area < " << minArea
1130 const scalar warnDet,
1131 const primitiveMesh&
mesh,
1140 scalar minDet = GREAT;
1141 scalar sumDet = 0.0;
1145 for (label celli : affectedCells)
1147 const cell& cFaces =
cells[celli];
1150 scalar magAreaSum = 0;
1154 label facei = cFaces[cFacei];
1156 scalar magArea =
mag(faceAreas[facei]);
1158 magAreaSum += magArea;
1159 areaSum += faceAreas[facei]*(faceAreas[facei]/magArea);
1162 scalar scaledDet =
det(areaSum/magAreaSum)/0.037037037037037;
1164 minDet =
min(minDet, scaledDet);
1165 sumDet += scaledDet;
1168 if (scaledDet < warnDet)
1175 label facei = cFaces[cFacei];
1176 setPtr->insert(facei);
1183 reduce(minDet, minOp<scalar>());
1184 reduce(sumDet, sumOp<scalar>());
1185 reduce(nSumDet, sumOp<label>());
1186 reduce(nWarnDet, sumOp<label>());
1192 Info<<
"Cell determinant (1 = uniform cube) : average = " 1193 << sumDet / nSumDet <<
" min = " << minDet <<
endl;
1198 Info<<
"There are " << nWarnDet
1199 <<
" cells with determinant < " << warnDet <<
'.' <<
nl 1204 Info<<
"All faces have determinant > " << warnDet <<
'.' <<
nl 1214 << nWarnDet <<
" cells with determinant < " << warnDet
1229 const scalar orthWarn,
1234 return checkFaceDotProduct
1250 const scalar minPyrVol,
1256 return checkFacePyramids
1272 const scalar internalSkew,
1273 const scalar boundarySkew,
1278 return checkFaceSkewness
1296 const scalar warnWeight,
1301 return checkFaceWeights
1318 const scalar maxDeg,
1324 return checkFaceAngles
1363 const scalar minTwist,
1369 return checkFaceTwist
1386 const scalar minArea,
1391 return checkFaceArea
1406 const scalar warnDet,
1412 return checkCellDeterminant
void size(const label n)
Older name for setAddressableSize.
dimensionedScalar acos(const dimensionedScalar &ds)
List< cell > cellList
List of cell.
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
virtual const labelList & faceNeighbour() const
Return face neighbour.
Cell-face mesh analysis engine.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Unit conversion functions.
constexpr char nl
The newline '\n' character (0x0a)
T & first()
Access first element of the list, position [0].
Ostream & endl(Ostream &os)
Add newline and flush stream.
const cellList & cells() const
dimensionedScalar det(const dimensionedSphericalTensor &dt)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
static bool checkFacePyramids(const bool report, const scalar minPyrVol, const primitiveMesh &, const vectorField &cellCentres, const pointField &p, const labelList &checkFaces, labelHashSet *)
dimensionedScalar asin(const dimensionedScalar &ds)
labelList affectedCells(const labelList &changedFaces) const
Helper function: get affected cells from faces.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
#define forAll(list, i)
Loop across all elements in list.
pyramid< point, const point &, const face & > pyramidPointFaceRef
A pyramid using referred point and face.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
static bool checkFaceArea(const bool report, const scalar minArea, const primitiveMesh &, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
List< face > faceList
List of 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 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)
vectorField pointField
pointField is a vectorField.
dimensionedScalar cos(const dimensionedScalar &ds)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static bool checkFaceAngles(const bool report, const scalar maxDeg, const primitiveMesh &mesh, const vectorField &faceAreas, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
virtual const labelList & faceOwner() const
Return face owner.
Point centre() const
Return centre (centroid)
constexpr scalar pi(M_PI)
primitiveMeshGeometry(const primitiveMesh &)
Construct from mesh.
virtual const faceList & faces() const
Return raw faces.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
void correct()
Take over properties from mesh.
dimensionedScalar sin(const dimensionedScalar &ds)
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(psi *p+alphal *rhol0+((alphav *psiv+alphal *psil) - psi) *pSat, rhoMin);# 1 "/home/chef2/andy/OpenFOAM/release/v2306/OpenFOAM-v2306/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=clamp((rho - rholSat)/(rhovSat - rholSat), zero_one{});alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
defineTypeNameAndDebug(combustionModel, 0)
T & last()
Access last element of the list, position [size()-1].
static bool checkFaceDotProduct(const bool report, const scalar orthWarn, const primitiveMesh &, const vectorField &cellCentres, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
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)
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
constexpr scalar radToDeg(const scalar rad) noexcept
Conversion from radians to degrees.
virtual const faceList & faces() const =0
Return faces.
const dimensionedScalar c
Speed of light in a vacuum.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< label > labelList
A List of labels.
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))
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
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)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
vector areaNormal() const
The area normal - with magnitude equal to area of triangle.
static constexpr const zero Zero
Global zero (0)