Newer
Older
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
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 "primitiveMesh.H"
#include "pyramidPointFaceRef.H"
#include "ListOps.H"
#include "unitConversion.H"
#include "SortableList.H"
#include "primitiveMeshTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
Foam::scalar Foam::primitiveMesh::closedThreshold_ = 1.0e-6;
Foam::scalar Foam::primitiveMesh::aspectThreshold_ = 1000;
Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70; // deg
Foam::scalar Foam::primitiveMesh::skewThreshold_ = 4;
Foam::scalar Foam::primitiveMesh::planarCosAngle_ = 1.0e-6;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::primitiveMesh::checkClosedBoundary
(
const vectorField& areas,
const bool report,
const bitSet& internalOrCoupledFaces
DebugInFunction << "Checking if boundary is closed" << endl;
// Loop through all boundary faces and sum up the face area vectors.
// For a closed boundary, this should be zero in all vector components
Vector<solveScalar> sumClosed(Zero);
solveScalar sumMagClosedBoundary = 0;
for (label facei = nInternalFaces(); facei < areas.size(); facei++)
if (!internalOrCoupledFaces.size() || !internalOrCoupledFaces[facei])
sumClosed += Vector<solveScalar>(areas[facei]);
sumMagClosedBoundary += mag(areas[facei]);
reduce(sumClosed, sumOp<Vector<solveScalar>>());
reduce(sumMagClosedBoundary, sumOp<solveScalar>());
Vector<solveScalar> openness = sumClosed/(sumMagClosedBoundary + VSMALL);
if (cmptMax(cmptMag(openness)) > closedThreshold_)
{
if (debug || report)
{
Info<< " ***Boundary openness " << openness
<< " possible hole in boundary description."
<< endl;
}
return true;
}
if (debug || report)
{
Info<< " Boundary openness " << openness << " OK."
<< endl;
bool Foam::primitiveMesh::checkClosedCells
const vectorField& faceAreas,
const scalarField& cellVolumes,
const bool report,
labelHashSet* setPtr,
labelHashSet* aspectSetPtr,
const Vector<label>& meshD
DebugInFunction << "Checking if cells are closed" << endl;
// Check that all cells labels are valid
const cellList& c = cells();
label nErrorClosed = 0;
forAll(c, cI)
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
{
const cell& curCell = c[cI];
if (min(curCell) < 0 || max(curCell) > nFaces())
{
if (setPtr)
{
setPtr->insert(cI);
}
nErrorClosed++;
}
}
if (nErrorClosed > 0)
{
if (debug || report)
{
Info<< " ***Cells with invalid face labels found, number of cells "
<< nErrorClosed << endl;
}
return true;
}
scalarField openness;
scalarField aspectRatio;
primitiveMeshTools::cellClosedness
(
*this,
meshD,
faceAreas,
cellVolumes,
openness,
aspectRatio
);
scalar maxOpennessCell = max(openness);
scalar maxAspectRatio = max(aspectRatio);
forAll(openness, celli)
if (openness[celli] > closedThreshold_)
setPtr->insert(celli);
}
nOpen++;
}
if (aspectRatio[celli] > aspectThreshold_)
{
if (aspectSetPtr)
{
aspectSetPtr->insert(celli);
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
}
nAspect++;
}
}
reduce(nOpen, sumOp<label>());
reduce(maxOpennessCell, maxOp<scalar>());
reduce(nAspect, sumOp<label>());
reduce(maxAspectRatio, maxOp<scalar>());
if (nOpen > 0)
{
if (debug || report)
{
Info<< " ***Open cells found, max cell openness: "
<< maxOpennessCell << ", number of open cells " << nOpen
<< endl;
}
return true;
}
if (nAspect > 0)
{
if (debug || report)
{
Info<< " ***High aspect ratio cells found, Max aspect ratio: "
<< maxAspectRatio
<< ", number of cells " << nAspect
<< endl;
}
return true;
}
if (debug || report)
{
Info<< " Max cell openness = " << maxOpennessCell << " OK." << nl
<< " Max aspect ratio = " << maxAspectRatio << " OK."
<< endl;
}
return false;
}
bool Foam::primitiveMesh::checkFaceAreas
const vectorField& faceAreas,
const bool detailedReport,
labelHashSet* setPtr
) const
{
DebugInFunction << "Checking face area magnitudes" << endl;
const scalarField magFaceAreas(mag(faceAreas));
scalar minArea = GREAT;
scalar maxArea = -GREAT;
forAll(magFaceAreas, facei)
if (magFaceAreas[facei] < VSMALL)
setPtr->insert(facei);
if (detailedReport)
{
if (isInternalFace(facei))
{
Pout<< "Zero or negative face area detected for "
<< "internal face "<< facei << " between cells "
<< faceOwner()[facei] << " and "
<< faceNeighbour()[facei]
<< ". Face area magnitude = " << magFaceAreas[facei]
<< endl;
}
else
{
Pout<< "Zero or negative face area detected for "
<< "boundary face " << facei << " next to cell "
<< faceOwner()[facei] << ". Face area magnitude = "
<< magFaceAreas[facei] << endl;
minArea = min(minArea, magFaceAreas[facei]);
maxArea = max(maxArea, magFaceAreas[facei]);
}
reduce(minArea, minOp<scalar>());
reduce(maxArea, maxOp<scalar>());
if (minArea < VSMALL)
{
if (debug || report)
{
Info<< " ***Zero or negative face area detected. "
"Minimum area: " << minArea << endl;
}
return true;
}
if (debug || report)
{
Info<< " Minimum face area = " << minArea
<< ". Maximum face area = " << maxArea
<< ". Face area magnitudes OK." << endl;
bool Foam::primitiveMesh::checkCellVolumes
const scalarField& vols,
const bool detailedReport,
labelHashSet* setPtr
) const
{
DebugInFunction << "Checking cell volumes" << endl;
scalar minVolume = GREAT;
scalar maxVolume = -GREAT;
label nNegVolCells = 0;
if (vols[celli] < VSMALL)
setPtr->insert(celli);
if (detailedReport)
{
Pout<< "Zero or negative cell volume detected for cell "
<< celli << ". Volume = " << vols[celli] << endl;
nNegVolCells++;
}
minVolume = min(minVolume, vols[celli]);
maxVolume = max(maxVolume, vols[celli]);
}
reduce(minVolume, minOp<scalar>());
reduce(maxVolume, maxOp<scalar>());
reduce(nNegVolCells, sumOp<label>());
if (minVolume < VSMALL)
{
if (debug || report)
{
Info<< " ***Zero or negative cell volume detected. "
<< "Minimum negative volume: " << minVolume
<< ", Number of negative volume cells: " << nNegVolCells
<< endl;
}
return true;
}
if (debug || report)
{
Info<< " Min volume = " << minVolume
<< ". Max volume = " << maxVolume
<< ". Total volume = " << gSum(vols)
<< ". Cell volumes OK." << endl;
bool Foam::primitiveMesh::checkFaceOrthogonality
const vectorField& fAreas,
const vectorField& cellCtrs,
const bool report,
labelHashSet* setPtr
) const
{
DebugInFunction << "Checking mesh non-orthogonality" << endl;
tmp<scalarField> tortho = primitiveMeshTools::faceOrthogonality
(
*this,
fAreas,
cellCtrs
);
const scalarField& ortho = tortho();
// Severe nonorthogonality threshold
const scalar severeNonorthogonalityThreshold =
::cos(degToRad(nonOrthThreshold_));
scalar minDDotS = min(ortho);
scalar sumDDotS = sum(ortho);
label severeNonOrth = 0;
label errorNonOrth = 0;
if (ortho[facei] < severeNonorthogonalityThreshold)
if (ortho[facei] > SMALL)
setPtr->insert(facei);
}
severeNonOrth++;
}
else
{
if (setPtr)
{
setPtr->insert(facei);
}
errorNonOrth++;
}
}
}
reduce(minDDotS, minOp<scalar>());
reduce(sumDDotS, sumOp<scalar>());
reduce(severeNonOrth, sumOp<label>());
reduce(errorNonOrth, sumOp<label>());
if (debug || report)
{
label neiSize = ortho.size();
reduce(neiSize, sumOp<label>());
if (neiSize > 0)
{
if (debug || report)
{
Info<< " Mesh non-orthogonality Max: "
<< radToDeg(::acos(minDDotS))
<< " average: " << radToDeg(::acos(sumDDotS/neiSize))
<< 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;
}
if (debug || report)
{
Info<< " Non-orthogonality check OK." << endl;
bool Foam::primitiveMesh::checkFacePyramids
const pointField& points,
const vectorField& ctrs,
const bool detailedReport,
const scalar minPyrVol,
labelHashSet* setPtr
) const
{
DebugInFunction << "Checking face orientation" << endl;
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
const faceList& f = faces();
scalarField ownPyrVol;
scalarField neiPyrVol;
primitiveMeshTools::facePyramidVolume
(
*this,
points,
ctrs,
ownPyrVol,
neiPyrVol
);
label nErrorPyrs = 0;
forAll(ownPyrVol, facei)
if (ownPyrVol[facei] < minPyrVol)
setPtr->insert(facei);
if (detailedReport)
{
Pout<< "Negative pyramid volume: " << ownPyrVol[facei]
<< " for face " << facei << " " << f[facei]
<< " and owner cell: " << own[facei] << endl
<< "Owner cell vertex labels: "
<< cells()[own[facei]].labels(faces())
nErrorPyrs++;
}
if (isInternalFace(facei))
if (neiPyrVol[facei] < minPyrVol)
setPtr->insert(facei);
if (detailedReport)
{
Pout<< "Negative pyramid volume: " << neiPyrVol[facei]
<< " for face " << facei << " " << f[facei]
<< " and neighbour cell: " << nei[facei] << nl
<< "Neighbour cell vertex labels: "
<< cells()[nei[facei]].labels(faces())
nErrorPyrs++;
}
}
}
reduce(nErrorPyrs, sumOp<label>());
if (nErrorPyrs > 0)
{
if (debug || report)
{
Info<< " ***Error in face pyramids: "
<< nErrorPyrs << " faces are incorrectly oriented."
<< endl;
}
return true;
}
if (debug || report)
{
Info<< " Face pyramids OK." << endl;
bool Foam::primitiveMesh::checkFaceSkewness
const pointField& points,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs,
const bool report,
labelHashSet* setPtr
) const
{
DebugInFunction << "Checking face skewness" << endl;
// Warn if the skew correction vector is more than skewWarning times
// larger than the face area vector
tmp<scalarField> tskewness = primitiveMeshTools::faceSkewness
(
*this,
points,
fCtrs,
fAreas,
cellCtrs
);
const scalarField& skewness = tskewness();
scalar maxSkew = max(skewness);
label nWarnSkew = 0;
forAll(skewness, 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 (skewness[facei] > skewThreshold_)
setPtr->insert(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;
}
if (debug || report)
{
Info<< " Max skewness = " << maxSkew << " OK." << endl;
bool Foam::primitiveMesh::checkFaceAngles
const pointField& points,
const vectorField& faceAreas,
const bool report,
const scalar maxDeg,
labelHashSet* setPtr
) const
{
DebugInFunction << "Checking face angles" << endl;
if (maxDeg < -SMALL || maxDeg > 180+SMALL)
{
FatalErrorInFunction
<< "maxDeg should be [0..180] but is now " << maxDeg
<< exit(FatalError);
}
const scalar maxSin = Foam::sin(degToRad(maxDeg));
tmp<scalarField> tfaceAngles = primitiveMeshTools::faceConcavity
(
maxSin,
*this,
points,
faceAreas
);
const scalarField& faceAngles = tfaceAngles();
scalar maxEdgeSin = max(faceAngles);
label nConcave = 0;
forAll(faceAngles, facei)
if (faceAngles[facei] > SMALL)
setPtr->insert(facei);
}
}
}
reduce(nConcave, sumOp<label>());
reduce(maxEdgeSin, maxOp<scalar>());
if (nConcave > 0)
{
scalar maxConcaveDegr =
radToDeg(Foam::asin(Foam::min(1.0, maxEdgeSin)));
if (debug || report)
{
Info<< " *There are " << nConcave
<< " faces with concave angles between consecutive"
<< " edges. Max concave angle = " << maxConcaveDegr
<< " degrees." << endl;
}
return true;
}
if (debug || report)
{
Info<< " All angles in faces OK." << endl;
bool Foam::primitiveMesh::checkFaceFlatness
const pointField& points,
const vectorField& faceCentres,
const vectorField& faceAreas,
const bool report,
const scalar warnFlatness,
labelHashSet* setPtr
) const
{
DebugInFunction << "Checking face flatness" << endl;
if (warnFlatness < 0 || warnFlatness > 1)
{
FatalErrorInFunction
<< "warnFlatness should be [0..1] but is " << warnFlatness
<< exit(FatalError);
}
const faceList& fcs = faces();
tmp<scalarField> tfaceFlatness = primitiveMeshTools::faceFlatness
(
*this,
points,
faceCentres,
faceAreas
);
const scalarField& faceFlatness = tfaceFlatness();
scalarField magAreas(mag(faceAreas));
scalar minFlatness = GREAT;
scalar sumFlatness = 0;
label nSummed = 0;
label nWarped = 0;
forAll(faceFlatness, facei)
if (fcs[facei].size() > 3 && magAreas[facei] > VSMALL)
sumFlatness += faceFlatness[facei];
minFlatness = min(minFlatness, faceFlatness[facei]);
if (faceFlatness[facei] < warnFlatness)
{
nWarped++;
if (setPtr)
{
setPtr->insert(facei);
}
}
}
}
reduce(nWarped, sumOp<label>());
reduce(minFlatness, minOp<scalar>());
reduce(nSummed, sumOp<label>());
reduce(sumFlatness, sumOp<scalar>());
if (debug || report)
{
if (nSummed > 0)
{
Info<< " Face flatness (1 = flat, 0 = butterfly) : min = "
<< minFlatness << " average = " << sumFlatness / nSummed
<< endl;
{
if (debug || report)
{
Info<< " *There are " << nWarped
<< " faces with ratio between projected and actual area < "
<< warnFlatness << endl;
Info<< " Minimum ratio (minimum flatness, maximum warpage) = "
<< minFlatness << endl;
}
return true;
}
if (debug || report)
{
Info<< " All face flatness OK." << endl;
bool Foam::primitiveMesh::checkConcaveCells
const vectorField& fAreas,
const pointField& fCentres,
const bool report,
labelHashSet* setPtr
) const
{
DebugInFunction << "Checking for concave cells" << endl;
const cellList& c = cells();
const labelList& fOwner = faceOwner();
label nConcaveCells = 0;
const cell& cFaces = c[celli];
bool concave = false;
forAll(cFaces, i)
{
if (concave)
{
break;
}
label fI = cFaces[i];
const point& fC = fCentres[fI];
vector fN = fAreas[fI];
fN /= max(mag(fN), VSMALL);
// Flip normal if required so that it is always pointing out of
// the cell
if (fOwner[fI] != celli)
{
fN *= -1;
}
// Is the centre of any other face of the cell on the
// wrong side of the plane of this face?
forAll(cFaces, j)
{
if (j != i)
{
label fJ = cFaces[j];
const point& pt = fCentres[fJ];
// If the cell is concave, the point will be on the
// positive normal side of the plane of f, defined by
// its centre and normal, and the angle between (pt -
// fC) and fN will be less than 90 degrees, so the dot
// product will be positive.
vector pC = (pt - fC);
pC /= max(mag(pC), VSMALL);
if ((pC & fN) > -planarCosAngle_)
// Concave or planar face
concave = true;
if (setPtr)
setPtr->insert(celli);
nConcaveCells++;
break;
reduce(nConcaveCells, sumOp<label>());
if (nConcaveCells > 0)
Info<< " ***Concave cells (using face planes) found,"
<< " number of cells: " << nConcaveCells << endl;
if (debug || report)
{
Info<< " Concave cell check OK." << endl;
return false;
bool Foam::primitiveMesh::checkUpperTriangular
(
const bool report,
labelHashSet* setPtr
) const
{
DebugInFunction << "Checking face ordering" << endl;
// Check whether internal faces are ordered in the upper triangular order
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
const cellList& c = cells();
label internal = nInternalFaces();
// Has error occurred?
bool error = false;
// Have multiple faces been detected?
label nMultipleCells = false;
// Loop through faceCells once more and make sure that for internal cell
// the first label is smaller
for (label facei = 0; facei < internal; facei++)
if (own[facei] >= nei[facei])
{
error = true;
if (setPtr)
{
setPtr->insert(facei);
}
}
}
// Loop through all cells. For each cell, find the face that is internal
// and add it to the check list (upper triangular order).
// Once the list is completed, check it against the faceCell list
const labelList& curFaces = c[celli];
// Neighbouring cells
SortableList<label> nbr(curFaces.size());
forAll(curFaces, i)
{
label facei = curFaces[i];
if (facei >= nInternalFaces())
{
// Sort last
nbr[i] = labelMax;
}
else
{
label nbrCelli = nei[facei];
if (nbrCelli == celli)
nbrCelli = own[facei];
if (celli < nbrCelli)
nbr[i] = nbrCelli;
}
else
{
// nbrCell is master. Let it handle this face.
nbr[i] = labelMax;
}
}
}
nbr.sort();
// Now nbr holds the cellCells in incremental order. Check:
// - neighbouring cells appear only once. Since nbr is sorted this
// is simple check on consecutive elements
// - faces indexed in same order as nbr are incrementing as well.
label prevCell = nbr[0];