41 void Foam::primitiveMeshGeometry::updateFaceCentresAndAreas
51 label facei = changedFaces[i];
54 label nPoints = f.size();
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) =
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) =
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();
235 faceAreas_ = mesh_.faceAreas();
236 faceCentres_ = mesh_.faceCentres();
237 cellCentres_ = mesh_.cellCentres();
238 cellVolumes_ = mesh_.cellVolumes();
249 updateFaceCentresAndAreas(p, changedFaces);
251 updateCellCentresAndVols(affectedCells(changedFaces), changedFaces);
258 const scalar orthWarn,
272 const scalar severeNonorthogonalityThreshold =
::cos(
degToRad(orthWarn));
274 scalar minDDotS = great;
278 label severeNonOrth = 0;
280 label errorNonOrth = 0;
284 label facei = checkFaces[i];
288 vector d = cellCentres[nei[facei]] - cellCentres[own[facei]];
289 const vector&
s = faceAreas[facei];
291 scalar dDotS = (d &
s)/(
mag(d)*
mag(s) + vSmall);
293 if (dDotS < severeNonorthogonalityThreshold)
300 Pout<<
"Severe non-orthogonality for face " << facei
301 <<
" between cells " << own[facei]
302 <<
" and " << nei[facei]
320 <<
"Severe non-orthogonality detected for face " 322 <<
" between cells " << own[facei] <<
" and " 337 if (dDotS < minDDotS)
357 if (report && minDDotS < severeNonorthogonalityThreshold)
359 Info<<
"Number of non-orthogonality errors: " << errorNonOrth
360 <<
". Number of severely non-orthogonal faces: " 361 << severeNonOrth <<
"." <<
endl;
369 Info<<
"Mesh non-orthogonality Max: " 376 if (errorNonOrth > 0)
381 <<
"Error in non-orthogonality detected" <<
endl;
390 Info<<
"Non-orthogonality check OK.\n" <<
endl;
401 const scalar minPyrVol,
415 label nErrorPyrs = 0;
419 label facei = checkFaces[i];
425 cellCentres[own[facei]]
428 if (pyrVol > -minPyrVol)
432 Pout<<
"bool primitiveMeshGeometry::checkFacePyramids(" 433 <<
"const bool, const scalar, const pointField&" 434 <<
", const labelList&, labelHashSet*): " 435 <<
"face " << facei <<
" points the wrong way. " <<
endl 436 <<
"Pyramid volume: " << -pyrVol
437 <<
" Face " << f[facei] <<
" area: " << f[facei].mag(p)
438 <<
" Owner cell: " << own[facei] <<
endl 439 <<
"Owner cell vertex labels: " 440 << mesh.
cells()[own[facei]].labels(f)
459 if (pyrVol < minPyrVol)
463 Pout<<
"bool primitiveMeshGeometry::checkFacePyramids(" 464 <<
"const bool, const scalar, const pointField&" 465 <<
", const labelList&, labelHashSet*): " 466 <<
"face " << facei <<
" points the wrong way. " <<
endl 467 <<
"Pyramid volume: " << -pyrVol
468 <<
" Face " << f[facei] <<
" area: " << f[facei].mag(p)
469 <<
" Neighbour cell: " << nei[facei] <<
endl 470 <<
"Neighbour cell vertex labels: " 471 << mesh.
cells()[nei[facei]].labels(f)
492 <<
"Error in face pyramids: faces pointing the wrong way!" 502 Info<<
"Face pyramids OK.\n" <<
endl;
513 const scalar internalSkew,
514 const scalar boundarySkew,
535 label facei = checkFaces[i];
539 scalar dOwn =
mag(faceCentres[facei] - cellCentres[own[facei]]);
540 scalar dNei =
mag(faceCentres[facei] - cellCentres[nei[facei]]);
542 point faceIntersection =
543 cellCentres[own[facei]]*dNei/(dOwn+dNei)
544 + cellCentres[nei[facei]]*dOwn/(dOwn+dNei);
547 mag(faceCentres[facei] - faceIntersection)
549 mag(cellCentres[nei[facei]]-cellCentres[own[facei]])
556 if (skewness > internalSkew)
560 Pout<<
"Severe skewness for face " << facei
561 <<
" skewness = " << skewness <<
endl;
572 if (skewness > maxSkew)
582 vector faceNormal = faceAreas[facei];
583 faceNormal /=
mag(faceNormal) + vSmall;
585 vector dOwn = faceCentres[facei] - cellCentres[own[facei]];
587 vector dWall = faceNormal*(faceNormal & dOwn);
589 point faceIntersection = cellCentres[own[facei]] + dWall;
592 mag(faceCentres[facei] - faceIntersection)
593 /(2*
mag(dWall) + vSmall);
598 if (skewness > boundarySkew)
602 Pout<<
"Severe skewness for boundary face " << facei
603 <<
" skewness = " << skewness <<
endl;
614 if (skewness > maxSkew)
630 <<
" percent.\nThis may impair the quality of the result." <<
nl 631 << nWarnSkew <<
" highly skew faces detected." 641 Info<<
"Max skewness = " << 100*maxSkew
642 <<
" percent. Face skewness OK.\n" <<
endl;
653 const scalar warnWeight,
667 scalar minWeight = great;
669 label nWarnWeight = 0;
673 label facei = checkFaces[i];
677 const point& fc = faceCentres[facei];
679 scalar dOwn =
mag(faceAreas[facei] & (fc-cellCentres[own[facei]]));
680 scalar dNei =
mag(faceAreas[facei] & (cellCentres[nei[facei]]-fc));
682 scalar weight =
min(dNei,dOwn)/(dNei+dOwn);
684 if (weight < warnWeight)
688 Pout<<
"Small weighting factor for face " << facei
689 <<
" weight = " << weight <<
endl;
700 minWeight =
min(minWeight, weight);
707 if (minWeight < warnWeight)
712 << minWeight <<
'.' <<
nl 713 << nWarnWeight <<
" faces with small weights detected." 723 Info<<
"Min weight = " << minWeight
724 <<
" percent. Weights OK.\n" <<
endl;
747 if (maxDeg < -small || maxDeg > 180+small)
750 <<
"maxDeg should be [0..180] but is now " << maxDeg
758 scalar maxEdgeSin = 0.0;
762 label errorFacei = -1;
766 label facei = checkFaces[i];
768 const face& f = fcs[facei];
770 vector faceNormal = faceAreas[facei];
771 faceNormal /=
mag(faceNormal) + vSmall;
775 scalar magEPrev =
mag(ePrev);
776 ePrev /= magEPrev + vSmall;
784 vector e10(p[f[fp1]] - p[f[fp0]]);
785 scalar magE10 =
mag(e10);
786 e10 /= magE10 + vSmall;
788 if (magEPrev > small && magE10 > small)
790 vector edgeNormal = ePrev ^ e10;
791 scalar magEdgeNormal =
mag(edgeNormal);
793 if (magEdgeNormal < maxSin)
800 edgeNormal /= magEdgeNormal;
802 if ((edgeNormal & faceNormal) < small)
804 if (facei != errorFacei)
816 maxEdgeSin =
max(maxEdgeSin, magEdgeNormal);
831 if (maxEdgeSin > small)
833 scalar maxConcaveDegr =
836 Info<<
"There are " << nConcave
837 <<
" faces with concave angles between consecutive" 838 <<
" edges. Max concave angle = " << maxConcaveDegr
839 <<
" degrees.\n" <<
endl;
843 Info<<
"All angles in faces are convex or less than " << maxDeg
844 <<
" degrees concave.\n" <<
endl;
853 << nConcave <<
" face points with severe concave angle (> " 854 << maxDeg <<
" deg) found.\n" 1002 const scalar minTwist,
1011 if (minTwist < -1-small || minTwist > 1+small)
1014 <<
"minTwist should be [-1..1] but is now " << minTwist
1028 label facei = checkFaces[i];
1030 const face& f = fcs[facei];
1032 scalar magArea =
mag(faceAreas[facei]);
1034 if (f.
size() > 3 && magArea > vSmall)
1036 const vector nf = faceAreas[facei] / magArea;
1038 const point& fc = faceCentres[facei];
1052 scalar magTri =
mag(triArea);
1054 if (magTri > vSmall && ((nf & triArea/magTri) < minTwist))
1074 Info<<
"There are " << nWarped
1075 <<
" faces with cosine of the angle" 1076 <<
" between triangle normal and face normal less than " 1077 << minTwist <<
nl <<
endl;
1081 Info<<
"All faces are flat in that the cosine of the angle" 1082 <<
" between triangle normal and face normal less than " 1083 << minTwist <<
nl <<
endl;
1092 << nWarped <<
" faces with severe warpage " 1093 <<
"(cosine of the angle between triangle normal and " 1094 <<
"face normal < " << minTwist <<
") found.\n" 1110 const scalar minArea,
1117 label nZeroArea = 0;
1121 label facei = checkFaces[i];
1123 if (
mag(faceAreas[facei]) < minArea)
1140 Info<<
"There are " << nZeroArea
1141 <<
" faces with area < " << minArea <<
'.' <<
nl <<
endl;
1145 Info<<
"All faces have area > " << minArea <<
'.' <<
nl <<
endl;
1154 << nZeroArea <<
" faces with area < " << minArea
1171 const scalar warnDet,
1181 scalar minDet = great;
1182 scalar sumDet = 0.0;
1188 const cell& cFaces = cells[affectedCells[i]];
1191 scalar magAreaSum = 0;
1195 label facei = cFaces[cFacei];
1197 scalar magArea =
mag(faceAreas[facei]);
1199 magAreaSum += magArea;
1200 areaSum += faceAreas[facei]*(faceAreas[facei]/magArea);
1203 scalar scaledDet =
det(areaSum/magAreaSum)/0.037037037037037;
1205 minDet =
min(minDet, scaledDet);
1206 sumDet += scaledDet;
1209 if (scaledDet < warnDet)
1216 label facei = cFaces[cFacei];
1233 Info<<
"Cell determinant (1 = uniform cube) : average = " 1234 << sumDet / nSumDet <<
" min = " << minDet <<
endl;
1239 Info<<
"There are " << nWarnDet
1240 <<
" cells with determinant < " << warnDet <<
'.' <<
nl 1245 Info<<
"All faces have determinant > " << warnDet <<
'.' <<
nl 1255 << nWarnDet <<
" cells with determinant < " << warnDet
1272 const scalar orthWarn,
1277 return checkFaceDotProduct
1293 const scalar minPyrVol,
1299 return checkFacePyramids
1315 const scalar internalSkew,
1316 const scalar boundarySkew,
1321 return checkFaceSkewness
1339 const scalar warnWeight,
1344 return checkFaceWeights
1361 const scalar maxDeg,
1367 return checkFaceAngles
1406 const scalar minTwist,
1412 return checkFaceTwist
1429 const scalar minArea,
1434 return checkFaceArea
1449 const scalar warnDet,
1455 return checkCellDeterminant
dimensionedScalar acos(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
A triangle primitive used to calculate face areas and swept volumes.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
A face is a list of labels corresponding to mesh vertices.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Cell-face mesh analysis engine.
Unit conversion functions.
void size(const label)
Override size to be inconsistent with allocated storage.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Ostream & endl(Ostream &os)
Add newline and flush stream.
label nextLabel(const label i) const
Next vertex on face.
const cellList & cells() const
Vector< scalar > vector
A scalar version of the templated Vector.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
bool insert(const Key &key)
Insert a new entry.
T & first()
Return the first element of the list.
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.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
static bool checkFaceArea(const bool report, const scalar minArea, const primitiveMesh &, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
virtual const labelList & faceNeighbour() const =0
Face face-neighbour addressing.
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)
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().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))
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)
List< label > labelList
A List of labels.
primitiveMeshGeometry(const primitiveMesh &)
Construct from mesh.
pyramid< point, const point &, const face & > pyramidPointFaceRef
errorManip< error > abort(error &err)
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
void correct()
Take over properties from mesh.
dimensionedScalar sin(const dimensionedScalar &ds)
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/ubuntu/OpenFOAM-6/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
defineTypeNameAndDebug(combustionModel, 0)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
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.
A cell is defined as a list of faces with extra functionality.
prefixOSstream Pout(cout, "Pout")
virtual const faceList & faces() const =0
Return faces.
const dimensionedScalar c
Speed of light in a vacuum.
virtual const labelList & faceOwner() const =0
Face face-owner addressing.
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
T & last()
Return the last element of the list.
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)