Commit d01adb74 authored by mattijs's avatar mattijs
Browse files

ENH: Added tet volume check to checkMesh and snappyHexMesh

parent c521e80c
......@@ -327,14 +327,17 @@ meshQualityControls
// Set to 180 to disable.
maxConcave 80;
//- Minimum projected area v.s. actual area. Set to -1 to disable.
minFlatness 0.5;
//- Minimum pyramid volume. Is absolute volume of cell pyramid.
// Set to a sensible fraction of the smallest cell volume expected.
// Set to very negative number (e.g. -1E30) to disable.
minVol 1e-13;
//- Minimum tet volume. Is absolute volume of the tet formed by the
// face-centre decomposition triangle and the cell centre.
// Set to a sensible fraction of the smallest cell volume expected.
// Set to very negative number (e.g. -1E30) to disable.
minTetVol 1e-20;
//- Minimum face area. Set to <0 to disable.
minArea -1;
......
......@@ -159,6 +159,25 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
}
}
{
faceSet faces(mesh, "wrongOrientedTriangleFaces", mesh.nFaces()/100 + 1);
if (mesh.checkFaceTets(true, 0, &faces))
{
noFailedChecks++;
label nFaces = returnReduce(faces.size(), sumOp<label>());
if (nFaces > 0)
{
Info<< " <<Writing " << nFaces
<< " faces with incorrectly orientated triangles to set "
<< faces.name() << endl;
faces.instance() = mesh.pointsInstance();
faces.write();
}
}
}
{
faceSet faces(mesh, "skewFaces", mesh.nFaces()/100 + 1);
if (mesh.checkFaceSkewness(true, &faces))
......
......@@ -556,7 +556,6 @@ public:
//- Check for negative cell volumes
bool checkCellVolumes
(
const bool report = false,
labelHashSet* setPtr = NULL
) const;
......@@ -576,6 +575,14 @@ public:
labelHashSet* setPtr = NULL
) const;
//- Check face-decomposition tet volume
bool checkFaceTets
(
const bool report = false,
const scalar minTetVol = 0,
labelHashSet* setPtr = NULL
) const;
//- Check face skewness
bool checkFaceSkewness
(
......
......@@ -26,6 +26,7 @@ License
#include "primitiveMesh.H"
#include "pyramidPointFaceRef.H"
#include "tetPointRef.H"
#include "ListOps.H"
#include "unitConversion.H"
#include "SortableList.H"
......@@ -401,7 +402,7 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
<< "checking mesh non-orthogonality" << endl;
}
// for all internal faces check theat the d dot S product is positive
// for all internal faces check that the d dot S product is positive
const vectorField& centres = cellCentres();
const vectorField& areas = faceAreas();
......@@ -591,6 +592,114 @@ bool Foam::primitiveMesh::checkFacePyramids
}
bool Foam::primitiveMesh::checkFaceTets
(
const bool report,
const scalar minTetVol,
labelHashSet* setPtr
) const
{
if (debug)
{
Info<< "bool primitiveMesh::checkFaceTets("
<< "const bool, const scalar, labelHashSet*) const: "
<< "checking face orientation" << endl;
}
// check whether face area vector points to the cell with higher label
const vectorField& cc = cellCentres();
const vectorField& fc = faceCentres();
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
const faceList& fcs = faces();
const pointField& p = points();
label nErrorPyrs = 0;
forAll (fcs, faceI)
{
// Create the owner tets - they will have negative volume
const face& f = fcs[faceI];
forAll(f, fp)
{
scalar tetVol = tetPointRef
(
p[f[fp]],
p[f.nextLabel(fp)],
fc[faceI],
cc[own[faceI]]
).mag();
if (tetVol > -minTetVol)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorPyrs++;
break; // no need to check other tets
}
}
if (isInternalFace(faceI))
{
// Create the neighbour tet - it will have positive volume
const face& f = fcs[faceI];
forAll(f, fp)
{
scalar tetVol = tetPointRef
(
p[f[fp]],
p[f.nextLabel(fp)],
fc[faceI],
cc[nei[faceI]]
).mag();
if (tetVol < minTetVol)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorPyrs++;
break;
}
}
}
}
reduce(nErrorPyrs, sumOp<label>());
if (nErrorPyrs > 0)
{
if (debug || report)
{
Info<< " ***Error in face tets: "
<< nErrorPyrs << " faces have incorrectly oriented face"
<< " decomposition triangles." << endl;
}
return true;
}
else
{
if (debug || report)
{
Info<< " Face tets OK." << endl;
}
return false;
}
}
bool Foam::primitiveMesh::checkFaceSkewness
(
const bool report,
......
......@@ -69,6 +69,10 @@ bool Foam::motionSmoother::checkMesh
(
readScalar(dict.lookup("minVol", true))
);
const scalar minTetVol
(
readScalar(dict.lookup("minTetVol", true))
);
const scalar maxConcave
(
readScalar(dict.lookup("maxConcave", true))
......@@ -105,7 +109,6 @@ bool Foam::motionSmoother::checkMesh
(
readScalar(dict.lookup("minDeterminant", true))
);
label nWrongFaces = 0;
Info<< "Checking faces in error :" << endl;
......@@ -158,6 +161,30 @@ bool Foam::motionSmoother::checkMesh
nWrongFaces = nNewWrongFaces;
}
if (minTetVol > -GREAT)
{
polyMeshGeometry::checkFaceTets
(
report,
minTetVol,
mesh,
mesh.cellCentres(),
mesh.faceCentres(),
mesh.points(),
checkFaces,
baffles,
&wrongFaces
);
label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
Info<< " faces with face-decomposition tet volume < "
<< setw(5) << minTetVol << " : "
<< nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
}
if (maxConcave < 180.0-SMALL)
{
polyMeshGeometry::checkFaceAngles
......@@ -417,6 +444,10 @@ bool Foam::motionSmoother::checkMesh
(
readScalar(dict.lookup("minVol", true))
);
const scalar minTetVol
(
readScalar(dict.lookup("minTetVol", true))
);
const scalar maxConcave
(
readScalar(dict.lookup("maxConcave", true))
......@@ -501,6 +532,27 @@ bool Foam::motionSmoother::checkMesh
nWrongFaces = nNewWrongFaces;
}
if (minTetVol > -GREAT)
{
meshGeom.checkFaceTets
(
report,
minTetVol,
meshGeom.mesh().points(),
checkFaces,
baffles,
&wrongFaces
);
label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
Info<< " faces with face-decomposition tet volume < "
<< setw(5) << minTetVol << " : "
<< nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
}
if (maxConcave < 180.0-SMALL)
{
meshGeom.checkFaceAngles
......
......@@ -26,6 +26,7 @@ License
#include "polyMeshGeometry.H"
#include "pyramidPointFaceRef.H"
#include "tetPointRef.H"
#include "syncTools.H"
#include "unitConversion.H"
......@@ -308,6 +309,59 @@ Foam::scalar Foam::polyMeshGeometry::calcSkewness
}
// Create the neighbour pyramid - it will have positive volume
bool Foam::polyMeshGeometry::checkFaceTet
(
const polyMesh& mesh,
const bool report,
const scalar minTetVol,
const pointField& p,
const label faceI,
const point& fc, // face centre
const point& cc, // cell centre
labelHashSet* setPtr
)
{
const face& f = mesh.faces()[faceI];
forAll(f, fp)
{
scalar tetVol = tetPointRef
(
p[f[fp]],
p[f.nextLabel(fp)],
fc,
cc
).mag();
if (tetVol < minTetVol)
{
if (report)
{
Pout<< "bool polyMeshGeometry::checkFaceTets("
<< "const bool, const scalar, const pointField&"
<< ", const pointField&, const labelList&,"
<< " labelHashSet*): "
<< "face " << faceI
<< " has a triangle that points the wrong way."
<< endl
<< "Tet volume: " << tetVol
<< " Face " << faceI
<< endl;
}
if (setPtr)
{
setPtr->insert(faceI);
}
return true;
}
}
return false;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
......@@ -719,6 +773,139 @@ bool Foam::polyMeshGeometry::checkFacePyramids
}
bool Foam::polyMeshGeometry::checkFaceTets
(
const bool report,
const scalar minTetVol,
const polyMesh& mesh,
const vectorField& cellCentres,
const vectorField& faceCentres,
const pointField& p,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet* setPtr
)
{
// check whether face area vector points to the cell with higher label
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
label nErrorPyrs = 0;
forAll(checkFaces, i)
{
label faceI = checkFaces[i];
// Create the owner pyramid - note: exchange cell and face centre
// to get positive volume.
bool tetError = checkFaceTet
(
mesh,
report,
minTetVol,
p,
faceI,
cellCentres[own[faceI]], // face centre
faceCentres[faceI], // cell centre
setPtr
);
if (tetError)
{
nErrorPyrs++;
}
if (mesh.isInternalFace(faceI))
{
// Create the neighbour pyramid - it will have positive volume
bool tetError = checkFaceTet
(
mesh,
report,
minTetVol,
p,
faceI,
faceCentres[faceI], // face centre
cellCentres[nei[faceI]], // cell centre
setPtr
);
if (tetError)
{
nErrorPyrs++;
}
}
}
forAll(baffles, i)
{
label face0 = baffles[i].first();
label face1 = baffles[i].second();
bool tetError = checkFaceTet
(
mesh,
report,
minTetVol,
p,
face0,
cellCentres[own[face0]], // face centre
faceCentres[face0], // cell centre
setPtr
);
if (tetError)
{
nErrorPyrs++;
}
// Create the neighbour pyramid - it will have positive volume
tetError = checkFaceTet
(
mesh,
report,
minTetVol,
p,
face0,
faceCentres[face0], // face centre
cellCentres[own[face1]], // cell centre
setPtr
);
if (tetError)
{
nErrorPyrs++;
}
}
reduce(nErrorPyrs, sumOp<label>());
if (nErrorPyrs > 0)
{
if (report)
{
SeriousErrorIn
(
"polyMeshGeometry::checkFaceTets("
"const bool, const scalar, const pointField&, const pointField&"
", const labelList&, labelHashSet*)"
) << "Error in face pyramids: faces pointing the wrong way!"
<< endl;
}
return true;
}
else
{
if (report)
{
Info<< "Face tets OK.\n" << endl;
}
return false;
}
}
bool Foam::polyMeshGeometry::checkFaceSkewness
(
const bool report,
......@@ -1946,6 +2133,31 @@ bool Foam::polyMeshGeometry::checkFacePyramids
}
bool Foam::polyMeshGeometry::checkFaceTets
(
const bool report,
const scalar minTetVol,
const pointField& p,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet* setPtr
) const
{
return checkFaceTets
(
report,
minTetVol,
mesh_,
cellCentres_,
faceCentres_,
p,
checkFaces,
baffles,
setPtr
);
}
bool Foam::polyMeshGeometry::checkFaceSkewness
(
const bool report,
......
......@@ -109,6 +109,19 @@ class polyMeshGeometry
const point& fc
);
//- Detect&report incorrect face-triangle orientation
static bool checkFaceTet
(
const polyMesh&,
const bool report,
const scalar minTetVol,
const pointField& p,
const label faceI,
const point& fc, // face centre
const point& cc, // cell centre
labelHashSet* setPtr
);
public:
ClassName("polyMeshGeometry");
......@@ -195,6 +208,20 @@ public:
labelHashSet*
);
//- See primitiveMesh
static bool checkFaceTets
(
const bool report,
const scalar minPyrVol,
const polyMesh&,
const vectorField& cellCentres,
const vectorField& faceCentres,
const pointField& p,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet*
);
//- See primitiveMesh
static bool checkFaceSkewness
(
......@@ -321,6 +348,16 @@ public:
labelHashSet* setPtr
) const;
bool checkFaceTets
(
const bool report,
const scalar minTetVol,
const pointField& p,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet* setPtr
) const;
bool checkFaceSkewness
(
const bool report,
......
......@@ -353,14 +353,17 @@ meshQualityControls
// Set to 180 to disable.
maxConcave 80;
//- Minimum projected area v.s. actual area. Set to -1 to disable.
minFlatness 0.5;
//- Minimum pyramid volume. Is absolute volume of cell pyramid.
// Set to a sensible fraction of the smallest cell volume expected.
// Set to very negative number (e.g. -1E30) to disable.
minVol 1e-13;
//- Minimum tet volume. Is absolute volume of the tet formed by the
// face-centre decomposition triangle and the cell centre.
// Set to a sensible fraction of the smallest cell volume expected.
// Set to very negative number (e.g. -1E30) to disable.