polyMeshGenChecksGeometry.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | cfMesh: A library for mesh generation
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by the origina author
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of cfMesh.

    cfMesh 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.

    cfMesh 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 cfMesh.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "polyMeshGenChecks.H"
#include "polyMeshGenAddressing.H"
#include "pyramidPointFaceRef.H"
#include "tetrahedron.H"

# ifdef USE_OMP
#include <omp.h>
# endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace polyMeshGenChecks
{

bool checkClosedBoundary(const polyMeshGen& mesh, const bool report)
{
    // 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 = mesh.addressingData().faceAreas();

    for(label faceI = mesh.nInternalFaces();faceI<areas.size();++faceI)
    {
        sumClosed += areas[faceI];
        sumMagClosedBoundary += mag(areas[faceI]);
    }

    // Check the openness in the maximum direction (this is APPROXIMATE!)
    scalar maxOpen = max
    (
        sumClosed.component(vector::X),
        max
        (
            sumClosed.component(vector::Y),
            sumClosed.component(vector::Z)
        )
    );

    reduce(sumClosed, sumOp<vector>());
    reduce(maxOpen, maxOp<scalar>());

    if( maxOpen > SMALL*max(1.0, sumMagClosedBoundary) )
    {
        SeriousErrorIn
        (
            "bool checkClosedBoundary(const polyMeshGen&, const bool report)"
        )   << "Possible hole in boundary description" << endl;

        Info<< "Boundary openness in x-direction = "
            << sumClosed.component(vector::X) << endl;

        Info<< "Boundary openness in y-direction = "
            << sumClosed.component(vector::Y) << endl;

        Info<< "Boundary openness in z-direction = "
            << sumClosed.component(vector::Z) << endl;

        return true;
    }
    else
    {
        if( report )
        {
            Info<< "Boundary openness in x-direction = "
                << sumClosed.component(vector::X) << endl;

            Info<< "Boundary openness in y-direction = "
                << sumClosed.component(vector::Y) << endl;

            Info<< "Boundary openness in z-direction = "
                << sumClosed.component(vector::Z) << endl;

            Info<< "Boundary closed (OK)." << endl;
        }

        return false;
    }
}

bool checkClosedCells
(
    const polyMeshGen& mesh,
    const bool report,
    const scalar aspectWarn,
    labelHashSet* setPtr
)
{
    // Check that all cells labels are valid
    const cellListPMG& cells = mesh.cells();
    const label nFaces = mesh.faces().size();

    label nErrorClosed = 0;

    # ifdef USE_OMP
    # pragma omp parallel for schedule(guided) reduction(+ : nErrorClosed)
    # endif
    forAll(cells, cI)
    {
        const cell& curCell = cells[cI];

        if( min(curCell) < 0 || max(curCell) > nFaces )
        {
            WarningIn
            (
                "bool checkClosedCells("
                "const polyMeshGen&, const bool, const scalar, labelHashSet*)"
            )   << "Cell " << cI << " contains face labels out of range: "
                << curCell << " Max face index = " << nFaces << endl;

            if( setPtr )
            {
                # ifdef USE_OMP
                # pragma omp critical
                # endif
                setPtr->insert(cI);
            }

            ++nErrorClosed;
        }
    }

    if( nErrorClosed > 0 )
    {
        SeriousErrorIn
        (
            "bool checkClosedCells("
            "const polyMeshGen&, const bool, const scalar, labelHashSet*)"
        )  << nErrorClosed << " cells with invalid face labels found"
            << 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(cells.size(), vector::zero);

    scalarField sumMagClosed(cells.size(), 0.0);

    const labelList& own = mesh.owner();
    const labelList& nei = mesh.neighbour();
    const vectorField& areas = mesh.addressingData().faceAreas();

    forAll(own, faceI)
    {
        // Add to owner
        sumClosed[own[faceI]] += areas[faceI];
        sumMagClosed[own[faceI]] += mag(areas[faceI]);

        if( nei[faceI] == -1 )
            continue;

        // Subtract from neighbour
        sumClosed[nei[faceI]] -= areas[faceI];
        sumMagClosed[nei[faceI]] += mag(areas[faceI]);
    }

    label nOpen = 0;
    scalar maxOpenCell = 0;

    label nAspect = 0;
    scalar maxAspectRatio = 0;

    const scalarField& vols = mesh.addressingData().cellVolumes();

    // Check the sums
    forAll(sumClosed, cellI)
    {
        scalar maxOpen = max
        (
            sumClosed[cellI].component(vector::X),
            max
            (
                sumClosed[cellI].component(vector::Y),
                sumClosed[cellI].component(vector::Z)
            )
        );

        maxOpenCell = max(maxOpenCell, maxOpen);

        if( maxOpen > SMALL*max(1.0, sumMagClosed[cellI]) )
        {
            if( report )
            {
                Pout<< "Cell " << cellI << " is not closed. "
                    << "Face area vectors sum up to " << mag(sumClosed[cellI])
                    << " directionwise " << sumClosed[cellI] << " or "
                    << mag(sumClosed[cellI])
                       /(mag(sumMagClosed[cellI]) + VSMALL)
                    << " of the area of all faces of the cell. " << endl
                    << "    Area magnitudes sum up to "
                    << sumMagClosed[cellI] << endl;
            }

            if( setPtr )
            {
                setPtr->insert(cellI);
            }

            ++nOpen;
        }

        scalar aspectRatio =
            1.0/6.0*sumMagClosed[cellI]/pow(vols[cellI], 2.0/3.0);

        maxAspectRatio = max(maxAspectRatio, aspectRatio);

        if( aspectRatio > aspectWarn )
        {
            if( report )
            {
                Pout<< "High aspect ratio for cell " << cellI
                    << ": " << aspectRatio << endl;
            }

            ++nAspect;
        }
    }

    reduce(nOpen, sumOp<label>());
    reduce(maxOpenCell, maxOp<scalar>());

    reduce(nAspect, sumOp<label>());
    reduce(maxAspectRatio, maxOp<scalar>());

    if( nOpen > 0 )
    {
        SeriousErrorIn
        (
            "bool checkClosedCells("
            "const polyMeshGen&, const bool, const scalar, labelHashSet*)"
        )   << nOpen << " open cells found. Max cell openness: "
            << maxOpenCell << endl;

        return true;
    }

    if( nAspect > 0 )
    {
        SeriousErrorIn
        (
            "bool checkClosedCells("
            "const polyMeshGen&, const bool, const scalar, labelHashSet*)"
        )   << nAspect << " high aspect ratio cells found.  "
            << "Max aspect ratio: " << maxAspectRatio
            << endl;

        return true;
    }

    if( report )
    {
        Info<< "Max cell openness = " << maxOpenCell
            << "  Max aspect ratio = " << maxAspectRatio
            << ".  All cells OK.\n" << endl;
    }

    return false;
}

bool checkCellVolumes
(
    const polyMeshGen& mesh,
    const bool report,
    labelHashSet* setPtr
)
{
    const scalarField& vols = mesh.addressingData().cellVolumes();

    scalar minVolume = GREAT;
    scalar maxVolume = -GREAT;

    label nNegVolCells = 0;

    forAll(vols, cellI)
    {
        if( vols[cellI] < VSMALL )
        {
            if( report )
                SeriousErrorIn
                (
                    "bool checkCellVolumes("
                    "const polyMeshGen&, const bool, labelHashSet*)"
                )   << "Zero or negative cell volume detected for cell "
                    << cellI << ".  Volume = " << vols[cellI] << endl;

            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 )
    {
        SeriousErrorIn
        (
            "bool checkCellVolumes("
            "const polyMeshGen&, const bool, labelHashSet*)"
        )   << "Zero or negative cell volume detected.  "
            << "Minimum negative volume: "
            << minVolume << ".\nNumber of negative volume cells: "
            << nNegVolCells << ".  This mesh is invalid"
            << endl;

        return true;
    }
    else
    {
        if( report )
        {
            Info<< "Min volume = " << minVolume
                << ". Max volume = " << maxVolume
                << ".  Total volume = " << sum(vols)
                << ".  Cell volumes OK.\n" << endl;
        }

        return false;
    }
}

bool checkFaceAreas
(
    const polyMeshGen& mesh,
    const bool report,
    const scalar minFaceArea,
    labelHashSet* setPtr,
    const boolList* changedFacePtr
)
{
    const vectorField& faceAreas = mesh.addressingData().faceAreas();

    const labelList& own = mesh.owner();
    const labelList& nei = mesh.neighbour();

    scalar minArea = VGREAT;
    scalar maxArea = -VGREAT;

    # ifdef USE_OMP
    # pragma omp parallel if( own.size() > 100 )
    # endif
    {
        scalar localMaxArea(-VGREAT), localMinArea(VGREAT);

        # ifdef USE_OMP
        # pragma omp for schedule(guided)
        # endif
        forAll(faceAreas, faceI)
        {
            if( changedFacePtr && !changedFacePtr->operator[](faceI) )
                continue;

            const scalar magFaceArea = mag(faceAreas[faceI]);

            if( magFaceArea < minFaceArea )
            {
                if( report )
                {
                    if( nei[faceI] != -1 )
                    {
                        Pout<< "Zero or negative face area detected for "
                            << "internal face " << faceI << " between cells "
                            << own[faceI] << " and " << nei[faceI]
                            << ".  Face area magnitude = "
                            << magFaceArea << endl;
                    }
                    else
                    {
                        Pout<< "Zero or negative face area detected for "
                            << "boundary face " << faceI << " next to cell "
                            << own[faceI] << ".  Face area magnitude = "
                            << magFaceArea << endl;
                    }
                }

                if( setPtr )
                {
                    # ifdef USE_OMP
                    # pragma omp critical
                    # endif
                    setPtr->insert(faceI);
                }
            }

            localMinArea = Foam::min(localMinArea, magFaceArea);
            localMaxArea = Foam::max(localMaxArea, magFaceArea);
        }

        # ifdef USE_OMP
        # pragma omp critical
        # endif
        {
            minArea = Foam::min(minArea, localMinArea);
            maxArea = Foam::max(maxArea, localMaxArea);
        }
    }

    reduce(minArea, minOp<scalar>());
    reduce(maxArea, maxOp<scalar>());

    if( minArea < VSMALL )
    {
        SeriousErrorIn
        (
            "bool checkFaceAreas("
            "const polyMeshGen&, const bool, const scalar,"
            " labelHashSet*, const boolList*)"
        )   << "Zero or negative face area detected.  Minimum negative area: "
            << minArea << ". This mesh is invalid"
            << endl;

        return true;
    }
    else
    {
        if( report )
        {
            Info<< "Minumum face area = " << minArea
                << ". Maximum face area = " << maxArea
                << ".  Face area magnitudes OK.\n" << endl;
        }

        return false;
    }
}

bool checkCellPartTetrahedra
(
    const polyMeshGen& mesh,
    const bool report,
    const scalar minPartTet,
    labelHashSet* setPtr,
    const boolList* changedFacePtr
)
{
    const pointFieldPMG& points = mesh.points();
    const faceListPMG& faces = mesh.faces();
    const labelList& owner = mesh.owner();
    const labelList& neighbour = mesh.neighbour();

    const vectorField& fCentres = mesh.addressingData().faceCentres();
    const vectorField& cCentres = mesh.addressingData().cellCentres();

    label nNegVolCells = 0;

    # ifdef USE_OMP
    # pragma omp parallel for if( owner.size() > 100 ) \
    schedule(guided) reduction(+ : nNegVolCells)
    # endif
    forAll(owner, faceI)
    {
        if( changedFacePtr && !(*changedFacePtr)[faceI] )
            continue;

        const face& f = faces[faceI];

        bool badFace(false);

        forAll(f, eI)
        {
            const tetrahedron<point, point> tetOwn
            (
                fCentres[faceI],
                points[f.nextLabel(eI)],
                points[f[eI]],
                cCentres[owner[faceI]]
            );

            if( tetOwn.mag() < minPartTet )
            {
                if( report )
                {
                    # ifdef USE_OMP
                    # pragma omp critical
                    # endif
                    Pout<< "Zero or negative cell volume detected for cell "
                        << owner[faceI] << "." << endl;
                }

                badFace = true;
            }

            if( neighbour[faceI] < 0 )
                continue;

            const tetrahedron<point, point> tetNei
            (
                fCentres[faceI],
                points[f[eI]],
                points[f.nextLabel(eI)],
                cCentres[neighbour[faceI]]
            );

            if( tetNei.mag() < minPartTet )
            {
                if( report )
                {
                    # ifdef USE_OMP
                    # pragma omp critical
                    # endif
                    Pout<< "Zero or negative cell volume detected for cell "
                        << neighbour[faceI] << "." << endl;
                }

                badFace = true;
            }
        }

        if( badFace )
        {
            if( setPtr )
            {
                # ifdef USE_OMP
                # pragma omp critical
                # endif
                setPtr->insert(faceI);
            }

            ++nNegVolCells;
        }
    }

    if( setPtr )
    {
        //- ensure that faces are selected at both sides
        const PtrList<processorBoundaryPatch>& procBnd = mesh.procBoundaries();
        forAll(procBnd, patchI)
        {
            const label start = procBnd[patchI].patchStart();
            const label size = procBnd[patchI].patchSize();

            labelLongList sendData;
            for(label faceI=0;faceI<size;++faceI)
            {
                if( setPtr->found(faceI+start) )
                    sendData.append(faceI);
            }

            OPstream toOtherProc
            (
                Pstream::blocking,
                procBnd[patchI].neiProcNo(),
                sendData.byteSize()
            );

            toOtherProc << sendData;
        }

        forAll(procBnd, patchI)
        {
            labelList receivedData;

            IPstream fromOtherProc
            (
                Pstream::blocking,
                procBnd[patchI].neiProcNo()
            );

            fromOtherProc >> receivedData;

            const label start = procBnd[patchI].patchStart();
            forAll(receivedData, i)
                setPtr->insert(start+receivedData[i]);
        }
    }

    reduce(nNegVolCells, sumOp<label>());

    if( nNegVolCells != 0 )
    {
        WarningIn
        (
            "bool checkCellPartTetrahedra("
            "const polyMeshGen&, const bool, const scalar,"
            " labelHashSet*, const boolList*)"
        )   << nNegVolCells << " zero or negative part tetrahedra detected."
            << endl;

        return true;
    }
    else
    {
        if( report )
            Info<< "Part cell tetrahedra OK.\n" << endl;

        return false;
    }
}

void checkFaceDotProduct
(
    const polyMeshGen& mesh,
    scalarField& faceDotProduct,
    const boolList* changedFacePtr
)
{
    //- for all internal faces check theat the d dot S product is positive
    const vectorField& centres = mesh.addressingData().cellCentres();
    const vectorField& areas = mesh.addressingData().faceAreas();

    const labelList& own = mesh.owner();
    const labelList& nei = mesh.neighbour();
    const label nInternalFaces = mesh.nInternalFaces();

    faceDotProduct.setSize(own.size());
    faceDotProduct = 1.0;

    # ifdef USE_OMP
    # pragma omp parallel if( nInternalFaces > 1000 )
    # endif
    {
        # ifdef USE_OMP
        # pragma omp for schedule(dynamic, 100)
        # endif
        for(label faceI=0;faceI<nInternalFaces;++faceI)
        {
            if( changedFacePtr && !(*changedFacePtr)[faceI] )
                continue;

            const vector d = centres[nei[faceI]] - centres[own[faceI]];
            const vector& s = areas[faceI];

            faceDotProduct[faceI] = (d & s)/(mag(d)*mag(s) + VSMALL);
        }
    }

    if( Pstream::parRun() )
    {
        const PtrList<processorBoundaryPatch>& procBoundaries =
            mesh.procBoundaries();

        forAll(procBoundaries, patchI)
        {
            const label start = procBoundaries[patchI].patchStart();

            vectorField cCentres(procBoundaries[patchI].patchSize());
            forAll(cCentres, faceI)
                cCentres[faceI] = centres[own[start+faceI]];

            OPstream toOtherProc
            (
                Pstream::blocking,
                procBoundaries[patchI].neiProcNo(),
                cCentres.byteSize()
            );

            toOtherProc << cCentres;
        }

        forAll(procBoundaries, patchI)
        {
            vectorField otherCentres;
            IPstream fromOtherProc
            (
                Pstream::blocking,
                procBoundaries[patchI].neiProcNo()
            );

            fromOtherProc >> otherCentres;

            //- calculate skewness at processor faces
            const label start = procBoundaries[patchI].patchStart();
            # ifdef USE_OMP
            # pragma omp parallel
            # endif
            {
                # ifdef USE_OMP
                # pragma omp for schedule(dynamic, 100)
                # endif
                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 d = cNei - cOwn;
                    const vector& s = areas[faceI];

                    faceDotProduct[faceI] = (d & s)/(mag(d)*mag(s) + VSMALL);
                }
            }
        }
    }
}

bool checkFaceDotProduct
(
    const polyMeshGen& mesh,
    const bool report,
    const scalar nonOrthWarn,
    labelHashSet* setPtr,
    const boolList* changedFacePtr
)
{
    //- for all internal faces check theat the d dot S product is positive
    scalarField faceDotProduct;
    checkFaceDotProduct(mesh, faceDotProduct, changedFacePtr);

    const labelList& own = mesh.owner();
    const labelList& nei = mesh.neighbour();
    const label nInternalFaces = mesh.nInternalFaces();

    // Severe nonorthogonality threshold
    const scalar severeNonorthogonalityThreshold =
        ::cos(nonOrthWarn/180.0*M_PI);

    scalar minDDotS = VGREAT;
    scalar sumDDotS = 0.0;

    label severeNonOrth = 0;
    label errorNonOrth = 0;
    label counter = 0;

    # ifdef USE_OMP
    # pragma omp parallel if( nInternalFaces > 1000 ) \
    reduction(+ : severeNonOrth, errorNonOrth, sumDDotS)
    # endif
    {
        scalar localMinDDotS(VGREAT);
        # ifdef USE_OMP
        # pragma omp for schedule(guided)
        # endif
        for(label faceI=0;faceI<nInternalFaces;++faceI)
        {
            if( changedFacePtr && !(*changedFacePtr)[faceI] )
                continue;

            const scalar dDotS = faceDotProduct[faceI];

            if( dDotS < severeNonorthogonalityThreshold )
            {
                if( dDotS > SMALL )
                {
                    if( report )
                    {
                        // Severe non-orthogonality but mesh still OK
                        # ifdef USE_OMP
                        # pragma omp critical
                        # endif
                        Pout<< "Severe non-orthogonality for face " << faceI
                            << " between cells " << own[faceI]
                            << " and " << nei[faceI]
                            << ": Angle = "
                            << ::acos(dDotS)/M_PI*180.0
                            << " deg." << endl;
                    }

                    if( setPtr )
                    {
                        # ifdef USE_OMP
                        # pragma omp critical
                        # endif
                        setPtr->insert(faceI);
                    }

                    ++severeNonOrth;
                }
                else
                {
                    ++errorNonOrth;

                    if( setPtr )
                    {
                        # ifdef USE_OMP
                        # pragma omp critical
                        # endif
                        setPtr->insert(faceI);
                    }
                }
            }

            localMinDDotS = Foam::min(dDotS, localMinDDotS);
            sumDDotS += dDotS;
        }

        # ifdef USE_OMP
        # pragma omp critical
        # endif
        minDDotS = Foam::min(minDDotS, localMinDDotS);
    }

    counter += nInternalFaces;

    if( Pstream::parRun() )
    {
        const PtrList<processorBoundaryPatch>& procBoundaries =
            mesh.procBoundaries();

        const label start = procBoundaries[0].patchStart();

        # ifdef USE_OMP
        # pragma omp parallel \
        reduction(+ : severeNonOrth, errorNonOrth, sumDDotS, counter)
        # endif
        {
            scalar localMinDDotS(VGREAT);

            # ifdef USE_OMP
            # pragma omp for schedule(dynamic, 100)
            # endif
            for(label faceI=start;faceI<own.size();++faceI)
            {
                if( changedFacePtr && !(*changedFacePtr)[start+faceI] )
                            continue;

                const scalar dDotS = faceDotProduct[faceI];

                if( dDotS < severeNonorthogonalityThreshold )
                {
                    if( dDotS > SMALL )
                    {
                        if( report )
                        {
                            // Severe non-orthogonality but mesh still OK
                            # ifdef USE_OMP
                            # pragma omp critical
                            # endif
                            {
                                const scalar angle
                                (
                                    Foam::acos(dDotS) /
                                    M_PI * 180.0
                                );
                                Pout<< "Severe non-orthogonality for face "
                                    << start+faceI
                                    << ": Angle = "
                                    << angle
                                    << " deg." << endl;
                            }
                        }

                        if( setPtr )
                        {
                            # ifdef USE_OMP
                            # pragma omp critical
                            # endif
                            setPtr->insert(start+faceI);
                        }

                        ++severeNonOrth;
                    }
                    else
                    {
                        ++errorNonOrth;

                        if( setPtr )
                        {
                            # ifdef USE_OMP
                            # pragma omp critical
                            # endif
                            setPtr->insert(start+faceI);
                        }
                    }
                }

                localMinDDotS = Foam::min(dDotS, localMinDDotS);
                sumDDotS += 0.5 * dDotS;

                ++counter;
            }

            # ifdef USE_OMP
            # pragma omp critical
            # endif
            minDDotS = Foam::min(minDDotS, localMinDDotS);
        }
    }

    reduce(minDDotS, minOp<scalar>());
    reduce(sumDDotS, sumOp<scalar>());
    reduce(severeNonOrth, sumOp<label>());
    reduce(errorNonOrth, sumOp<label>());

    reduce(counter, sumOp<label>());

    // Only report if there are some internal faces
    if( counter > 0 )
    {
        if( minDDotS < severeNonorthogonalityThreshold )
        {
            Info<< "Number of non-orthogonality errors: " << errorNonOrth
                << ". Number of severely non-orthogonal faces: "
                << severeNonOrth  << "." << endl;
        }
    }

    if( report )
    {
        if( counter > 0 )
        {
            Info<< "Mesh non-orthogonality Max: "
                << ::acos(minDDotS)/M_PI*180.0
                << " average: " <<
                   ::acos(sumDDotS/counter)/M_PI*180.0
                << endl;
        }
    }

    if( errorNonOrth > 0 )
    {
        WarningIn
        (
            "checkFaceDotProduct("
            "const polyMeshGen&, const bool, const scalar,"
            " labelHashSet*, const boolList*)"
        )   << "Error in non-orthogonality detected" << endl;

        return true;
    }
    else
    {
        if( report )
            Info<< "Non-orthogonality check OK.\n" << endl;

        return false;
    }
}

bool checkFacePyramids
(
    const polyMeshGen& mesh,
    const bool report,
    const scalar minPyrVol,
    labelHashSet* setPtr,
    const boolList* changedFacePtr
)
{
    // check whether face area vector points to the cell with higher label
    const vectorField& ctrs = mesh.addressingData().cellCentres();

    const labelList& owner = mesh.owner();
    const labelList& neighbour = mesh.neighbour();

    const faceListPMG& faces = mesh.faces();
    const pointFieldPMG& points = mesh.points();

    label nErrorPyrs = 0;

    # ifdef USE_OMP
    # pragma omp parallel for schedule(guided) reduction(+ : nErrorPyrs)
    # endif
    forAll(faces, faceI)
    {
        if( changedFacePtr && !(*changedFacePtr)[faceI] )
            continue;

        // Create the owner pyramid - it will have negative volume
        const scalar pyrVol = pyramidPointFaceRef
        (
            faces[faceI],
            ctrs[owner[faceI]]
        ).mag(points);

        bool badFace(false);

        if( pyrVol > -minPyrVol )
        {
            if( report )
            {
                # ifdef USE_OMP
                # pragma omp critical
                # endif
                Pout<< "bool checkFacePyramids("
                    << "const bool, const scalar, labelHashSet*) : "
                    << "face " << faceI << " points the wrong way. " << endl
                    << "Pyramid volume: " << -pyrVol
                    << " Face " << faces[faceI] << " area: "
                    << faces[faceI].mag(points)
                    << " Owner cell: " << owner[faceI] << endl
                    << "Owner cell vertex labels: "
                    << mesh.cells()[owner[faceI]].labels(faces)
                    << endl;
            }

            badFace = true;
        }

        if( neighbour[faceI] != -1 )
        {
            // Create the neighbour pyramid - it will have positive volume
            const scalar pyrVol =
                pyramidPointFaceRef
                (
                    faces[faceI],
                    ctrs[neighbour[faceI]]
                ).mag(points);

            if( pyrVol < minPyrVol )
            {
                if( report )
                {
                    # ifdef USE_OMP
                    # pragma omp critical
                    # endif
                    Pout<< "bool checkFacePyramids("
                        << "const bool, const scalar, labelHashSet*) : "
                        << "face " << faceI << " points the wrong way. " << endl
                        << "Pyramid volume: " << -pyrVol
                        << " Face " << faces[faceI] << " area: "
                        << faces[faceI].mag(points)
                        << " Neighbour cell: " << neighbour[faceI] << endl
                        << "Neighbour cell vertex labels: "
                        << mesh.cells()[neighbour[faceI]].labels(faces)
                        << endl;
                }

                badFace = true;
            }
        }

        if( badFace )
        {
            if( setPtr )
            {
                # ifdef USE_OMP
                # pragma omp critical
                # endif
                setPtr->insert(faceI);
            }

            ++nErrorPyrs;
        }
    }

    reduce(nErrorPyrs, sumOp<label>());

    if( 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( nErrorPyrs > 0 )
    {
        if( Pstream::master() )
            WarningIn
            (
                "bool checkFacePyramids("
                "const polyMeshGen&, const bool, const scalar,"
                " labelHashSet*, const boolList*)"
            )   << "Error in face pyramids: " << nErrorPyrs
                << " faces pointing the wrong way!"
                << endl;

        return true;
    }
    else
    {
        if( report )
            Info<< "Face pyramids OK.\n" << endl;

        return false;
    }
}

void checkFaceSkewness
(
    const polyMeshGen& mesh,
    scalarField& faceSkewness,
    const boolList* changedFacePtr
)
{
    //- Warn if the skew correction vector is more than skewWarning times
    //- larger than the face area vector
    const labelList& own = mesh.owner();
    const labelList& nei = mesh.neighbour();
    const label nInternalFaces = mesh.nInternalFaces();

    const vectorField& centres = mesh.addressingData().cellCentres();
    const vectorField& fCentres = mesh.addressingData().faceCentres();

    faceSkewness.setSize(own.size());
    faceSkewness = 0.0;

    //- check internal faces
    # ifdef USE_OMP
    # pragma omp parallel for schedule(dynamic, 100)
    # endif
    for(label faceI=0;faceI<nInternalFaces;++faceI)
    {
        if( changedFacePtr && !changedFacePtr->operator[](faceI) )
            continue;

        const scalar dOwn = mag(fCentres[faceI] - centres[own[faceI]]);
        const scalar dNei = mag(fCentres[faceI] - centres[nei[faceI]]);

        const point faceIntersection =
            centres[own[faceI]]*dNei/(dOwn+dNei)
          + centres[nei[faceI]]*dOwn/(dOwn+dNei);

        faceSkewness[faceI] =
            mag(fCentres[faceI] - faceIntersection)
            /(mag(centres[nei[faceI]] - centres[own[faceI]]) + VSMALL);
    }

    if( Pstream::parRun() )
    {
        //- check parallel boundaries
        const PtrList<processorBoundaryPatch>& procBoundaries =
            mesh.procBoundaries();

        forAll(procBoundaries, patchI)
        {
            const label start = procBoundaries[patchI].patchStart();

            vectorField cCentres(procBoundaries[patchI].patchSize());
            forAll(cCentres, faceI)
                cCentres[faceI] = centres[own[start+faceI]];

            OPstream toOtherProc
            (
                Pstream::blocking,
                procBoundaries[patchI].neiProcNo(),
                cCentres.byteSize()
            );

            toOtherProc << cCentres;
        }

        forAll(procBoundaries, patchI)
        {
            const label start = procBoundaries[patchI].patchStart();

            vectorField otherCentres;
            IPstream fromOtheProc
            (
                Pstream::blocking,
                procBoundaries[patchI].neiProcNo()
            );

            fromOtheProc >> otherCentres;

            //- calculate skewness at processor faces
            # ifdef USE_OMP
            # pragma omp parallel
            # endif
            {
                # ifdef USE_OMP
                # pragma omp for schedule(dynamic, 100)
                # endif
                forAll(otherCentres, fI)
                {
                    const label faceI = start + fI;