Commit 34e52b8c authored by Franjo's avatar Franjo
Browse files

Merge branch 'task-qualityControls' into development

parents 87b7c279 771a8434
*.dep
*lnInclude*
*linux*
cfMesh.*
cfMesh*
**/lnInclude/*
**/Make/**/*
log
log.*
*~
......@@ -764,6 +764,54 @@ void checkMeshDict::checkRenameBoundary() const
}
}
//Mesh quality criteria specified by the user
void checkMeshDict::checkQualitySettings() const
{
if( meshDict_.found("meshQualitySettings") )
{
const dictionary& qualityDict = meshDict_.subDict("meshQualitySettings");
//- read maximum non-orthogonality defined by the user
if( qualityDict.found("maxNonOrthogonality") )
{
readScalar(qualityDict.lookup("maxNonOrthogonality"));
}
//- read maximum skewness defined by the user
if( qualityDict.found("maxSkewness") )
{
readScalar(qualityDict.lookup("maxSkewness"));
}
//- read minimum volume of the face pyramid defined by the user
if( qualityDict.found("minPyramidVolume") )
{
readScalar(qualityDict.lookup("minPyramidVolume"));
}
//- read face flatness defined by the user
if( qualityDict.found("faceFlatness") )
{
readScalar(qualityDict.lookup("faceFlatness"));
}
//- read minimum tetrahedral part of a cell defined by the user
if( qualityDict.found("minCellPartTetrahedra") )
{
readScalar(qualityDict.lookup("minCellPartTetrahedra"));
}
//- read minimum area of a face defined by the user
if( qualityDict.found("minimumFaceArea") )
{
readScalar(qualityDict.lookup("minimumFaceArea"));
}
}
}
void checkMeshDict::checkEntries() const
{
checkBasicSettings();
......@@ -785,6 +833,8 @@ void checkMeshDict::checkEntries() const
checkBoundaryLayers();
checkRenameBoundary();
checkQualitySettings();
}
void checkMeshDict::updatePatchCellSize
......
......@@ -92,6 +92,9 @@ class checkMeshDict
//- check renameBoundary entry
void checkRenameBoundary() const;
//- check entry for mesh quality
void checkQualitySettings() const;
//- perform all checks
void checkEntries() const;
......
......@@ -289,6 +289,9 @@ namespace help
DynList<bool>& OkPoints
);
//- calculate quality metric of a tetrahedron
inline scalar tetQuality(const tetrahedron<point, point>& tet);
//- check if the vertex is on the positive side of the face plane
inline bool isVertexVisible(const point& p, const plane& pl);
......
......@@ -1612,6 +1612,17 @@ inline bool isFaceConvexAndOk
return valid;
}
inline scalar tetQuality(const tetrahedron<point, point>& tet)
{
return
tet.mag()
/(
8.0/(9.0*sqrt(3.0))
*pow3(min(tet.circumRadius(), GREAT))
+ ROOTVSMALL
);
}
inline point nearestPointOnTheEdge
(
const point& edgePoint0,
......
......@@ -82,6 +82,66 @@ bool checkFaceAreas
const boolList* changedFacePtr = NULL
);
//- Check quality of tetrahedra
bool checkTetQuality
(
const polyMeshGen&,
const bool report = false,
const scalar minTetQuality = VSMALL,
labelHashSet* setPtr = NULL,
const boolList* changedFacePtr = NULL
);
//- Check minimum face twist
bool checkMinTwist
(
const polyMeshGen&,
const bool report = false,
const scalar minTwist = VSMALL,
labelHashSet* setPtr = NULL,
const boolList* changedFacePtr = NULL
);
//- Check the area of internal faces versus boundary faces
bool checkCellDeterminant
(
const polyMeshGen&,
const bool report = false,
const scalar warnDet = 1e-3,
labelHashSet* setPtr = NULL,
const boolList* changedFacePtr = NULL
);
//- Check volume ratio
void checkMinVolRatio
(
const polyMeshGen&,
scalarField&,
const boolList* changedFacePtr = NULL
);
bool checkMinVolRatio
(
const polyMeshGen&,
const bool report = false,
const scalar warnVolRatio = 0.01,
labelHashSet* setPtr = NULL,
const boolList* changedFacePtr = NULL
);
//- Check face triangle twist
bool checkTriangleTwist
(
const polyMeshGen&,
const bool report = false,
const scalar minTwist = VSMALL,
labelHashSet* setPtr = NULL,
const boolList* changedFacePtr = NULL
);
//- Check for negative part tetrahedra
//- Cells are decomposed into tetrahedra formed by
//- the cell centre, face centre and the edge vertices
......@@ -238,6 +298,15 @@ bool checkGeometry(const polyMeshGen&, const bool report = false);
//- Check mesh for correctness. Returns false for no error.
bool checkMesh(const polyMeshGen&, const bool report = false);
//- Read the user defined mesh quality settings
label findBadFacesAdditionalChecks
(
const polyMeshGen& mesh,
const bool report,
labelHashSet& badFaces,
const boolList* activeFacePtr = NULL
);
//- checks for bad faces making the mesh unusable
//- checks for negative pyramids and zero area faces
label findBadFacesRelaxed
......
......@@ -26,7 +26,9 @@ License
#include "polyMeshGenChecks.H"
#include "polyMeshGenAddressing.H"
#include "pyramidPointFaceRef.H"
#include "helperFunctions.H"
#include "tetrahedron.H"
#include "syncTools.H"
# ifdef USE_OMP
#include <omp.h>
......@@ -463,6 +465,993 @@ bool checkFaceAreas
}
}
bool checkTetQuality
(
const polyMeshGen& mesh,
const bool report,
const scalar minTetQuality,
labelHashSet* setPtr,
const boolList* changedFacePtr
)
{
const pointFieldPMG& points = mesh.points();
const faceListPMG& faces = mesh.faces();
const labelList& owner = mesh.owner();
const labelList& neighbour = mesh.neighbour();
const vectorField& cCentres = mesh.addressingData().cellCentres();
label nBadFaces(0);
# ifdef USE_OMP
# pragma omp parallel for schedule(dynamic, 100) reduction(+:nBadFaces)
# endif
forAll(owner, faceI)
{
if( changedFacePtr && !(*changedFacePtr)[faceI] )
continue;
const face& f = faces[faceI];
const label nDecomposed = f.size() - 1;
forAll(f, pI)
{
bool badQualityFace(false);
for(label j=1;j<nDecomposed;++j)
{
const label fpJ = f[(pI+j) % f.size()];
const label nfpJ = f[(pI+j+1) % f.size()];
//- check the tet on the neighbour side
tetrahedron<point, point> tetOwn
(
cCentres[owner[faceI]],
points[f[pI]],
points[fpJ],
points[nfpJ]
);
const scalar tetQualityOwn = help::tetQuality(tetOwn);
if( tetQualityOwn < minTetQuality )
{
++nBadFaces;
if( report )
{
# ifdef USE_OMP
# pragma omp critical(output)
# endif
Pout<< "Face " << faceI
<< " has a triangle that points the wrong way."
<< endl
<< "Tet quality: " << tetQualityOwn
<< " Face " << faceI
<< endl;
}
if( setPtr )
{
# ifdef USE_OMP
# pragma omp critical(insertingBadFaces)
# endif
setPtr->insert(faceI);
}
//- found a problematic face. Do not search further
badQualityFace = true;
continue;
}
if( neighbour[faceI] < 0 )
continue;
//- check the tet on the neighbour side
tetrahedron<point, point> tetNei
(
cCentres[neighbour[faceI]],
points[f[pI]],
points[nfpJ],
points[fpJ]
);
const scalar tetQualityNei = help::tetQuality(tetNei);
if( tetQualityNei < minTetQuality )
{
++nBadFaces;
if( report )
{
# ifdef USE_OMP
# pragma omp critical(output)
# endif
Pout<< "Face " << faceI
<< " has a triangle that points the wrong way."
<< endl
<< "Tet quality: " << tetQualityNei
<< " Face " << faceI
<< endl;
}
if( setPtr )
{
# ifdef USE_OMP
# pragma omp critical(insertingBadFaces)
# endif
setPtr->insert(faceI);
}
//- found a problematic face. Do not search further
badQualityFace = true;
continue;
}
}
if( badQualityFace )
continue;
}
}
reduce(nBadFaces, sumOp<label>());
if( Pstream::parRun() && nBadFaces && setPtr )
{
//- make sure that processor faces are marked on both sides
const PtrList<processorBoundaryPatch>& procBoundaries =
mesh.procBoundaries();
//- send and receive data where needed
forAll(procBoundaries, patchI)
{
const label start = procBoundaries[patchI].patchStart();
const label size = procBoundaries[patchI].patchSize();
labelLongList markedFaces;
for(label faceI=0;faceI<size;++faceI)
{
if( setPtr->found(start+faceI) )
markedFaces.append(faceI);
}
OPstream toOtherProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo(),
markedFaces.byteSize()
);
toOtherProc << markedFaces;
}
forAll(procBoundaries, patchI)
{
labelList receivedData;
IPstream fromOtheProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo()
);
fromOtheProc >> receivedData;
const label start = procBoundaries[patchI].patchStart();
forAll(receivedData, i)
setPtr->insert(start+receivedData[i]);
}
}
if( nBadFaces != 0 )
{
WarningIn
(
"bool checkTetQuality("
"const polyMeshGen&, const bool, const scalar,"
" labelHashSet*, const boolList*)"
) << "Found " << nBadFaces
<< " faces with negative tet decomposition (minTetQuality < "
<< minTetQuality << ")." << endl;
return true;
}
return false;
}
//- check twist of faces
bool checkMinTwist
(
const polyMeshGen& mesh,
const bool report,
const scalar minTwist,
labelHashSet* setPtr,
const boolList* changedFacePtr
)
{
if( (minTwist < (-1.-SMALL)) || (minTwist > (1.+SMALL)) )
{
FatalErrorIn
(
"bool checkMinTwist("
"const polyMeshGen&, const bool, const scalar,"
" labelHashSet*, const boolList*"
) << "minTwist should be [-1..1] but is now " << minTwist
<< abort(FatalError);
}
const faceListPMG& faces = mesh.faces();
label nWarped = 0;
const labelList& own = mesh.owner();
const labelList& nei = mesh.neighbour();
const vectorField& fCentres = mesh.addressingData().faceCentres();
const vectorField& centres = mesh.addressingData().cellCentres();
const pointFieldPMG& points = mesh.points();
const label nInternalFaces = mesh.nInternalFaces();
# ifdef USE_OMP
# pragma omp parallel if ( nInternalFaces > 1000 ) reduction(+:nWarped)
# endif
{
# ifdef USE_OMP
# pragma omp for schedule(guided, 100)
# endif
for(label faceI=0; faceI<nInternalFaces;++faceI)
{
if( changedFacePtr && !(*changedFacePtr)[faceI] )
continue;
const face& f = faces[faceI];
if( f.size() > 3 )
{
vector nf = centres[nei[faceI]] - centres[own[faceI]];
const scalar magNf = mag(nf) + VSMALL;
nf /= magNf;
forAll(f, fpI)
{
const triangle<point, point> triangle
(
points[f[fpI]],
points[f.nextLabel(fpI)],
fCentres[faceI]
);
const vector triArea = triangle.normal();
const scalar magTri = (mag(triArea) + VSMALL);
if( magTri > VSMALL &&
((nf & (triArea/magTri)) < minTwist) )
{
++nWarped;
if( setPtr )
{
# ifdef USE_OMP
# pragma omp critical(badFace)
# endif
setPtr->insert(faceI);
}
break;
}
}
}
}
//- boundary faces
const label start = nInternalFaces;
label end = faces.size();
if( Pstream::parRun() )
{
const PtrList<processorBoundaryPatch>& procBoundaries =
mesh.procBoundaries();
end = procBoundaries[0].patchStart();
}
# ifdef USE_OMP
# pragma omp parallel for schedule(static, 1) reduction(+:nWarped)
# endif
for(label faceI = start; faceI < end;++faceI)
{
const face& f = faces[faceI];
vector nf(vector::zero);
if( f.size() > 3 )
{
nf = fCentres[faceI] - centres[own[faceI]];
nf /= mag(nf) + VSMALL;
}
forAll(f, fpI)
{
const triangle<point, point> triangle
(
points[f[fpI]],
points[f.nextLabel(fpI)],
fCentres[faceI]
);
const vector triArea = triangle.normal();
const scalar magTri = mag(triArea);
if( magTri > VSMALL &&
((nf & triArea/magTri) < minTwist) )
{
nWarped++;
if( setPtr )
{
# ifdef USE_OMP
# pragma omp critical(badFace)
# endif
setPtr->insert(faceI);
}
break;
}
}
}
}
//- if running in parallel
if( Pstream::parRun() )
{
//- check parallel boundaries
const PtrList<processorBoundaryPatch>& procBoundaries =
mesh.procBoundaries();
forAll(procBoundaries, patchI)
{
const label start = procBoundaries[patchI].patchStart();
vectorField cCentres(procBoundaries[patchI].patchSize());
forAll(cCentres, faceI)
cCentres[faceI] = centres[own[start+faceI]];
OPstream toOtherProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo(),
cCentres.byteSize()
);
toOtherProc << cCentres;
}
forAll(procBoundaries, patchI)
{
const label start = procBoundaries[patchI].patchStart();
vectorField otherCentres;
IPstream fromOtherProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo()
);
fromOtherProc >> otherCentres;