zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

primitiveMesh,polyMesh: Initial refactoring of mesh checking

Browse Source
This commit is contained in:
Henry Weller
2023-06-21 19:10:39 +01:00
parent cead8bb02e
commit c40f918c0d
5 changed files with 57 additions and 405 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
src/OpenFOAM/meshes/polyMesh/polyMesh.C
View File

@ -1515,21 +1515,6 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
points_ = newPoints;
bool moveError = false;
if (debug)
{
// Check mesh motion
if (checkMeshMotion(points_, true))
{
moveError = true;
InfoInFunction
<< "Moving the mesh with given points will "
<< "invalidate the mesh." << nl
<< "Mesh motion should not be executed." << endl;
}
}
setPointsInstance(time().name());
tmp<scalarField> sweptVols = primitiveMesh::movePoints
@ -1563,13 +1548,6 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
meshObjects::movePoints<polyMesh>(*this);
meshObjects::movePoints<pointMesh>(*this);
if (debug && moveError)
{
// Write mesh to ease debugging. Note we want to avoid calling
// e.g. fvMesh::write since meshPhi not yet complete.
polyMesh::write();
}
return sweptVols;
}

67
src/OpenFOAM/meshes/polyMesh/polyMesh.H
View File

@ -239,65 +239,6 @@ private:
);
// Geometry checks
//- Check non-orthogonality
bool checkFaceOrthogonality
(
const vectorField& fAreas,
const vectorField& cellCtrs,
const bool report,
const bool detailedReport,
labelHashSet* setPtr
) const;
//- Check face skewness
bool checkFaceSkewness
(
const pointField& points,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs,
const bool report,
const bool detailedReport,
labelHashSet* setPtr
) const;
bool checkEdgeAlignment
(
const pointField& p,
const bool report,
const Vector<label>& directions,
labelHashSet* setPtr
) const;
bool checkCellDeterminant
(
const vectorField& faceAreas,
const bool report,
labelHashSet* setPtr,
const Vector<label>& meshD
) const;
bool checkFaceWeight
(
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs,
const bool report,
const scalar minWeight,
labelHashSet* setPtr
) const;
bool checkVolRatio
(
const scalarField& cellVols,
const bool report,
const scalar minRatio,
labelHashSet* setPtr
) const;
protected:
// Protected Member Data
@ -707,14 +648,6 @@ public:
labelHashSet* setPtr
) const;
//- Check mesh motion for correctness given motion points
virtual bool checkMeshMotion
(
const pointField& newPoints,
const bool report = false,
const bool detailedReport = false
) const;
//- Check for face weights
virtual bool checkFaceWeight
(

257
src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2012-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -32,10 +32,7 @@ License
bool Foam::polyMesh::checkFaceOrthogonality
(
const vectorField& fAreas,
const vectorField& cellCtrs,
const bool report,
const bool detailedReport,
labelHashSet* setPtr
) const
{
@ -44,8 +41,8 @@ bool Foam::polyMesh::checkFaceOrthogonality
InfoInFunction << "Checking mesh non-orthogonality" << endl;
}
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
const vectorField& fAreas = faceAreas();
const vectorField& cellCtrs = cellCentres();
// Calculate orthogonality for all internal and coupled boundary faces
// (1 for uncoupled boundary faces)
@ -59,7 +56,7 @@ bool Foam::polyMesh::checkFaceOrthogonality
// Severe nonorthogonality threshold
const scalar severeNonorthogonalityThreshold =
::cos(degToRad(primitiveMesh::nonOrthThreshold_));
::cos(degToRad(primitiveMeshCheck::nonOrthThreshold));
scalar minDDotS = great;
@ -92,18 +89,6 @@ bool Foam::polyMesh::checkFaceOrthogonality
{
setPtr->insert(facei);
}
if (detailedReport && errorNonOrth == 0)
{
// Non-orthogonality greater than 90 deg
WarningInFunction
<< "Severe non-orthogonality for face "
<< facei
<< " between cells " << own[facei]
<< " and " << nei[facei]
<< ": Angle = "
<< radToDeg(::acos(min(1.0, max(-1.0, ortho[facei]))))
<< " deg." << endl;
}
errorNonOrth++;
}
@ -140,7 +125,7 @@ bool Foam::polyMesh::checkFaceOrthogonality
if (severeNonOrth > 0)
{
Info<< " *Number of severely non-orthogonal (> "
<< primitiveMesh::nonOrthThreshold_ << " degrees) faces: "
<< primitiveMeshCheck::nonOrthThreshold << " degrees) faces: "
<< severeNonOrth << "." << endl;
}
}
@ -169,12 +154,7 @@ bool Foam::polyMesh::checkFaceOrthogonality
bool Foam::polyMesh::checkFaceSkewness
(
const pointField& points,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs,
const bool report,
const bool detailedReport,
labelHashSet* setPtr
) const
{
@ -183,8 +163,10 @@ bool Foam::polyMesh::checkFaceSkewness
InfoInFunction << "Checking face skewness" << endl;
}
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
const pointField& points = this->points();
const vectorField& fCtrs = faceCentres();
const vectorField& fAreas = faceAreas();
const vectorField& cellCtrs = cellCentres();
// Warn if the skew correction vector is more than skewWarning times
// larger than the face area vector
@ -209,31 +191,12 @@ bool Foam::polyMesh::checkFaceSkewness
{
// Check if the skewness vector is greater than the PN vector.
// This does not cause trouble but is a good indication of a poor mesh.
if (skew[facei] > skewThreshold_)
if (skew[facei] > primitiveMeshCheck::skewThreshold)
{
if (setPtr)
{
setPtr->insert(facei);
}
if (detailedReport && nWarnSkew == 0)
{
// Non-orthogonality greater than 90 deg
if (isInternalFace(facei))
{
WarningInFunction
<< "Severe skewness " << skew[facei]
<< " for face " << facei
<< " between cells " << own[facei]
<< " and " << nei[facei];
}
else
{
WarningInFunction
<< "Severe skewness " << skew[facei]
<< " for boundary face " << facei
<< " on cell " << own[facei];
}
}
if (isMasterFace[facei])
{
@ -271,7 +234,6 @@ bool Foam::polyMesh::checkFaceSkewness
bool Foam::polyMesh::checkEdgeAlignment
(
const pointField& p,
const bool report,
const Vector<label>& directions,
labelHashSet* setPtr
@ -289,6 +251,8 @@ bool Foam::polyMesh::checkEdgeAlignment
InfoInFunction << "Checking edge alignment" << endl;
}
const pointField& p = points();
label nDirs = 0;
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
@ -407,12 +371,13 @@ bool Foam::polyMesh::checkEdgeAlignment
bool Foam::polyMesh::checkCellDeterminant
(
const vectorField& faceAreas,
const bool report,
labelHashSet* setPtr,
const Vector<label>& meshD
labelHashSet* setPtr
) const
{
const vectorField& faceAreas = this->faceAreas();
const Vector<label>& meshD = geometricD();
const scalar warnDet = 1e-3;
if (debug)
@ -489,9 +454,6 @@ bool Foam::polyMesh::checkCellDeterminant
bool Foam::polyMesh::checkFaceWeight
(
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs,
const bool report,
const scalar minWeight,
labelHashSet* setPtr
@ -502,6 +464,10 @@ bool Foam::polyMesh::checkFaceWeight
InfoInFunction << "Checking for low face interpolation weights" << endl;
}
const vectorField& fCtrs = faceCentres();
const vectorField& fAreas = faceAreas();
const vectorField& cellCtrs = this->cellCentres();
tmp<scalarField> tfaceWght = polyMeshTools::faceWeights
(
*this,
@ -584,7 +550,6 @@ bool Foam::polyMesh::checkFaceWeight
bool Foam::polyMesh::checkVolRatio
(
const scalarField& cellVols,
const bool report,
const scalar minRatio,
labelHashSet* setPtr
@ -595,6 +560,8 @@ bool Foam::polyMesh::checkVolRatio
InfoInFunction << "Checking for volume ratio < " << minRatio << endl;
}
const scalarField& cellVols = cellVolumes();
tmp<scalarField> tvolRatio = polyMeshTools::volRatio(*this, cellVols);
scalarField& volRatio = tvolRatio.ref();
@ -669,184 +636,4 @@ bool Foam::polyMesh::checkVolRatio
}
bool Foam::polyMesh::checkFaceOrthogonality
(
const bool report,
labelHashSet* setPtr
) const
{
return checkFaceOrthogonality
(
faceAreas(),
cellCentres(),
report,
false, // detailedReport
setPtr
);
}
bool Foam::polyMesh::checkFaceSkewness
(
const bool report,
labelHashSet* setPtr
) const
{
return checkFaceSkewness
(
points(),
faceCentres(),
faceAreas(),
cellCentres(),
report,
false, // detailedReport
setPtr
);
}
bool Foam::polyMesh::checkEdgeAlignment
(
const bool report,
const Vector<label>& directions,
labelHashSet* setPtr
) const
{
return checkEdgeAlignment
(
points(),
report,
directions,
setPtr
);
}
bool Foam::polyMesh::checkCellDeterminant
(
const bool report,
labelHashSet* setPtr
) const
{
return checkCellDeterminant
(
faceAreas(),
report,
setPtr,
geometricD()
);
}
bool Foam::polyMesh::checkFaceWeight
(
const bool report,
const scalar minWeight,
labelHashSet* setPtr
) const
{
return checkFaceWeight
(
faceCentres(),
faceAreas(),
cellCentres(),
report,
minWeight,
setPtr
);
}
bool Foam::polyMesh::checkVolRatio
(
const bool report,
const scalar minRatio,
labelHashSet* setPtr
) const
{
return checkVolRatio(cellVolumes(), report, minRatio, setPtr);
}
bool Foam::polyMesh::checkMeshMotion
(
const pointField& newPoints,
const bool report,
const bool detailedReport
) const
{
if (debug || report)
{
Pout<< "bool polyMesh::checkMeshMotion("
<< "const pointField&, const bool, const bool) const: "
<< "checking mesh motion" << endl;
}
vectorField fCtrs(nFaces());
vectorField fAreas(nFaces());
scalarField magfAreas(nFaces());
makeFaceCentresAndAreas(newPoints, fCtrs, fAreas, magfAreas);
// Check cell volumes and calculate new cell centres
vectorField cellCtrs(nCells());
scalarField cellVols(nCells());
makeCellCentresAndVols(fCtrs, fAreas, cellCtrs, cellVols);
// Check cell volumes
bool error = checkCellVolumes
(
cellVols, // vols
report, // report
detailedReport, // detailedReport
nullptr // setPtr
);
// Check face areas
bool areaError = checkFaceAreas
(
fAreas,
report, // report
detailedReport, // detailedReport,
nullptr // setPtr
);
error = error || areaError;
// Check pyramid volumes
bool pyrVolError = checkFacePyramids
(
newPoints,
cellCtrs,
report, // report,
detailedReport, // detailedReport,
-small, // minPyrVol
nullptr // setPtr
);
error = error || pyrVolError;
// Check face non-orthogonality
bool nonOrthoError = checkFaceOrthogonality
(
fAreas,
cellCtrs,
report, // report
detailedReport, // detailedReport
nullptr // setPtr
);
error = error || nonOrthoError;
if (!error && (debug || report))
{
Pout<< "Mesh motion check OK." << endl;
}
return error;
}
// ************************************************************************* //

56
src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -222,25 +222,6 @@ protected:
// Static Data Members
//- Static data to control mesh checking
//- Cell closedness warning threshold
// set as the fraction of un-closed area to closed area
static scalar closedThreshold_;
//- Aspect ratio warning threshold
static scalar aspectThreshold_;
//- Non-orthogonality warning threshold in deg
static scalar nonOrthThreshold_;
//- Skewness warning threshold
static scalar skewThreshold_;
//- Threshold where faces are considered coplanar
static scalar planarCosAngle_;
// Geometrical calculations
//- Calculate face centres and areas
@ -751,19 +732,6 @@ public:
//- Check mesh for correctness. Returns false for no error.
virtual bool checkMesh(const bool report = false) const;
//- Set the closedness ratio warning threshold
static scalar setClosedThreshold(const scalar);
//- Set the aspect ratio warning threshold
static scalar setAspectThreshold(const scalar);
//- Set the non-orthogonality warning threshold in degrees
static scalar setNonOrthThreshold(const scalar);
//- Set the skewness warning threshold as percentage
// of the face area vector
static scalar setSkewThreshold(const scalar);
// Useful derived info
@ -907,6 +875,28 @@ public:
};
namespace primitiveMeshCheck
{
//- Data to control mesh checking
//- Cell closedness warning threshold
// set as the fraction of un-closed area to closed area
extern scalar closedThreshold;
//- Aspect ratio warning threshold
extern scalar aspectThreshold;
//- Non-orthogonality warning threshold in deg
extern scalar nonOrthThreshold;
//- Skewness warning threshold
extern scalar skewThreshold;
//- Threshold where faces are considered coplanar
extern scalar planarCosAngle;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

60
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -33,11 +33,11 @@ License
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
Foam::scalar Foam::primitiveMesh::closedThreshold_ = 1.0e-6;
Foam::scalar Foam::primitiveMesh::aspectThreshold_ = 1000;
Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70; // deg
Foam::scalar Foam::primitiveMesh::skewThreshold_ = 4;
Foam::scalar Foam::primitiveMesh::planarCosAngle_ = 1.0e-6;
Foam::scalar Foam::primitiveMeshCheck::closedThreshold = 1.0e-6;
Foam::scalar Foam::primitiveMeshCheck::aspectThreshold = 1000;
Foam::scalar Foam::primitiveMeshCheck::nonOrthThreshold = 70; // deg
Foam::scalar Foam::primitiveMeshCheck::skewThreshold = 4;
Foam::scalar Foam::primitiveMeshCheck::planarCosAngle = 1.0e-6;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -75,7 +75,7 @@ bool Foam::primitiveMesh::checkClosedBoundary
vector openness = sumClosed/(sumMagClosedBoundary + vSmall);
if (cmptMax(cmptMag(openness)) > closedThreshold_)
if (cmptMax(cmptMag(openness)) > primitiveMeshCheck::closedThreshold)
{
if (debug || report)
{
@ -167,7 +167,7 @@ bool Foam::primitiveMesh::checkClosedCells
// Check the sums
forAll(openness, celli)
{
if (openness[celli] > closedThreshold_)
if (openness[celli] > primitiveMeshCheck::closedThreshold)
{
if (setPtr)
{
@ -177,7 +177,7 @@ bool Foam::primitiveMesh::checkClosedCells
nOpen++;
}
if (aspectRatio[celli] > aspectThreshold_)
if (aspectRatio[celli] > primitiveMeshCheck::aspectThreshold)
{
if (aspectSetPtr)
{
@ -403,7 +403,7 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
// Severe nonorthogonality threshold
const scalar severeNonorthogonalityThreshold =
::cos(degToRad(nonOrthThreshold_));
::cos(degToRad(primitiveMeshCheck::nonOrthThreshold));
scalar minDDotS = min(ortho);
@ -626,7 +626,7 @@ bool Foam::primitiveMesh::checkFaceSkewness
{
// Check if the skewness vector is greater than the PN vector.
// This does not cause trouble but is a good indication of a poor mesh.
if (skewness[facei] > skewThreshold_)
if (skewness[facei] > primitiveMeshCheck::skewThreshold)
{
if (setPtr)
{
@ -916,7 +916,7 @@ bool Foam::primitiveMesh::checkConcaveCells
pC /= max(mag(pC), vSmall);
if ((pC & fN) > -planarCosAngle_)
if ((pC & fN) > -primitiveMeshCheck::planarCosAngle)
{
// Concave or planar face
@ -1942,40 +1942,4 @@ bool Foam::primitiveMesh::checkMesh(const bool report) const
}
Foam::scalar Foam::primitiveMesh::setClosedThreshold(const scalar val)
{
scalar prev = closedThreshold_;
closedThreshold_ = val;
return prev;
}
Foam::scalar Foam::primitiveMesh::setAspectThreshold(const scalar val)
{
scalar prev = aspectThreshold_;
aspectThreshold_ = val;
return prev;
}
Foam::scalar Foam::primitiveMesh::setNonOrthThreshold(const scalar val)
{
scalar prev = nonOrthThreshold_;
nonOrthThreshold_ = val;
return prev;
}
Foam::scalar Foam::primitiveMesh::setSkewThreshold(const scalar val)
{
scalar prev = skewThreshold_;
skewThreshold_ = val;
return prev;
}
// ************************************************************************* //