Commit 771a8434 authored by Franjo's avatar Franjo
Browse files

Updated mesh checks

parent 264432b6
...@@ -479,7 +479,6 @@ bool checkTetQuality ...@@ -479,7 +479,6 @@ bool checkTetQuality
const labelList& owner = mesh.owner(); const labelList& owner = mesh.owner();
const labelList& neighbour = mesh.neighbour(); const labelList& neighbour = mesh.neighbour();
const vectorField& fCentres = mesh.addressingData().faceCentres();
const vectorField& cCentres = mesh.addressingData().cellCentres(); const vectorField& cCentres = mesh.addressingData().cellCentres();
label nBadFaces(0); label nBadFaces(0);
...@@ -494,15 +493,24 @@ bool checkTetQuality ...@@ -494,15 +493,24 @@ bool checkTetQuality
const face& f = faces[faceI]; const face& f = faces[faceI];
forAll(f, eI) 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 //- check the tet on the neighbour side
tetrahedron<point, point> tetOwn tetrahedron<point, point> tetOwn
( (
fCentres[faceI], cCentres[owner[faceI]],
points[f[eI]], points[f[pI]],
points[f.nextLabel(eI)], points[fpJ],
cCentres[owner[faceI]] points[nfpJ]
); );
const scalar tetQualityOwn = help::tetQuality(tetOwn); const scalar tetQualityOwn = help::tetQuality(tetOwn);
...@@ -514,7 +522,7 @@ bool checkTetQuality ...@@ -514,7 +522,7 @@ bool checkTetQuality
if( report ) if( report )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(output)
# endif # endif
Pout<< "Face " << faceI Pout<< "Face " << faceI
<< " has a triangle that points the wrong way." << " has a triangle that points the wrong way."
...@@ -531,6 +539,10 @@ bool checkTetQuality ...@@ -531,6 +539,10 @@ bool checkTetQuality
# endif # endif
setPtr->insert(faceI); setPtr->insert(faceI);
} }
//- found a problematic face. Do not search further
badQualityFace = true;
continue;
} }
if( neighbour[faceI] < 0 ) if( neighbour[faceI] < 0 )
...@@ -539,10 +551,10 @@ bool checkTetQuality ...@@ -539,10 +551,10 @@ bool checkTetQuality
//- check the tet on the neighbour side //- check the tet on the neighbour side
tetrahedron<point, point> tetNei tetrahedron<point, point> tetNei
( (
fCentres[faceI], cCentres[neighbour[faceI]],
points[f.nextLabel(eI)], points[f[pI]],
points[f[eI]], points[nfpJ],
cCentres[neighbour[faceI]] points[fpJ]
); );
const scalar tetQualityNei = help::tetQuality(tetNei); const scalar tetQualityNei = help::tetQuality(tetNei);
...@@ -554,7 +566,7 @@ bool checkTetQuality ...@@ -554,7 +566,7 @@ bool checkTetQuality
if( report ) if( report )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(output)
# endif # endif
Pout<< "Face " << faceI Pout<< "Face " << faceI
<< " has a triangle that points the wrong way." << " has a triangle that points the wrong way."
...@@ -571,12 +583,78 @@ bool checkTetQuality ...@@ -571,12 +583,78 @@ bool checkTetQuality
# endif # endif
setPtr->insert(faceI); 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 ) 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 true;
}
return false; return false;
} }
...@@ -622,7 +700,7 @@ bool checkMinTwist ...@@ -622,7 +700,7 @@ bool checkMinTwist
# endif # endif
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp for schedule(static, 1) # pragma omp for schedule(guided, 100)
# endif # endif
for(label faceI=0; faceI<nInternalFaces;++faceI) for(label faceI=0; faceI<nInternalFaces;++faceI)
{ {
...@@ -658,7 +736,7 @@ bool checkMinTwist ...@@ -658,7 +736,7 @@ bool checkMinTwist
if( setPtr ) if( setPtr )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(badFace)
# endif # endif
setPtr->insert(faceI); setPtr->insert(faceI);
} }
...@@ -715,7 +793,7 @@ bool checkMinTwist ...@@ -715,7 +793,7 @@ bool checkMinTwist
if( setPtr ) if( setPtr )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(badFace)
# endif # endif
setPtr->insert(faceI); setPtr->insert(faceI);
} }
...@@ -906,17 +984,17 @@ bool checkCellDeterminant ...@@ -906,17 +984,17 @@ bool checkCellDeterminant
{ {
if( affectedCells[cI] ) if( affectedCells[cI] )
{ {
const cell& cFaces = cells[cI]; const cell& c = cells[cI];
const vectorField& areas = mesh.addressingData().faceAreas(); const vectorField& areas = mesh.addressingData().faceAreas();
tensor areaSum(tensor::zero); tensor areaSum(tensor::zero);
scalar magAreaSum = 0.0; scalar magAreaSum = 0.0;
forAll(cFaces, cFaceI) forAll(c, fI)
{ {
label faceI = cFaces[cFaceI]; const label faceI = c[fI];
scalar magArea = mag(areas[faceI]) + VSMALL; const scalar magArea = mag(areas[faceI]) + VSMALL;
magAreaSum += magArea; magAreaSum += magArea;
areaSum += areas[faceI]*(areas[faceI]/magArea); areaSum += areas[faceI]*(areas[faceI]/magArea);
...@@ -934,12 +1012,12 @@ bool checkCellDeterminant ...@@ -934,12 +1012,12 @@ bool checkCellDeterminant
if( setPtr ) if( setPtr )
{ {
//Insert all faces of the cell //Insert all faces of the cell
forAll(cFaces, cFaceI) forAll(c, fI)
{ {
label faceI = cFaces[cFaceI]; const label faceI = c[fI];
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(badFace)
# endif # endif
setPtr->insert(faceI); setPtr->insert(faceI);
} }
...@@ -950,7 +1028,7 @@ bool checkCellDeterminant ...@@ -950,7 +1028,7 @@ bool checkCellDeterminant
} }
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(minDet)
# endif # endif
minDet = min(minDet, localMinDet); minDet = min(minDet, localMinDet);
} }
...@@ -1017,13 +1095,14 @@ void checkMinVolRatio ...@@ -1017,13 +1095,14 @@ void checkMinVolRatio
const scalarField& vols = mesh.addressingData().cellVolumes(); const scalarField& vols = mesh.addressingData().cellVolumes();
volRatio.setSize(own.size()); volRatio.setSize(own.size());
volRatio = 1.0;
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp parallel for schedule (dynamic, 100) # pragma omp parallel for schedule (dynamic, 100)
# endif # endif
for(label faceI=0; faceI<nInternalFaces;++faceI) for(label faceI=0; faceI<nInternalFaces;++faceI)
{ {
volRatio[faceI] = 1.0;
if( changedFacePtr && !changedFacePtr->operator[](faceI) ) if( changedFacePtr && !changedFacePtr->operator[](faceI) )
continue; continue;
...@@ -1142,7 +1221,7 @@ bool checkMinVolRatio ...@@ -1142,7 +1221,7 @@ bool checkMinVolRatio
if( setPtr ) if( setPtr )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(badFace)
# endif # endif
setPtr->insert(faceI); setPtr->insert(faceI);
} }
...@@ -1156,7 +1235,7 @@ bool checkMinVolRatio ...@@ -1156,7 +1235,7 @@ bool checkMinVolRatio
} }
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(minVolRatio)
# endif # endif
{ {
maxVolRatio = Foam::max(maxVolRatio , localMaxVolRatio); maxVolRatio = Foam::max(maxVolRatio , localMaxVolRatio);
...@@ -1321,7 +1400,7 @@ bool checkTriangleTwist ...@@ -1321,7 +1400,7 @@ bool checkTriangleTwist
if( setPtr ) if( setPtr )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(badFace)
# endif # endif
setPtr->insert(faceI); setPtr->insert(faceI);
} }
...@@ -1420,7 +1499,7 @@ bool checkCellPartTetrahedra ...@@ -1420,7 +1499,7 @@ bool checkCellPartTetrahedra
if( report ) if( report )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(report)
# endif # endif
Pout<< "Zero or negative cell volume detected for cell " Pout<< "Zero or negative cell volume detected for cell "
<< owner[faceI] << "." << endl; << owner[faceI] << "." << endl;
...@@ -1445,7 +1524,7 @@ bool checkCellPartTetrahedra ...@@ -1445,7 +1524,7 @@ bool checkCellPartTetrahedra
if( report ) if( report )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(report)
# endif # endif
Pout<< "Zero or negative cell volume detected for cell " Pout<< "Zero or negative cell volume detected for cell "
<< neighbour[faceI] << "." << endl; << neighbour[faceI] << "." << endl;
...@@ -1460,7 +1539,7 @@ bool checkCellPartTetrahedra ...@@ -1460,7 +1539,7 @@ bool checkCellPartTetrahedra
if( setPtr ) if( setPtr )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(badFace)
# endif # endif
setPtr->insert(faceI); setPtr->insert(faceI);
} }
...@@ -1687,7 +1766,7 @@ bool checkFaceDotProduct ...@@ -1687,7 +1766,7 @@ bool checkFaceDotProduct
{ {
// Severe non-orthogonality but mesh still OK // Severe non-orthogonality but mesh still OK
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(report)
# endif # endif
Pout<< "Severe non-orthogonality for face " << faceI Pout<< "Severe non-orthogonality for face " << faceI
<< " between cells " << own[faceI] << " between cells " << own[faceI]
...@@ -1700,7 +1779,7 @@ bool checkFaceDotProduct ...@@ -1700,7 +1779,7 @@ bool checkFaceDotProduct
if( setPtr ) if( setPtr )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(badFace)
# endif # endif
setPtr->insert(faceI); setPtr->insert(faceI);
} }
...@@ -1714,7 +1793,7 @@ bool checkFaceDotProduct ...@@ -1714,7 +1793,7 @@ bool checkFaceDotProduct
if( setPtr ) if( setPtr )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(badFace)
# endif # endif
setPtr->insert(faceI); setPtr->insert(faceI);
} }
...@@ -1726,7 +1805,7 @@ bool checkFaceDotProduct ...@@ -1726,7 +1805,7 @@ bool checkFaceDotProduct
} }
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(minDDotS)
# endif # endif
minDDotS = Foam::min(minDDotS, localMinDDotS); minDDotS = Foam::min(minDDotS, localMinDDotS);
} }
...@@ -1765,7 +1844,7 @@ bool checkFaceDotProduct ...@@ -1765,7 +1844,7 @@ bool checkFaceDotProduct
{ {
// Severe non-orthogonality but mesh still OK // Severe non-orthogonality but mesh still OK
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(report)
# endif # endif
{ {
const scalar angle const scalar angle
...@@ -1784,7 +1863,7 @@ bool checkFaceDotProduct ...@@ -1784,7 +1863,7 @@ bool checkFaceDotProduct
if( setPtr ) if( setPtr )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(badFace)
# endif # endif
setPtr->insert(start+faceI); setPtr->insert(start+faceI);
} }
...@@ -1798,7 +1877,7 @@ bool checkFaceDotProduct ...@@ -1798,7 +1877,7 @@ bool checkFaceDotProduct
if( setPtr ) if( setPtr )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(badFace)
# endif # endif
setPtr->insert(start+faceI); setPtr->insert(start+faceI);
} }
...@@ -1812,7 +1891,7 @@ bool checkFaceDotProduct ...@@ -1812,7 +1891,7 @@ bool checkFaceDotProduct
} }
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(minDDotS)
# endif # endif
minDDotS = Foam::min(minDDotS, localMinDDotS); minDDotS = Foam::min(minDDotS, localMinDDotS);
} }
...@@ -1910,7 +1989,7 @@ bool checkFacePyramids ...@@ -1910,7 +1989,7 @@ bool checkFacePyramids
if( report ) if( report )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(report)
# endif # endif
Pout<< "bool checkFacePyramids(" Pout<< "bool checkFacePyramids("
<< "const bool, const scalar, labelHashSet*) : " << "const bool, const scalar, labelHashSet*) : "
...@@ -1942,7 +2021,7 @@ bool checkFacePyramids ...@@ -1942,7 +2021,7 @@ bool checkFacePyramids
if( report ) if( report )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(report)
# endif # endif
Pout<< "bool checkFacePyramids(" Pout<< "bool checkFacePyramids("
<< "const bool, const scalar, labelHashSet*) : " << "const bool, const scalar, labelHashSet*) : "
...@@ -1965,7 +2044,7 @@ bool checkFacePyramids ...@@ -1965,7 +2044,7 @@ bool checkFacePyramids
if( setPtr ) if( setPtr )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(badFace)
# endif # endif
setPtr->insert(faceI); setPtr->insert(faceI);
} }
...@@ -2229,7 +2308,7 @@ bool checkFaceSkewness ...@@ -2229,7 +2308,7 @@ bool checkFaceSkewness
if( report ) if( report )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(report)
# endif # endif
Pout<< " Severe skewness for face " << faceI Pout<< " Severe skewness for face " << faceI
<< " skewness = " << skewness << endl; << " skewness = " << skewness << endl;
...@@ -2238,7 +2317,7 @@ bool checkFaceSkewness ...@@ -2238,7 +2317,7 @@ bool checkFaceSkewness
if( setPtr ) if( setPtr )
{ {
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(badFace)
# endif # endif
setPtr->insert(faceI); setPtr->insert(faceI);
} }
...@@ -2251,7 +2330,7 @@ bool checkFaceSkewness ...@@ -2251,7 +2330,7 @@ bool checkFaceSkewness
} }
# ifdef USE_OMP # ifdef USE_OMP
# pragma omp critical # pragma omp critical(maxSkew)
# endif # endif
maxSkew = Foam::max(maxSkew, localMaxSkew); maxSkew = Foam::max(maxSkew, localMaxSkew);
} }
...@@ -2420,7 +2499,7 @@ bool checkFaceUniformity ...@@ -2420,7 +2499,7 @@ bool checkFaceUniformity
if( setPtr ) if( setPtr )
{