Commit 1886d252 authored by Alen Cukrov's avatar Alen Cukrov
Browse files

Implemented all the necessary quality control settings.

parent 16bd8ddf
......@@ -104,6 +104,44 @@ bool checkMinTwist
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
......@@ -260,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
......
......@@ -465,7 +465,6 @@ bool checkFaceAreas
}
}
bool checkTetQuality
(
const polyMeshGen& mesh,
......@@ -583,7 +582,7 @@ bool checkTetQuality
}
//- check twist of faces
bool checkMinTwist
(
const polyMeshGen& mesh,
......@@ -594,7 +593,7 @@ bool checkMinTwist
)
{
if( minTwist < -1.-SMALL || minTwist > 1.+SMALL )
if( (minTwist < (-1.-SMALL)) || (minTwist > (1.+SMALL)) )
{
FatalErrorIn
(
......@@ -627,198 +626,187 @@ bool checkMinTwist
# endif
for(label faceI=0; faceI<nInternalFaces;++faceI)
{
if( changedFacePtr && !(*changedFacePtr)[faceI] )
continue;
const face& f = faces[faceI];
if( f.size() > 3 )
{
vector nf(vector::zero);
nf = centres[nei[faceI]] - centres[own[faceI]];
nf /= mag(nf) + VSMALL;
vector nf = centres[nei[faceI]] - centres[own[faceI]];
const scalar magNf = mag(nf) + VSMALL;
nf /= magNf;
if( magSqr(nf) > VSMALL )
forAll(f, fpI)
{
forAll(f, fpI)
{
const triangle<point, point> triangle
(
points[f[fpI]],
points[f.nextLabel(fpI)],
fCentres[faceI]
);
const vector triArea = triangle.normal();
const triangle<point, point> triangle
(
points[f[fpI]],
points[f.nextLabel(fpI)],
fCentres[faceI]
);
const scalar magTri = mag(triArea);
const vector triArea = triangle.normal();
if( magTri > VSMALL &&
((nf & triArea/magTri) < minTwist) )
{
nWarped++;
const scalar magTri = (mag(triArea) + VSMALL);
if( setPtr )
{
# ifdef USE_OMP
# pragma omp critical
# endif
setPtr->insert(faceI);
}
if( magTri > VSMALL &&
((nf & (triArea/magTri)) < minTwist) )
{
++nWarped;
break;
if( setPtr )
{
# ifdef USE_OMP
# pragma omp critical
# endif
setPtr->insert(faceI);
}
break;
}
}
}
}
//If running in parallel
//- boundary faces
const label start = nInternalFaces;
label end = faces.size();
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]];
mesh.procBoundaries();
end = procBoundaries[0].patchStart();
}
OPstream toOtherProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo(),
cCentres.byteSize()
);
# 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];
toOtherProc << cCentres;
}
vector nf(vector::zero);
forAll(procBoundaries, patchI)
if( f.size() > 3 )
{
const label start = procBoundaries[patchI].patchStart();
vectorField otherCentres;
nf = fCentres[faceI] - centres[own[faceI]];
nf /= mag(nf) + VSMALL;
}
IPstream fromOtherProc
forAll(f, fpI)
{
const triangle<point, point> triangle
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo()
points[f[fpI]],
points[f.nextLabel(fpI)],
fCentres[faceI]
);
fromOtherProc >> otherCentres;
forAll(otherCentres, fI)
{
const label faceI = start + fI;
if( changedFacePtr && !(*changedFacePtr)[faceI] )
continue;
const point& cOwn = centres[own[faceI]];
const point& cNei = otherCentres[fI];
const vector triArea = triangle.normal();
const face& f = faces[faceI];
const scalar magTri = mag(triArea);
vector nf(vector::zero);
if( magTri > VSMALL &&
((nf & triArea/magTri) < minTwist) )
{
nWarped++;
if( f.size() > 3 )
if( setPtr )
{
nf = cNei - cOwn;
nf /= mag(nf) + VSMALL;
# ifdef USE_OMP
# pragma omp critical
# endif
setPtr->insert(faceI);
}
break;
}
}
}
}
if( magSqr(nf) > VSMALL )
{
forAll(f, fpI)
{
const triangle<point, point> triangle
(
points[f[fpI]],
points[f.nextLabel(fpI)],
fCentres[faceI]
);
//- if running in parallel
if( Pstream::parRun() )
{
//- check parallel boundaries
const PtrList<processorBoundaryPatch>& procBoundaries =
mesh.procBoundaries();
const vector triArea = triangle.normal();
forAll(procBoundaries, patchI)
{
const label start = procBoundaries[patchI].patchStart();
const scalar magTri = mag(triArea);
vectorField cCentres(procBoundaries[patchI].patchSize());
if( magTri > VSMALL &&
((nf & triArea/magTri) < minTwist) )
{
nWarped++;
forAll(cCentres, faceI)
cCentres[faceI] = centres[own[start+faceI]];
if( setPtr )
{
setPtr->insert(faceI);
}
OPstream toOtherProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo(),
cCentres.byteSize()
);
break;
}
}
}
}
}
toOtherProc << cCentres;
}
//- boundary faces
const label start = nInternalFaces;
label end = faces.size();
if( Pstream::parRun() )
forAll(procBoundaries, patchI)
{
const PtrList<processorBoundaryPatch>& procBoundaries =
mesh.procBoundaries();
end = procBoundaries[0].patchStart();
}
const label start = procBoundaries[patchI].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];
vectorField otherCentres;
vector nf(vector::zero);
IPstream fromOtherProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo()
);
if( f.size() > 3 )
fromOtherProc >> otherCentres;
forAll(otherCentres, fI)
{
nf = fCentres[faceI] - centres[own[faceI]];
nf /= mag(nf) + VSMALL;
}
const label faceI = start + fI;
if( magSqr(nf) > VSMALL )
{
if( changedFacePtr && !(*changedFacePtr)[faceI] )
continue;
forAll(f, fpI)
{
const triangle<point, point> triangle
(
points[f[fpI]],
points[f.nextLabel(fpI)],
fCentres[faceI]
);
const point& cOwn = centres[own[faceI]];
const point& cNei = otherCentres[fI];
const vector triArea = triangle.normal();
const face& f = faces[faceI];
const scalar magTri = mag(triArea);
if( f.size() > 3 )
{
vector nf = cNei - cOwn;
nf /= (mag(nf) + VSMALL);
if( magTri > VSMALL &&
((nf & triArea/magTri) < minTwist) )
forAll(f, fpI)
{
nWarped++;
const triangle<point, point> triangle
(
points[f[fpI]],
points[f.nextLabel(fpI)],
fCentres[faceI]
);
if( setPtr )
const vector triArea = triangle.normal();
const scalar magTri = (mag(triArea) + VSMALL);
if( magTri > VSMALL &&
((nf & triArea/magTri) < minTwist) )
{
# ifdef USE_OMP
# pragma omp critical
# endif
setPtr->insert(faceI);
++nWarped;
if( setPtr )
{
setPtr->insert(faceI);
}
break;
}
break;
}
}
}
......@@ -868,402 +856,919 @@ bool checkMinTwist
}
bool checkCellPartTetrahedra
bool checkCellDeterminant
(
const polyMeshGen& mesh,
const bool report,
const scalar minPartTet,
const scalar warnDet,
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 cellListPMG& cells = mesh.cells();
const labelList& own = mesh.owner();
const labelList& nei = mesh.neighbour();
const label nInternalFaces = mesh.nInternalFaces();
const vectorField& fCentres = mesh.addressingData().faceCentres();
const vectorField& cCentres = mesh.addressingData().cellCentres();
scalar minDet = VGREAT;
scalar sumDet = 0.0;
label nSumDet = 0;
label nWarnDet = 0;
boolList affectedCells(cells.size(), false);
if( changedFacePtr )
{
const boolList& changedFaces = *changedFacePtr;
forAll(changedFaces, fI)
{
if( changedFaces[fI] )
{
affectedCells[own[fI]] = true;
label nNegVolCells = 0;
if( fI >= 0 && fI < nInternalFaces )
{
affectedCells[nei[fI]] = true;
}
}
}
}
# ifdef USE_OMP
# pragma omp parallel for if( owner.size() > 100 ) \
schedule(guided) reduction(+ : nNegVolCells)
# pragma omp parallel
# endif
forAll(owner, faceI)
{
if( changedFacePtr && !(*changedFacePtr)[faceI] )
continue;
const face& f = faces[faceI];
bool badFace(false);
scalar localMinDet(VGREAT);
forAll(f, eI)
# ifdef USE_OMP
# pragma omp for schedule(guided, 50) \
reduction(+: nSumDet, sumDet, nWarnDet)
# endif
forAll(cells, cI)
{
if( affectedCells[cI] )
{
const cell& cFaces = cells[cI];
const vectorField& areas = mesh.addressingData().faceAreas();
const tetrahedron<point, point> tetOwn
(
fCentres[faceI],
points[f.nextLabel(eI)],
points[f[eI]],
cCentres[owner[faceI]]
);
tensor areaSum(tensor::zero);
scalar magAreaSum = 0.0;
if( tetOwn.mag() < minPartTet )
{
if( report )
forAll(cFaces, cFaceI)
{
# ifdef USE_OMP
# pragma omp critical
# endif
Pout<< "Zero or negative cell volume detected for cell "
<< owner[faceI] << "." << endl;
}
label faceI = cFaces[cFaceI];
badFace = true;
}
scalar magArea = mag(areas[faceI]) + VSMALL;
if( neighbour[faceI] < 0 )
continue;
magAreaSum += magArea;
areaSum += areas[faceI]*(areas[faceI]/magArea);
}
const tetrahedron<point, point> tetNei
(
fCentres[faceI],
points[f[eI]],
points[f.nextLabel(eI)],
cCentres[neighbour[faceI]]
);
const scalar scaledDet =
det(areaSum/magAreaSum)/0.037037037037037;
if( tetNei.mag() < minPartTet )
{
if( report )
localMinDet = min(localMinDet, scaledDet);
sumDet += scaledDet;
++nSumDet;
if( scaledDet < warnDet )
{
# ifdef USE_OMP
# pragma omp critical
# endif
Pout<< "Zero or negative cell volume detected for cell "
<< neighbour[faceI] << "." << endl;
}
if( setPtr )
{
//Insert all faces of the cell
forAll(cFaces, cFaceI)
{
label faceI = cFaces[cFaceI];
badFace = true;
# ifdef USE_OMP
# pragma omp critical
# endif
setPtr->insert(faceI);
}
}
++nWarnDet;
}
}
}
if( badFace )
{
if( setPtr )
{
# ifdef USE_OMP
# pragma omp critical
# endif