Commit ceecb845 authored by mattijs's avatar mattijs
Browse files

ENH: checkMesh: moved some checking functionality to polyMesh

parent 5bbbcdd5
......@@ -803,7 +803,7 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
if (allGeometry)
{
cellSet cells(mesh, "underdeterminedCells", mesh.nCells()/100);
if (mesh.checkCellDeterminant(true, &cells, mesh.geometricD()))
if (mesh.checkCellDeterminant(true, &cells))
{
noFailedChecks++;
......
......@@ -425,6 +425,10 @@ $(polyMesh)/polyMeshInitMesh.C
$(polyMesh)/polyMeshClear.C
$(polyMesh)/polyMeshUpdate.C
polyMeshCheck = $(polyMesh)/polyMeshCheck
$(polyMeshCheck)/polyMeshCheck.C
$(polyMeshCheck)/polyMeshTools.C
primitiveMesh = meshes/primitiveMesh
$(primitiveMesh)/primitiveMesh.C
$(primitiveMesh)/primitiveMeshCellCells.C
......@@ -445,9 +449,9 @@ $(primitiveMesh)/primitiveMeshCalcCellShapes.C
primitiveMeshCheck = $(primitiveMesh)/primitiveMeshCheck
$(primitiveMeshCheck)/primitiveMeshCheck.C
$(primitiveMeshCheck)/primitiveMeshCheckMotion.C
$(primitiveMeshCheck)/primitiveMeshCheckPointNearness.C
$(primitiveMeshCheck)/primitiveMeshCheckEdgeLength.C
$(primitiveMeshCheck)/primitiveMeshTools.C
primitivePatch = $(primitiveMesh)/primitivePatch
$(primitivePatch)/patchZones.C
......
......@@ -1129,7 +1129,7 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
if (debug)
{
// Check mesh motion
if (primitiveMesh::checkMeshMotion(points_, true))
if (checkMeshMotion(points_, true))
{
Info<< "tmp<scalarField> polyMesh::movePoints"
<< "(const pointField&) : "
......
......@@ -34,6 +34,7 @@ SourceFiles
polyMeshFromShapeMesh.C
polyMeshIO.C
polyMeshUpdate.C
polyMeshCheck.C
\*---------------------------------------------------------------------------*/
......@@ -224,7 +225,45 @@ private:
cellList& cells
);
// 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;
public:
......@@ -535,6 +574,48 @@ public:
void removeFiles() const;
// Geometric checks. Selectively override primitiveMesh functionality.
//- Check boundary for closedness
virtual bool checkClosedBoundary(const bool report = false) const;
//- Check non-orthogonality
virtual bool checkFaceOrthogonality
(
const bool report = false,
labelHashSet* setPtr = NULL
) const;
//- Check face skewness
virtual bool checkFaceSkewness
(
const bool report = false,
labelHashSet* setPtr = NULL
) const;
//- Check edge alignment for 1D/2D cases
virtual bool checkEdgeAlignment
(
const bool report,
const Vector<label>& directions,
labelHashSet* setPtr
) const;
virtual bool checkCellDeterminant
(
const bool report,
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;
// Helper functions
//- Find the cell, tetFaceI and tetPtI for the given position
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "polyMesh.H"
#include "polyMeshTools.H"
#include "unitConversion.H"
#include "syncTools.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::polyMesh::checkFaceOrthogonality
(
const vectorField& fAreas,
const vectorField& cellCtrs,
const bool report,
const bool detailedReport,
labelHashSet* setPtr
) const
{
if (debug)
{
Info<< "bool polyMesh::checkFaceOrthogonality("
<< "const bool, labelHashSet*) const: "
<< "checking mesh non-orthogonality" << endl;
}
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
// Calculate orthogonality for all internal and coupled boundary faces
// (1 for uncoupled boundary faces)
tmp<scalarField> tortho = polyMeshTools::faceOrthogonality
(
*this,
fAreas,
cellCtrs
);
const scalarField& ortho = tortho();
// Severe nonorthogonality threshold
const scalar severeNonorthogonalityThreshold =
::cos(degToRad(primitiveMesh::nonOrthThreshold_));
scalar minDDotS = GREAT;
scalar sumDDotS = 0.0;
label nSummed = 0;
label severeNonOrth = 0;
label errorNonOrth = 0;
// Statistics only for internal and masters of coupled faces
PackedBoolList isMasterFace(syncTools::getInternalOrMasterFaces(*this));
forAll(ortho, faceI)
{
if (ortho[faceI] < severeNonorthogonalityThreshold)
{
if (ortho[faceI] > SMALL)
{
if (setPtr)
{
setPtr->insert(faceI);
}
severeNonOrth++;
}
else
{
// Error : non-ortho too large
if (setPtr)
{
setPtr->insert(faceI);
}
if (detailedReport && errorNonOrth == 0)
{
// Non-orthogonality greater than 90 deg
WarningIn
(
"polyMesh::checkFaceOrthogonality"
"(const pointField&, const bool) const"
) << "Severe non-orthogonality for face "
<< faceI
<< " between cells " << own[faceI]
<< " and " << nei[faceI]
<< ": Angle = " << radToDeg(::acos(ortho[faceI]))
<< " deg." << endl;
}
errorNonOrth++;
}
}
if (isMasterFace[faceI])
{
minDDotS = min(minDDotS, ortho[faceI]);
sumDDotS += ortho[faceI];
nSummed++;
}
}
reduce(minDDotS, minOp<scalar>());
reduce(sumDDotS, sumOp<scalar>());
reduce(nSummed, sumOp<label>());
reduce(severeNonOrth, sumOp<label>());
reduce(errorNonOrth, sumOp<label>());
if (debug || report)
{
if (nSummed > 0)
{
if (debug || report)
{
Info<< " Mesh non-orthogonality Max: "
<< radToDeg(::acos(minDDotS))
<< " average: " << radToDeg(::acos(sumDDotS/nSummed))
<< endl;
}
}
if (severeNonOrth > 0)
{
Info<< " *Number of severely non-orthogonal faces: "
<< severeNonOrth << "." << endl;
}
}
if (errorNonOrth > 0)
{
if (debug || report)
{
Info<< " ***Number of non-orthogonality errors: "
<< errorNonOrth << "." << endl;
}
return true;
}
else
{
if (debug || report)
{
Info<< " Non-orthogonality check OK." << endl;
}
return false;
}
}
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
{
if (debug)
{
Info<< "bool polyMesh::checkFaceSkewnesss("
<< "const bool, labelHashSet*) const: "
<< "checking face skewness" << endl;
}
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
// Warn if the skew correction vector is more than skewWarning times
// larger than the face area vector
tmp<scalarField> tskew = polyMeshTools::faceSkewness
(
*this,
points,
fCtrs,
fAreas,
cellCtrs
);
const scalarField& skew = tskew();
scalar maxSkew = max(skew);
label nWarnSkew = 0;
// Statistics only for all faces except slave coupled faces
PackedBoolList isMasterFace(syncTools::getMasterFaces(*this));
forAll(skew, faceI)
{
// 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 (setPtr)
{
setPtr->insert(faceI);
}
if (detailedReport && nWarnSkew == 0)
{
// Non-orthogonality greater than 90 deg
if (isInternalFace(faceI))
{
WarningIn
(
"polyMesh::checkFaceSkewnesss"
"(const pointField&, const bool) const"
) << "Severe skewness " << skew[faceI]
<< " for face " << faceI
<< " between cells " << own[faceI]
<< " and " << nei[faceI];
}
else
{
WarningIn
(
"polyMesh::checkFaceSkewnesss"
"(const pointField&, const bool) const"
) << "Severe skewness " << skew[faceI]
<< " for boundary face " << faceI
<< " on cell " << own[faceI];
}
}
if (isMasterFace[faceI])
{
nWarnSkew++;
}
}
}
reduce(maxSkew, maxOp<scalar>());
reduce(nWarnSkew, sumOp<label>());
if (nWarnSkew > 0)
{
if (debug || report)
{
Info<< " ***Max skewness = " << maxSkew
<< ", " << nWarnSkew << " highly skew faces detected"
" which may impair the quality of the results"
<< endl;
}
return true;
}
else
{
if (debug || report)
{
Info<< " Max skewness = " << maxSkew << " OK." << endl;
}
return false;
}
}
// Check 1D/2Dness of edges. Gets passed the non-empty directions and
// checks all edges in the mesh whether they:
// - have no component in a non-empty direction or
// - are only in a singe non-empty direction.
// Empty direction info is passed in as a vector of labels (synchronised)
// which are 1 if the direction is non-empty, 0 if it is.
bool Foam::polyMesh::checkEdgeAlignment
(
const pointField& p,
const bool report,
const Vector<label>& directions,
labelHashSet* setPtr
) const
{
if (debug)
{
Info<< "bool polyMesh::checkEdgeAlignment("
<< "const bool, const Vector<label>&, labelHashSet*) const: "
<< "checking edge alignment" << endl;
}
label nDirs = 0;
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
if (directions[cmpt] == 1)
{
nDirs++;
}
else if (directions[cmpt] != 0)
{
FatalErrorIn
(
"polyMesh::checkEdgeAlignment"
"(const bool, const Vector<label>&, labelHashSet*)"
) << "directions should contain 0 or 1 but is now " << directions
<< exit(FatalError);
}
}
if (nDirs == vector::nComponents)
{
return false;
}
const faceList& fcs = faces();
EdgeMap<label> edgesInError;
forAll(fcs, faceI)
{
const face& f = fcs[faceI];
forAll(f, fp)
{
label p0 = f[fp];
label p1 = f.nextLabel(fp);
if (p0 < p1)
{
vector d(p[p1]-p[p0]);
scalar magD = mag(d);
if (magD > ROOTVSMALL)
{
d /= magD;
// Check how many empty directions are used by the edge.
label nEmptyDirs = 0;
label nNonEmptyDirs = 0;
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
if (mag(d[cmpt]) > 1e-6)
{
if (directions[cmpt] == 0)
{
nEmptyDirs++;
}
else
{
nNonEmptyDirs++;
}
}
}
if (nEmptyDirs == 0)
{
// Purely in ok directions.
}
else if (nEmptyDirs == 1)
{
// Ok if purely in empty directions.
if (nNonEmptyDirs > 0)
{
edgesInError.insert(edge(p0, p1), faceI);
}
}
else if (nEmptyDirs > 1)
{
// Always an error
edgesInError.insert(edge(p0, p1), faceI);
}
}
}
}
}
label nErrorEdges = returnReduce(edgesInError.size(), sumOp<label>());
if (nErrorEdges > 0)
{
if (debug || report)
{
Info<< " ***Number of edges not aligned with or perpendicular to "
<< "non-empty directions: " << nErrorEdges << endl;
}
if (setPtr)
{
setPtr->resize(2*edgesInError.size());
forAllConstIter(EdgeMap<label>, edgesInError, iter)
{
setPtr->insert(iter.key()[0]);
setPtr->insert(iter.key()[1]);
}
}
return true;
}
else
{
if (debug || report)
{
Info<< " All edges aligned with or perpendicular to "
<< "non-empty directions." << endl;
}
return false;
}
}
bool Foam::polyMesh::checkCellDeterminant
(
const vectorField& faceAreas,
const bool report,
labelHashSet* setPtr,
const Vector<label>& meshD
) const
{
if (debug)
{
Info<< "bool polyMesh::checkCellDeterminant(const bool"
<< ", labelHashSet*) const: "
<< "checking for under-determined cells" << endl;
}
// Determine number of dimensions and (for 2D) missing dimension
label nDims = 0;
label twoD = -1;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (meshD[dir] == 1)
{
nDims++;
}
else
{
twoD = dir;
}
}
tmp<scalarField> tcellDeterminant = primitiveMeshTools::cellDeterminant
(
*this,
meshD,
faceAreas,
syncTools::getInternalOrCoupledFaces(*this)
);
scalarField& cellDeterminant = tcellDeterminant();
label nErrorCells = 0;
scalar minDet = min(cellDeterminant);