primitiveMeshCheck.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     3.2
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.
 
    foam-extend 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.
 
    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
 
\*---------------------------------------------------------------------------*/
 
#include "primitiveMesh.H"
#include "pyramidPointFaceRef.H"
#include "ListOps.H"
#include "mathematicalConstants.H"
#include "SortableList.H"
 
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::debug::tolerancesSwitch
Foam::primitiveMesh::closedThreshold_
(
    "primitiveMeshClosedThreshold",
    1e-6
);
 
const Foam::debug::tolerancesSwitch
Foam::primitiveMesh::aspectThreshold_
(
    "primitiveMeshAspectThreshold",
    1000
);
 
Foam::debug::tolerancesSwitch
Foam::primitiveMesh::nonOrthThreshold_ // deg
(
    "primitiveMeshNonOrthThreshold",
    70
);
 
const Foam::debug::tolerancesSwitch
Foam::primitiveMesh::skewThreshold_
(
    "primitiveMeshSkewThreshold",
    4
);
 
Foam::debug::tolerancesSwitch
Foam::primitiveMesh::faceAngleThreshold_
(
    "primitiveMeshFaceAngleThreshold",
    10
);
 
const Foam::debug::tolerancesSwitch
Foam::primitiveMesh::faceFlatnessThreshold_
(
    "primitiveMeshFaceFlatnessThreshold",
    0.8
);
 
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
bool Foam::primitiveMesh::checkClosedBoundary(const bool report) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkClosedBoundary("
            << "const bool) const: "
            << "checking whether the 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 sumClosed(vector::zero);
    scalar sumMagClosedBoundary = 0;
 
    const vectorField& areas = faceAreas();
 
    for (label faceI = nInternalFaces(); faceI < areas.size(); faceI++)
    {
        sumClosed += areas[faceI];
        sumMagClosedBoundary += mag(areas[faceI]);
    }
 
    reduce(sumClosed, sumOp<vector>());
    reduce(sumMagClosedBoundary, sumOp<scalar>());
 
    vector openness = sumClosed/(sumMagClosedBoundary + VSMALL);
 
    if (cmptMax(cmptMag(openness)) > closedThreshold_())
    {
        if (debug || report)
        {
            Info<< " ***Boundary openness " << openness
                << " Threshold = " << closedThreshold_()
                << " possible hole in boundary description."
                << endl;
        }
 
        return true;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Boundary openness " << openness
                << " Threshold = " << closedThreshold_()
                << " OK."
                << endl;
        }
 
        return false;
    }
}
 
 
bool Foam::primitiveMesh::checkClosedCells
(
    const bool report,
    labelHashSet* setPtr,
    labelHashSet* aspectSetPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkClosedCells("
            << "const bool, labelHashSet*, labelHashSet*) const: "
            << "checking whether cells are closed" << endl;
    }
 
    // Check that all cells labels are valid
    const cellList& c = cells();
 
    label nErrorClosed = 0;
 
    forAll (c, cI)
    {
        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;
    }
 
    // Check if a face has got the same cell as owner and neighbour
    // HJ, 25/May/2011
    const labelList& own = faceOwner();
    const labelList& nei = faceNeighbour();
 
    label nErrorOwnNei = 0;
 
    // Check internal faces only
    forAll (nei, faceI)
    {
        if (own[faceI] == nei[faceI])
        {
            if (setPtr)
            {
                setPtr->insert(own[faceI]);
            }
 
            nErrorOwnNei++;
        }
    }
 
    if (nErrorOwnNei > 0)
    {
        if (debug || report)
        {
            Info<< " ***Faces declaring same cell as owner and neighbour "
                << "found, number of faces "
                << nErrorOwnNei << endl;
        }
 
        return true;
    }
 
    // Loop through cell faces and sum up the face area vectors for each cell.
    // This should be zero in all vector components
 
    vectorField sumClosed(nCells(), vector::zero);
    vectorField sumMagClosed(nCells(), vector::zero);
 
    const vectorField& areas = faceAreas();
 
    forAll (own, faceI)
    {
        // Add to owner
        sumClosed[own[faceI]] += areas[faceI];
        sumMagClosed[own[faceI]] += cmptMag(areas[faceI]);
    }
 
    forAll (nei, faceI)
    {
        // Subtract from neighbour
        sumClosed[nei[faceI]] -= areas[faceI];
        sumMagClosed[nei[faceI]] += cmptMag(areas[faceI]);
    }
 
    label nOpen = 0;
    scalar maxOpennessCell = 0;
 
    label nAspect = 0;
    scalar maxAspectRatio = 0;
 
    const scalarField& vols = cellVolumes();
 
    // Check the sums
    forAll (sumClosed, cellI)
    {
        scalar maxOpenness = 0;
 
        for(direction cmpt=0; cmpt<vector::nComponents; cmpt++)
        {
            maxOpenness = max
            (
                maxOpenness,
                mag(sumClosed[cellI][cmpt])
               /(sumMagClosed[cellI][cmpt] + VSMALL)
            );
        }
 
        maxOpennessCell = max(maxOpennessCell, maxOpenness);
 
        if (maxOpenness > closedThreshold_())
        {
            if (setPtr)
            {
                setPtr->insert(cellI);
            }
 
            nOpen++;
        }
 
        // Calculate the aspect ration as the maximum of Cartesian component
        // aspect ratio to the total area hydraulic area aspect ratio
        scalar aspectRatio = max
        (
            cmptMax(sumMagClosed[cellI])
           /(cmptMin(sumMagClosed[cellI]) + VSMALL),
            1.0/6.0*cmptSum(sumMagClosed[cellI])/
            Foam::pow(Foam::max(vols[cellI], SMALL), 2.0/3.0)
        );
 
        maxAspectRatio = max(maxAspectRatio, aspectRatio);
 
        if (aspectRatio > aspectThreshold_())
        {
            if (aspectSetPtr)
            {
                aspectSetPtr->insert(cellI);
            }
 
            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
                << " Threshold = " << closedThreshold_()
                << endl;
        }
 
        return true;
    }
 
    if (nAspect > 0)
    {
        if (debug || report)
        {
            Info<< " ***High aspect ratio cells found, Max aspect ratio: "
                << maxAspectRatio
                << ", number of cells " << nAspect
                << " Threshold = " << aspectThreshold_()
                << 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 bool report,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkFaceAreas("
            << "const bool, labelHashSet*) const: "
            << "checking face area magnitudes" << endl;
    }
 
    const scalarField magFaceAreas = mag(faceAreas());
 
    scalar minArea = GREAT;
    scalar maxArea = -GREAT;
 
    forAll (magFaceAreas, faceI)
    {
        if (magFaceAreas[faceI] < VSMALL)
        {
            if (setPtr)
            {
                setPtr->insert(faceI);
            }
        }
 
        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;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Minumum face area = " << minArea
                << ". Maximum face area = " << maxArea
                << ".  Face area magnitudes OK." << endl;
        }
 
        return false;
    }
}
 
 
bool Foam::primitiveMesh::checkCellVolumes
(
    const bool report,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkCellVolumes("
            << "const bool, labelHashSet*) const: "
            << "checking cell volumes" << endl;
    }
 
    const scalarField& vols = cellVolumes();
 
    scalar minVolume = GREAT;
    scalar maxVolume = -GREAT;
 
    label nNegVolCells = 0;
 
    forAll (vols, cellI)
    {
        if (vols[cellI] < VSMALL)
        {
            if (setPtr)
            {
                setPtr->insert(cellI);
            }
 
            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;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Min volume = " << minVolume
                << ". Max volume = " << maxVolume
                << ".  Total volume = " << gSum(vols)
                << ".  Cell volumes OK." << endl;
        }
 
        return false;
    }
}
 
 
bool Foam::primitiveMesh::checkFaceOrthogonality
(
    const bool report,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkFaceOrthogonality("
            << "const bool, labelHashSet*) const: "
            << "checking mesh non-orthogonality" << endl;
    }
 
    // for all internal faces check theat the d dot S product is positive
    const vectorField& centres = cellCentres();
    const vectorField& areas = faceAreas();
 
    const labelList& own = faceOwner();
    const labelList& nei = faceNeighbour();
 
    // Severe nonorthogonality threshold
    const scalar severeNonorthogonalityThreshold =
        ::cos(nonOrthThreshold_()/180.0*mathematicalConstant::pi);
 
    scalar minDDotS = GREAT;
 
    scalar sumDDotS = 0;
 
    label severeNonOrth = 0;
 
    label errorNonOrth = 0;
 
    forAll (nei, faceI)
    {
        vector d = centres[nei[faceI]] - centres[own[faceI]];
        const vector& s = areas[faceI];
 
        scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
 
        if (dDotS < severeNonorthogonalityThreshold)
        {
            if (dDotS > SMALL)
            {
                if (setPtr)
                {
                    setPtr->insert(faceI);
                }
 
                severeNonOrth++;
            }
            else
            {
                if (setPtr)
                {
                    setPtr->insert(faceI);
                }
 
                errorNonOrth++;
            }
        }
 
        if (dDotS < minDDotS)
        {
            minDDotS = dDotS;
        }
 
        sumDDotS += dDotS;
    }
 
    reduce(minDDotS, minOp<scalar>());
    reduce(sumDDotS, sumOp<scalar>());
    reduce(severeNonOrth, sumOp<label>());
    reduce(errorNonOrth, sumOp<label>());
 
    if (debug || report)
    {
        label neiSize = nei.size();
        reduce(neiSize, sumOp<label>());
 
        if (neiSize > 0)
        {
            if (debug || report)
            {
                Info<< "    Mesh non-orthogonality Max: "
                    << ::acos(minDDotS)/mathematicalConstant::pi*180.0
                    << " average: " <<
                    ::acos(sumDDotS/neiSize)/mathematicalConstant::pi*180.0
                    << " Threshold = " << nonOrthThreshold_()
                    << 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;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Non-orthogonality check OK." << endl;
        }
 
        return false;
    }
}
 
 
bool Foam::primitiveMesh::checkFacePyramids
(
    const bool report,
    const scalar minPyrVol,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkFacePyramids("
            << "const bool, const scalar, labelHashSet*) const: "
            << "checking face orientation" << endl;
    }
 
    // check whether face area vector points to the cell with higher label
    const vectorField& ctrs = cellCentres();
 
    const labelList& own = faceOwner();
    const labelList& nei = faceNeighbour();
 
    const faceList& f = faces();
 
    const pointField& p = points();
 
    label nErrorPyrs = 0;
 
    forAll (f, faceI)
    {
        // Create the owner pyramid - it will have negative volume
        scalar pyrVol = pyramidPointFaceRef(f[faceI], ctrs[own[faceI]]).mag(p);
 
        if (pyrVol > -minPyrVol)
        {
            if (setPtr)
            {
                setPtr->insert(faceI);
            }
 
            nErrorPyrs++;
        }
 
        if (isInternalFace(faceI))
        {
            // Create the neighbour pyramid - it will have positive volume
            scalar pyrVol =
                pyramidPointFaceRef(f[faceI], ctrs[nei[faceI]]).mag(p);
 
            if (pyrVol < minPyrVol)
            {
                if (setPtr)
                {
                    setPtr->insert(faceI);
                }
 
                nErrorPyrs++;
            }
        }
    }
 
    reduce(nErrorPyrs, sumOp<label>());
 
    if (nErrorPyrs > 0)
    {
        if (debug || report)
        {
            Info<< " ***Error in face pyramids: "
                << nErrorPyrs << " faces are incorrectly oriented."
                << endl;
        }
 
        return true;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Face pyramids OK." << endl;
        }
 
        return false;
    }
}
 
 
bool Foam::primitiveMesh::checkFaceSkewness
(
    const bool report,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkFaceSkewnesss("
            << "const bool, labelHashSet*) const: "
            << "checking face skewness" << endl;
    }
 
    // Warn if the skew correction vector is more than skewWarning times
    // larger than the face area vector
 
    const pointField& p = points();
    const faceList& fcs = faces();
 
    const labelList& own = faceOwner();
    const labelList& nei = faceNeighbour();
    const vectorField& cellCtrs = cellCentres();
    const vectorField& faceCtrs = faceCentres();
    const vectorField& fAreas = faceAreas();
 
    scalar maxSkew = 0;
    label nWarnSkew = 0;
 
    forAll(nei, faceI)
    {
        vector Cpf = faceCtrs[faceI] - cellCtrs[own[faceI]];
        vector d = cellCtrs[nei[faceI]] - cellCtrs[own[faceI]];
 
        // Skewness vector
        vector sv =
            Cpf - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + SMALL))*d;
        vector svHat = sv/(mag(sv) + VSMALL);
 
        // Normalisation distance calculated as the approximate distance
        // from the face centre to the edge of the face in the direction of
        // the skewness
        scalar fd = 0.2*mag(d) + VSMALL;
        const face& f = fcs[faceI];
        forAll(f, pi)
        {
            fd = max(fd, mag(svHat & (p[f[pi]] - faceCtrs[faceI])));
        }
 
        // Normalised skewness
        scalar skewness = mag(sv)/fd;
 
        // 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 > skewThreshold_())
        {
            if (setPtr)
            {
                setPtr->insert(faceI);
            }
 
            nWarnSkew++;
        }
 
        if(skewness > maxSkew)
        {
            maxSkew = skewness;
        }
    }
 
 
    // Boundary faces: consider them to have only skewness error.
    // (i.e. treat as if mirror cell on other side)
 
    for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
    {
        vector Cpf = faceCtrs[faceI] - cellCtrs[own[faceI]];
 
        vector normal = fAreas[faceI];
        normal /= mag(normal) + VSMALL;
        vector d = normal*(normal & Cpf);
 
 
        // Skewness vector
        vector sv = Cpf - ((fAreas[faceI]&Cpf)/((fAreas[faceI]&d)+VSMALL))*d;
        vector svHat = sv/(mag(sv) + VSMALL);
 
        // Normalisation distance calculated as the approximate distance
        // from the face centre to the edge of the face in the direction of
        // the skewness
        scalar fd = 0.4*mag(d) + VSMALL;
        const face& f = fcs[faceI];
        forAll(f, pi)
        {
            fd = max(fd, mag(svHat & (p[f[pi]] - faceCtrs[faceI])));
        }
 
        // Normalised skewness
        scalar skewness = mag(sv)/fd;
 
        // 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 > skewThreshold_())
        {
            if (setPtr)
            {
                setPtr->insert(faceI);
            }
 
            nWarnSkew++;
        }
 
        if(skewness > maxSkew)
        {
            maxSkew = skewness;
        }
    }
 
 
    reduce(maxSkew, maxOp<scalar>());
    reduce(nWarnSkew, sumOp<label>());
 
    if (nWarnSkew > 0)
    {
        if (debug || report)
        {
            Info<< " ***Max skewness = " << maxSkew
                << ", " << nWarnSkew << " highly skew faces detected"
                << " Threshold = " << skewThreshold_()
                << endl;
        }
 
        return true;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Max skewness = " << maxSkew << " OK." << endl;
        }
 
        return false;
    }
}
 
 
bool Foam::primitiveMesh::checkPoints
(
    const bool report,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkPoints"
            << "(const bool, labelHashSet*) const: "
            << "checking points" << endl;
    }
 
    label nFaceErrors = 0;
    label nCellErrors = 0;
 
    const labelListList& pf = pointFaces();
 
    forAll (pf, pointI)
    {
        if (pf[pointI].empty())
        {
            if (setPtr)
            {
                setPtr->insert(pointI);
            }
 
            nFaceErrors++;
        }
    }
 
 
    forAll (pf, pointI)
    {
        const labelList& pc = pointCells(pointI);
 
        if (pc.empty())
        {
            if (setPtr)
            {
                setPtr->insert(pointI);
            }
 
            nCellErrors++;
        }
    }
 
    reduce(nFaceErrors, sumOp<label>());
    reduce(nCellErrors, sumOp<label>());
 
    if (nFaceErrors > 0 || nCellErrors > 0)
    {
        if (debug || report)
        {
            Info<< " ***Unused points found in the mesh, "
                   "number unused by faces: " << nFaceErrors
                << " number unused by cells: " << nCellErrors
                << endl;
        }
 
        return true;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Point usage OK." << endl;
        }
 
        return false;
    }
}
 
 
// Check convexity of angles in a face. Allow a slight non-convexity.
// E.g. maxDeg = 10 allows for angles < 190 (or 10 degrees concavity)
// (if truly concave and points not visible from face centre the face-pyramid
//  check in checkMesh will fail)
bool Foam::primitiveMesh::checkFaceAngles
(
    const bool report,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkFaceAngles"
            << "(const bool, labelHashSet*) const: "
            << "checking face angles" << endl;
    }
 
    if (faceAngleThreshold_() < -SMALL || faceAngleThreshold_() > 180 + SMALL)
    {
        FatalErrorIn
        (
            "primitiveMesh::checkFaceAngles(const bool, labelHashSet*)"
        )   << "faceAngleThreshold_ should be [0..180] but is now "
            << faceAngleThreshold_()
            << exit(FatalError);
    }
 
    const scalar maxSin =
        Foam::sin(faceAngleThreshold_()/180.0*mathematicalConstant::pi);
 
    const pointField& p = points();
    const faceList& fcs = faces();
    vectorField faceNormals(faceAreas());
    faceNormals /= mag(faceNormals) + VSMALL;
 
    scalar maxEdgeSin = 0.0;
 
    label nConcave = 0;
 
    label errorFaceI = -1;
 
    forAll(fcs, faceI)
    {
        const face& f = fcs[faceI];
 
        // Get edge from f[0] to f[size-1];
        vector ePrev(p[f[0]] - p[f[f.size()-1]]);
        scalar magEPrev = mag(ePrev);
        ePrev /= magEPrev + VSMALL;
 
        forAll(f, fp0)
        {
            // Get vertex after fp
            label fp1 = f.fcIndex(fp0);
 
            // Normalized vector between two consecutive points
            vector e10(p[f[fp1]] - p[f[fp0]]);
            scalar magE10 = mag(e10);
            e10 /= magE10 + VSMALL;
 
            if (magEPrev > SMALL && magE10 > SMALL)
            {
                vector edgeNormal = ePrev ^ e10;
                scalar magEdgeNormal = mag(edgeNormal);
 
                if (magEdgeNormal < maxSin)
                {
                    // Edges (almost) aligned -> face is ok.
                }
                else
                {
                    // Check normal
                    edgeNormal /= magEdgeNormal;
 
                    if ((edgeNormal & faceNormals[faceI]) < SMALL)
                    {
                        if (faceI != errorFaceI)
                        {
                            // Count only one error per face.
                            errorFaceI = faceI;
                            nConcave++;
                        }
 
                        if (setPtr)
                        {
                            setPtr->insert(faceI);
                        }
 
                        maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
                    }
                }
            }
 
            ePrev = e10;
            magEPrev = magE10;
        }
    }
 
    reduce(nConcave, sumOp<label>());
    reduce(maxEdgeSin, maxOp<scalar>());
 
    if (nConcave > 0)
    {
        scalar maxConcaveDegr =
            Foam::asin(Foam::min(1.0, maxEdgeSin))
           *180.0/mathematicalConstant::pi;
 
        if (debug || report)
        {
            Info<< "   *There are " << nConcave
                << " faces with concave angles between consecutive"
                << " edges. Max concave angle = " << maxConcaveDegr
                << " degrees." << endl;
        }
 
        return true;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    All angles in faces OK." << endl;
        }
 
        return false;
    }
}
 
 
// Check warpage of faces. Is calculated as the difference between areas of
// individual triangles and the overall area of the face (which ifself is
// is the average of the areas of the individual triangles).
bool Foam::primitiveMesh::checkFaceFlatness
(
    const bool report,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkFaceFlatness"
            << "(const bool, labelHashSet*) const: "
            << "checking face flatness" << endl;
    }
 
    if (faceFlatnessThreshold_() < 0 || faceFlatnessThreshold_() > 1)
    {
        FatalErrorIn
        (
            "primitiveMesh::checkFaceFlatness"
            "(const bool, labelHashSet*)"
        )   << "faceFlatnessThreshold_ should be [0..1] but is now "
            << faceFlatnessThreshold_()
            << exit(FatalError);
    }
 
 
    const pointField& p = points();
    const faceList& fcs = faces();
    const pointField& fctrs = faceCentres();
 
    // Areas are calculated as the sum of areas. (see
    // primitiveMeshFaceCentresAndAreas.C)
    scalarField magAreas(mag(faceAreas()));
 
    label nWarped = 0;
 
    scalar minFlatness = GREAT;
    scalar sumFlatness = 0;
    label nSummed = 0;
 
    forAll(fcs, faceI)
    {
        const face& f = fcs[faceI];
 
        if (f.size() > 3 && magAreas[faceI] > VSMALL)
        {
            const point& fc = fctrs[faceI];
 
            // Calculate the sum of magnitude of areas and compare to magnitude
            // of sum of areas.
 
            scalar sumA = 0.0;
 
            forAll(f, fp)
            {
                const point& thisPoint = p[f[fp]];
                const point& nextPoint = p[f.nextLabel(fp)];
 
                // Triangle around fc.
                vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
                sumA += mag(n);
            }
 
            scalar flatness = magAreas[faceI] / (sumA+VSMALL);
 
            sumFlatness += flatness;
            nSummed++;
 
            minFlatness = min(minFlatness, flatness);
 
            if (flatness < faceFlatnessThreshold_())
            {
                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) : average = "
                << sumFlatness / nSummed << "  min = " << minFlatness << endl;
        }
    }
 
 
    if (nWarped > 0)
    {
        if (debug || report)
        {
            Info<< "   *There are " << nWarped
                << " faces with ratio between projected and actual area < "
                << faceFlatnessThreshold_() << endl;
 
            Info<< "    Minimum ratio (minimum flatness, maximum warpage) = "
                << minFlatness << endl;
        }
 
        return true;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    All face flatness OK." << endl;
        }
 
        return false;
    }
}
 
 
// Check 1D/2Dness of edges. Gets passed the non-empty directions and
// checks all edges in the mesh whether they:
// - have no component in a non-empty direction or
// - are only in a singe non-empty direction.
// Empty direction info is passed in as a vector of labels (synchronised)
// which are 1 if the direction is non-empty, 0 if it is.
bool Foam::primitiveMesh::checkEdgeAlignment
(
    const bool report,
    const Vector<label>& directions,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkEdgeAlignment("
            << "const bool, const Vector<label>&, labelHashSet*) const: "
            << "checking edge alignment" << endl;
    }
 
    label nDirs = 0;
    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
    {
        if (directions[cmpt] == 1)
        {
            nDirs++;
        }
        else if (directions[cmpt] != 0)
        {
            FatalErrorIn
            (
                "primitiveMesh::checkEdgeAlignment"
                "(const bool, const Vector<label>&, labelHashSet*)"
            )   << "directions should contain 0 or 1 but is now " << directions
                << exit(FatalError);
        }
    }
 
    if (nDirs == vector::nComponents)
    {
        return false;
    }
 
 
 
    const pointField& p = points();
    const faceList& fcs = faces();
 
    EdgeMap<label> edgesInError;
 
    forAll(fcs, faceI)
    {
        const face& f = fcs[faceI];
 
        forAll(f, fp)
        {
            label p0 = f[fp];
            label p1 = f.nextLabel(fp);
            if (p0 < p1)
            {
                vector d(p[p1]-p[p0]);
                scalar magD = mag(d);
 
                if (magD > ROOTVSMALL)
                {
                    d /= magD;
 
                    // Check how many empty directions are used by the edge.
                    label nEmptyDirs = 0;
                    label nNonEmptyDirs = 0;
                    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
                    {
                        if (mag(d[cmpt]) > 1e-6)
                        {
                            if (directions[cmpt] == 0)
                            {
                                nEmptyDirs++;
                            }
                            else
                            {
                                nNonEmptyDirs++;
                            }
                        }
                    }
 
                    if (nEmptyDirs == 0)
                    {
                        // Purely in ok directions.
                    }
                    else if (nEmptyDirs == 1)
                    {
                        // Ok if purely in empty directions.
                        if (nNonEmptyDirs > 0)
                        {
                            edgesInError.insert(edge(p0, p1), faceI);
                        }
                    }
                    else if (nEmptyDirs > 1)
                    {
                        // Always an error
                        edgesInError.insert(edge(p0, p1), faceI);
                    }
                }
            }
        }
    }
 
    label nErrorEdges = returnReduce(edgesInError.size(), sumOp<label>());
 
    if (nErrorEdges > 0)
    {
        if (debug || report)
        {
            Info<< " ***Number of edges not aligned with or perpendicular to "
                << "non-empty directions: " << nErrorEdges << endl;
        }
 
        if (setPtr)
        {
            setPtr->resize(2*edgesInError.size());
            forAllConstIter(EdgeMap<label>, edgesInError, iter)
            {
                setPtr->insert(iter.key()[0]);
                setPtr->insert(iter.key()[1]);
            }
        }
 
        return true;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    All edges aligned with or perpendicular to "
                << "non-empty directions." << endl;
        }
        return false;
    }
}
 
 
bool Foam::primitiveMesh::checkUpperTriangular
(
    const bool report,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkUpperTriangular("
            << "const bool, labelHashSet*) const: "
            << "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
 
    forAll (c, cellI)
    {
        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)
                {
                    // cellI is master
                    nbr[i] = nbrCellI;
                }
                else
                {
                    // nbrCell is master. Let it handle this face.
                    nbr[i] = labelMax;
                }
            }
        }
 
        if (!nbr.empty())
        {
            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];
            label prevFace = curFaces[nbr.indices()[0]];
 
            bool hasMultipleFaces = false;
 
            for (label i = 1; i < nbr.size(); i++)
            {
                label thisCell = nbr[i];
                label thisFace = curFaces[nbr.indices()[i]];
 
                if (thisCell == labelMax)
                {
                    break;
                }
 
                if (thisCell == prevCell)
                {
                    hasMultipleFaces = true;
 
                    if (setPtr)
                    {
                        setPtr->insert(prevFace);
                        setPtr->insert(thisFace);
                    }
                }
                else if (thisFace < prevFace)
                {
                    error = true;
 
                    if (setPtr)
                    {
                        setPtr->insert(thisFace);
                    }
                }
 
                prevCell = thisCell;
                prevFace = thisFace;
            }
 
            if (hasMultipleFaces)
            {
                nMultipleCells++;
            }
        }
    }
 
    reduce(error, orOp<bool>());
    reduce(nMultipleCells, sumOp<label>());
 
    if ((debug || report) && nMultipleCells > 0)
    {
        Info<< "  <<Found " << nMultipleCells
            << " neighbouring cells with multiple inbetween faces." << endl;
    }
 
    if (error)
    {
        if (debug || report)
        {
            Info<< " ***Faces not in upper triangular order." << endl;
        }
 
        return true;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Upper triangular ordering OK." << endl;
        }
 
        return false;
    }
}
 
 
bool Foam::primitiveMesh::checkCellsZipUp
(
    const bool report,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkCellsZipUp("
            << "const bool, labelHashSet*) const: "
            << "checking topological cell openness" << endl;
    }
 
    label nOpenCells = 0;
 
    const faceList& f = faces();
    const cellList& c = cells();
 
    forAll (c, cellI)
    {
        const labelList& curFaces = c[cellI];
 
        const edgeList cellEdges = c[cellI].edges(f);
 
        labelList edgeUsage(cellEdges.size(), 0);
 
        forAll (curFaces, faceI)
        {
            edgeList curFaceEdges = f[curFaces[faceI]].edges();
 
            forAll (curFaceEdges, faceEdgeI)
            {
                const edge& curEdge = curFaceEdges[faceEdgeI];
 
                forAll (cellEdges, cellEdgeI)
                {
                    if (cellEdges[cellEdgeI] == curEdge)
                    {
                        edgeUsage[cellEdgeI]++;
                        break;
                    }
                }
            }
        }
 
        edgeList singleEdges(cellEdges.size());
        label nSingleEdges = 0;
 
        forAll (edgeUsage, edgeI)
        {
            if (edgeUsage[edgeI] == 1)
            {
                singleEdges[nSingleEdges] = cellEdges[edgeI];
                nSingleEdges++;
            }
            else if (edgeUsage[edgeI] != 2)
            {
                if (setPtr)
                {
                    setPtr->insert(cellI);
                }
            }
        }
 
        if (nSingleEdges > 0)
        {
            if (setPtr)
            {
                setPtr->insert(cellI);
            }
 
            nOpenCells++;
        }
    }
 
    reduce(nOpenCells, sumOp<label>());
 
    if (nOpenCells > 0)
    {
        if (debug || report)
        {
            Info<< " ***Open cells found, number of cells: " << nOpenCells
                << ". This problem may be fixable using the zipUpMesh utility."
                << endl;
        }
 
        return true;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Topological cell zip-up check OK." << endl;
        }
 
        return false;
    }
}
 
 
// Vertices of face within point range and unique.
bool Foam::primitiveMesh::checkFaceVertices
(
    const bool report,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkFaceVertices("
            << "const bool, labelHashSet*) const: "
            << "checking face vertices" << endl;
    }
 
    // Check that all vertex labels are valid
    const faceList& f = faces();
 
    label nErrorFaces = 0;
 
    forAll (f, fI)
    {
        const face& curFace = f[fI];
 
        if (min(curFace) < 0 || max(curFace) > nPoints())
        {
            if (setPtr)
            {
                setPtr->insert(fI);
            }
 
            nErrorFaces++;
        }
 
        // Uniqueness of vertices
        labelHashSet facePoints(2*curFace.size());
 
        forAll(curFace, fp)
        {
            bool inserted = facePoints.insert(curFace[fp]);
 
            if (!inserted)
            {
                if (setPtr)
                {
                    setPtr->insert(fI);
                }
 
                nErrorFaces++;
            }
        }
    }
 
    reduce(nErrorFaces, sumOp<label>());
 
    if (nErrorFaces > 0)
    {
        if (debug || report)
        {
            Info<< "    Faces with invalid vertex labels found, "
                << " number of faces: " << nErrorFaces << endl;
        }
 
        return true;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Face vertices OK." << endl;
        }
 
        return false;
    }
}
 
 
// Check if all points on face are shared between faces.
bool Foam::primitiveMesh::checkDuplicateFaces
(
    const label faceI,
    const Map<label>& nCommonPoints,
    label& nBaffleFaces,
    labelHashSet* setPtr
) const
{
    bool error = false;
 
    forAllConstIter(Map<label>, nCommonPoints, iter)
    {
        label nbFaceI = iter.key();
        label nCommon = iter();
 
        const face& curFace = faces()[faceI];
        const face& nbFace = faces()[nbFaceI];
 
        if (nCommon == nbFace.size() || nCommon == curFace.size())
        {
            if (nbFace.size() != curFace.size())
            {
                error = true;
            }
            else
            {
                nBaffleFaces++;
            }
 
            if (setPtr)
            {
                setPtr->insert(faceI);
                setPtr->insert(nbFaceI);
            }
        }
    }
 
    return error;
}
 
 
// Check that shared points are in consecutive order.
bool Foam::primitiveMesh::checkCommonOrder
(
    const label faceI,
    const Map<label>& nCommonPoints,
    labelHashSet* setPtr
) const
{
    bool error = false;
 
    forAllConstIter(Map<label>, nCommonPoints, iter)
    {
        label nbFaceI = iter.key();
        label nCommon = iter();
 
        const face& curFace = faces()[faceI];
        const face& nbFace = faces()[nbFaceI];
 
        if
        (
            nCommon >= 2
         && nCommon != nbFace.size()
         && nCommon != curFace.size()
        )
        {
            forAll(curFace, fp)
            {
                // Get the index in the neighbouring face shared with curFace
                label nb = findIndex(nbFace, curFace[fp]);
 
                if (nb != -1)
                {
 
                    // Check the whole face from nb onwards for shared vertices
                    // with neighbouring face. Rule is that any shared vertices
                    // should be consecutive on both faces i.e. if they are
                    // vertices fp,fp+1,fp+2 on one face they should be
                    // vertices nb, nb+1, nb+2 (or nb+2, nb+1, nb) on the
                    // other face.
 
 
                    // Vertices before and after on curFace
                    label fpPlus1 = curFace.fcIndex(fp);
                    label fpMin1  = curFace.rcIndex(fp);
 
                    // Vertices before and after on nbFace
                    label nbPlus1 = nbFace.fcIndex(nb);
                    label nbMin1  = nbFace.rcIndex(nb);
 
                    // Find order of walking by comparing next points on both
                    // faces.
                    label curInc = labelMax;
                    label nbInc = labelMax;
 
                    if (nbFace[nbPlus1] == curFace[fpPlus1])
                    {
                        curInc = 1;
                        nbInc = 1;
                    }
                    else if (nbFace[nbPlus1] == curFace[fpMin1])
                    {
                        curInc = -1;
                        nbInc = 1;
                    }
                    else if (nbFace[nbMin1] == curFace[fpMin1])
                    {
                        curInc = -1;
                        nbInc = -1;
                    }
                    else
                    {
                        curInc = 1;
                        nbInc = -1;
                    }
 
 
                    // Pass1: loop until start of common vertices found.
                    label curNb = nb;
                    label curFp = fp;
 
                    do
                    {
                        curFp += curInc;
 
                        if (curFp >= curFace.size())
                        {
                            curFp = 0;
                        }
                        else if (curFp < 0)
                        {
                            curFp = curFace.size()-1;
                        }
 
                        curNb += nbInc;
 
                        if (curNb >= nbFace.size())
                        {
                            curNb = 0;
                        }
                        else if (curNb < 0)
                        {
                            curNb = nbFace.size()-1;
                        }
                    } while (curFace[curFp] == nbFace[curNb]);
 
 
                    // Pass2: check equality walking from curFp, curNb
                    // in opposite order.
 
                    curInc = -curInc;
                    nbInc = -nbInc;
 
                    for (label commonI = 0; commonI < nCommon; commonI++)
                    {
                        curFp += curInc;
 
                        if (curFp >= curFace.size())
                        {
                            curFp = 0;
                        }
                        else if (curFp < 0)
                        {
                            curFp = curFace.size()-1;
                        }
 
                        curNb += nbInc;
 
                        if (curNb >= nbFace.size())
                        {
                            curNb = 0;
                        }
                        else if (curNb < 0)
                        {
                            curNb = nbFace.size()-1;
                        }
 
                        if (curFace[curFp] != nbFace[curNb])
                        {
                            if (setPtr)
                            {
                                setPtr->insert(faceI);
                                setPtr->insert(nbFaceI);
                            }
 
                            error = true;
 
                            break;
                        }
                    }
 
 
                    // Done the curFace - nbFace combination.
                    break;
                }
            }
        }
    }
 
    return error;
}
 
 
// Checks common vertices between faces. If more than 2 they should be
// consecutive on both faces.
bool Foam::primitiveMesh::checkFaceFaces
(
    const bool report,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkFaceFaces(const bool, labelHashSet*)"
            << " const: " << "checking face-face connectivity" << endl;
    }
 
    const labelListList& pf = pointFaces();
 
    label nBaffleFaces = 0;
    label nErrorDuplicate = 0;
    label nErrorOrder = 0;
    Map<label> nCommonPoints(100);
 
    for (label faceI = 0; faceI < nFaces(); faceI++)
    {
        const face& curFace = faces()[faceI];
 
        // Calculate number of common points between current faceI and
        // neighbouring face. Store on map.
        nCommonPoints.clear();
 
        forAll(curFace, fp)
        {
            label pointI = curFace[fp];
 
            const labelList& nbs = pf[pointI];
 
            forAll(nbs, nbI)
            {
                label nbFaceI = nbs[nbI];
 
                if (faceI < nbFaceI)
                {
                    // Only check once for each combination of two faces.
 
                    Map<label>::iterator fnd = nCommonPoints.find(nbFaceI);
 
                    if (fnd == nCommonPoints.end())
                    {
                        // First common vertex found.
                        nCommonPoints.insert(nbFaceI, 1);
                    }
                    else
                    {
                        fnd()++;
                    }
                }
            }
        }
 
        // Perform various checks on common points
 
        // Check all vertices shared (duplicate point)
        if (checkDuplicateFaces(faceI, nCommonPoints, nBaffleFaces, setPtr))
        {
            nErrorDuplicate++;
        }
 
        // Check common vertices are consecutive on both faces
        if (checkCommonOrder(faceI, nCommonPoints, setPtr))
        {
            nErrorOrder++;
        }
    }
 
    reduce(nBaffleFaces, sumOp<label>());
    reduce(nErrorDuplicate, sumOp<label>());
    reduce(nErrorOrder, sumOp<label>());
 
    if (nBaffleFaces)
    {
        Info<< "    Number of identical duplicate faces (baffle faces): "
            << nBaffleFaces << endl;
    }
 
    if (nErrorDuplicate > 0 || nErrorOrder > 0)
    {
        if (nErrorDuplicate > 0)
        {
            Info<< " ***Number of duplicate (not baffle) faces found: "
                << nErrorDuplicate << endl;
        }
 
        if (nErrorOrder > 0)
        {
            Info<< " ***Number of faces with non-consecutive shared points: "
                << nErrorOrder << endl;
        }
 
        return true;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Face-face connectivity OK." << endl;
        }
 
        return false;
    }
}
 
 
// Checks cells with 1 or less internal faces. Give numerical problems.
bool Foam::primitiveMesh::checkCellDeterminant
(
    const bool report,    // report,
    labelHashSet* setPtr  // setPtr
) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkCellDeterminant(const bool"
            << ", labelHashSet*) const: "
            << "checking for under-determined cells" << endl;
    }
 
    const cellList& c = cells();
 
    label nErrorCells = 0;
 
    scalar minDet = GREAT;
    scalar sumDet = 0;
    label nSummed = 0;
 
    forAll (c, cellI)
    {
        const labelList& curFaces = c[cellI];
 
        // Calculate local normalization factor
        scalar avgArea = 0;
 
        label nInternalFaces = 0;
 
        forAll(curFaces, i)
        {
            if (isInternalFace(curFaces[i]))
            {
                avgArea += mag(faceAreas()[curFaces[i]]);
 
                nInternalFaces++;
            }
        }
 
        if (nInternalFaces == 0)
        {
            if (setPtr)
            {
                setPtr->insert(cellI);
            }
 
            nErrorCells++;
        }
        else
        {
            avgArea /= nInternalFaces;
 
            symmTensor areaTensor(symmTensor::zero);
 
            forAll(curFaces, i)
            {
                if (isInternalFace(curFaces[i]))
                {
                    areaTensor += sqr(faceAreas()[curFaces[i]]/avgArea);
                }
            }
 
            scalar determinant = mag(det(areaTensor));
 
            minDet = min(determinant, minDet);
            sumDet += determinant;
            nSummed++;
 
            if (determinant < 1e-3)
            {
                if (setPtr)
                {
                    setPtr->insert(cellI);
                }
 
                nErrorCells++;
            }
        }
    }
 
    reduce(nErrorCells, sumOp<label>());
    reduce(minDet, minOp<scalar>());
    reduce(sumDet, sumOp<scalar>());
    reduce(nSummed, sumOp<label>());
 
    if (debug || report)
    {
        if (nSummed > 0)
        {
            Info<< "    Cell determinant (wellposedness) : minimum: " << minDet
                << " average: " << sumDet/nSummed
                << endl;
        }
    }
 
    if (nErrorCells > 0)
    {
        if (debug || report)
        {
            Info<< " ***Cells with small determinant found, number of cells: "
                << nErrorCells << endl;
        }
 
        return true;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Cell determinant check OK." << endl;
        }
 
        return false;
    }
 
    return false;
}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
bool Foam::primitiveMesh::checkTopology(const bool report) const
{
    label noFailedChecks = 0;
 
    if (checkPoints(report)) noFailedChecks++;
    if (checkUpperTriangular(report)) noFailedChecks++;
    if (checkCellsZipUp(report)) noFailedChecks++;
    if (checkFaceVertices(report)) noFailedChecks++;
    if (checkFaceFaces(report)) noFailedChecks++;
    //if (checkCellDeterminant(report)) noFailedChecks++;
 
    if (noFailedChecks == 0)
    {
        if (debug || report)
        {
            Info<< "    Mesh topology OK." << endl;
        }
 
        return false;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Failed " << noFailedChecks
                << " mesh topology checks." << endl;
        }
 
        return true;
    }
}
 
 
bool Foam::primitiveMesh::checkGeometry(const bool report) const
{
    label noFailedChecks = 0;
 
    if (checkClosedBoundary(report)) noFailedChecks++;
    if (checkClosedCells(report)) noFailedChecks++;
    if (checkFaceAreas(report)) noFailedChecks++;
    if (checkCellVolumes(report)) noFailedChecks++;
    if (checkFaceOrthogonality(report)) noFailedChecks++;
    if (checkFacePyramids(report)) noFailedChecks++;
    if (checkFaceSkewness(report)) noFailedChecks++;
 
    if (noFailedChecks == 0)
    {
        if (debug || report)
        {
            Info<< "    Mesh geometry OK." << endl;
        }
        return false;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Failed " << noFailedChecks
                << " mesh geometry checks." << endl;
        }
 
        return true;
    }
}
 
 
bool Foam::primitiveMesh::checkMesh(const bool report) const
{
    if (debug)
    {
        Info<< "bool primitiveMesh::checkMesh(const bool report) const: "
            << "checking primitiveMesh" << endl;
    }
 
    label noFailedChecks = checkTopology(report) + checkGeometry(report);
 
    if (noFailedChecks == 0)
    {
        if (debug || report)
        {
            Info<< "Mesh OK." << endl;
        }
 
        return false;
    }
    else
    {
        if (debug || report)
        {
            Info<< "    Failed " << noFailedChecks
                << " mesh checks." << endl;
        }
 
        return true;
    }
}
 
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 
// ************************************************************************* //