Skip to content
Snippets Groups Projects
polyMeshGenChecksGeometry.C 59.5 KiB
Newer Older
  • Learn to ignore specific revisions
  • franjo's avatar
    franjo committed
    /*---------------------------------------------------------------------------*\
      =========                 |
    
    Tomislav Lugaric's avatar
    Tomislav Lugaric committed
      \\      /  F ield         | cfMesh: A library for mesh generation
    
    franjo's avatar
    franjo committed
       \\    /   O peration     |
    
        \\  /    A nd           | Copyright held by the origina author
         \\/     M anipulation  |
    
    franjo's avatar
    franjo committed
    -------------------------------------------------------------------------------
    License
    
    Tomislav Lugaric's avatar
    Tomislav Lugaric committed
        This file is part of cfMesh.
    
    franjo's avatar
    franjo committed
    
    
    Tomislav Lugaric's avatar
    Tomislav Lugaric committed
        cfMesh is free software; you can redistribute it and/or modify it
    
    franjo's avatar
    franjo committed
        under the terms of the GNU General Public License as published by the
    
    Tomislav Lugaric's avatar
    Tomislav Lugaric committed
        Free Software Foundation; either version 3 of the License, or (at your
    
    franjo's avatar
    franjo committed
        option) any later version.
    
    
    Tomislav Lugaric's avatar
    Tomislav Lugaric committed
        cfMesh is distributed in the hope that it will be useful, but WITHOUT
    
    franjo's avatar
    franjo committed
        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
    
    Tomislav Lugaric's avatar
    Tomislav Lugaric committed
        along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
    
    franjo's avatar
    franjo committed
    
    \*---------------------------------------------------------------------------*/
    
    #include "polyMeshGenChecks.H"
    #include "polyMeshGenAddressing.H"
    #include "pyramidPointFaceRef.H"
    #include "tetrahedron.H"
    
    
    franjo's avatar
    franjo committed
    #include <omp.h>
    
    franjo's avatar
    franjo committed
    
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
    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();
    
    franjo's avatar
    franjo committed
    
        label nErrorClosed = 0;
    
    
    franjo's avatar
    franjo committed
        # pragma omp parallel for schedule(guided) reduction(+ : nErrorClosed)
    
    franjo's avatar
    franjo committed
        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 )
                {
    
    franjo's avatar
    franjo committed
                    # pragma omp critical
    
    franjo's avatar
    franjo committed
                    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;
    
    franjo's avatar
    franjo committed
            // 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: "
    
    franjo's avatar
    franjo committed
                << 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();
    
    franjo's avatar
    franjo committed
    
        const labelList& own = mesh.owner();
        const labelList& nei = mesh.neighbour();
    
        scalar minArea = VGREAT;
        scalar maxArea = -VGREAT;
    
    
    franjo's avatar
    franjo committed
        # pragma omp parallel if( own.size() > 100 )
    
    franjo's avatar
    franjo committed
        {
            scalar localMaxArea(-VGREAT), localMinArea(VGREAT);
    
    franjo's avatar
    franjo committed
            # pragma omp for schedule(guided)
    
            forAll(faceAreas, faceI)
    
    franjo's avatar
    franjo committed
            {
                if( changedFacePtr && !changedFacePtr->operator[](faceI) )
                    continue;
    
                const scalar magFaceArea = mag(faceAreas[faceI]);
    
                if( magFaceArea < minFaceArea )
    
    franjo's avatar
    franjo committed
                {
                    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;
    
    franjo's avatar
    franjo committed
                        }
                        else
                        {
                            Pout<< "Zero or negative face area detected for "
                                << "boundary face " << faceI << " next to cell "
                                << own[faceI] << ".  Face area magnitude = "
    
                                << magFaceArea << endl;
    
    franjo's avatar
    franjo committed
                        }
                    }
    
    franjo's avatar
    franjo committed
                    if( setPtr )
                    {
    
    franjo's avatar
    franjo committed
                        # pragma omp critical
    
    franjo's avatar
    franjo committed
                        setPtr->insert(faceI);
                    }
                }
    
                localMinArea = Foam::min(localMinArea, magFaceArea);
                localMaxArea = Foam::max(localMaxArea, magFaceArea);
    
    franjo's avatar
    franjo committed
            }
    
    franjo's avatar
    franjo committed
            # pragma omp critical
    
    franjo's avatar
    franjo committed
            {
                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: "
    
    franjo's avatar
    franjo committed
                << 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();
    
    franjo's avatar
    franjo committed
        const vectorField& fCentres = mesh.addressingData().faceCentres();
        const vectorField& cCentres = mesh.addressingData().cellCentres();
    
        label nNegVolCells = 0;
    
    
    franjo's avatar
    franjo committed
        # pragma omp parallel for if( owner.size() > 100 ) \
        schedule(guided) reduction(+ : nNegVolCells)
    
    franjo's avatar
    franjo committed
        forAll(owner, faceI)
        {
            if( changedFacePtr && !(*changedFacePtr)[faceI] )
                continue;
    
    franjo's avatar
    franjo committed
            const face& f = faces[faceI];
    
    franjo's avatar
    franjo committed
            bool badFace(false);
    
    franjo's avatar
    franjo committed
            forAll(f, eI)
            {
                const tetrahedron<point, point> tetOwn
                (
                    fCentres[faceI],
                    points[f.nextLabel(eI)],
                    points[f[eI]],
                    cCentres[owner[faceI]]
                );
    
    franjo's avatar
    franjo committed
                if( tetOwn.mag() < minPartTet )
                {
                    if( report )
                    {
    
    franjo's avatar
    franjo committed
                        # pragma omp critical
    
    franjo's avatar
    franjo committed
                        Pout<< "Zero or negative cell volume detected for cell "
                            << owner[faceI] << "." << endl;
                    }
    
    franjo's avatar
    franjo committed
                    badFace = true;
                }
    
    franjo's avatar
    franjo committed
                if( neighbour[faceI] < 0 )
                    continue;
    
    franjo's avatar
    franjo committed
                const tetrahedron<point, point> tetNei
                (
                    fCentres[faceI],
                    points[f[eI]],
                    points[f.nextLabel(eI)],
                    cCentres[neighbour[faceI]]
                );
    
    franjo's avatar
    franjo committed
                if( tetNei.mag() < minPartTet )
                {
                    if( report )
                    {
    
    franjo's avatar
    franjo committed
                        # pragma omp critical
    
    franjo's avatar
    franjo committed
                        Pout<< "Zero or negative cell volume detected for cell "
                            << neighbour[faceI] << "." << endl;
                    }
    
    franjo's avatar
    franjo committed
                    badFace = true;
                }
            }
    
    franjo's avatar
    franjo committed
            if( badFace )
            {
                if( setPtr )
                {
    
    franjo's avatar
    franjo committed
                    # pragma omp critical
    
    franjo's avatar
    franjo committed
                    setPtr->insert(faceI);
                }
    
    franjo's avatar
    franjo committed
                ++nNegVolCells;
            }
        }
    
    franjo's avatar
    franjo committed
        if( setPtr )
        {
            //- ensure that faces are selected at both sides
    
            const PtrList<processorBoundaryPatch>& procBnd = mesh.procBoundaries();
    
    franjo's avatar
    franjo committed
            forAll(procBnd, patchI)
            {
                const label start = procBnd[patchI].patchStart();
                const label size = procBnd[patchI].patchSize();
    
    franjo's avatar
    franjo committed
                for(label faceI=0;faceI<size;++faceI)
                {
                    if( setPtr->found(faceI+start) )
                        sendData.append(faceI);
                }
    
    franjo's avatar
    franjo committed
                OPstream toOtherProc
                (
                    Pstream::blocking,
                    procBnd[patchI].neiProcNo(),
                    sendData.byteSize()
                );
    
    franjo's avatar
    franjo committed
                toOtherProc << sendData;
            }
    
    franjo's avatar
    franjo committed
            forAll(procBnd, patchI)
            {
                labelList receivedData;
    
    franjo's avatar
    franjo committed
                IPstream fromOtherProc
                (
                    Pstream::blocking,
                    procBnd[patchI].neiProcNo()
                );
    
    franjo's avatar
    franjo committed
                fromOtherProc >> receivedData;
    
    franjo's avatar
    franjo committed
                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);
                    }
                }
            }
        }
    }
    
    
    franjo's avatar
    franjo committed
    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);
    
    franjo's avatar
    franjo committed
    
        const labelList& own = mesh.owner();
        const labelList& nei = mesh.neighbour();
    
        const label nInternalFaces = mesh.nInternalFaces();
    
    franjo's avatar
    franjo committed
    
        // Severe nonorthogonality threshold
        const scalar severeNonorthogonalityThreshold =
    
            ::cos(nonOrthWarn/180.0*M_PI);
    
    franjo's avatar
    franjo committed
    
        scalar minDDotS = VGREAT;
        scalar sumDDotS = 0.0;
    
        label severeNonOrth = 0;
        label errorNonOrth = 0;
        label counter = 0;
    
    
    franjo's avatar
    franjo committed
        # pragma omp parallel if( nInternalFaces > 1000 ) \
        reduction(+ : severeNonOrth, errorNonOrth, sumDDotS)
    
    franjo's avatar
    franjo committed
        {
            scalar localMinDDotS(VGREAT);
    
    franjo's avatar
    franjo committed
            # pragma omp for schedule(guided)
    
    franjo's avatar
    franjo committed
            for(label faceI=0;faceI<nInternalFaces;++faceI)
            {
    
                if( changedFacePtr && !(*changedFacePtr)[faceI] )
                    continue;
    
    
                const scalar dDotS = faceDotProduct[faceI];
    
    franjo's avatar
    franjo committed
                if( dDotS < severeNonorthogonalityThreshold )
                {
                    if( dDotS > SMALL )
                    {
                        if( report )
                        {
                            // Severe non-orthogonality but mesh still OK
    
    franjo's avatar
    franjo committed
                            # pragma omp critical
    
    franjo's avatar
    franjo committed
                            Pout<< "Severe non-orthogonality for face " << faceI
                                << " between cells " << own[faceI]
                                << " and " << nei[faceI]
                                << ": Angle = "
    
                                << ::acos(dDotS)/M_PI*180.0
    
    franjo's avatar
    franjo committed
                                << " deg." << endl;
                        }
    
    franjo's avatar
    franjo committed
                        if( setPtr )
                        {
    
    franjo's avatar
    franjo committed
                            # pragma omp critical
    
    franjo's avatar
    franjo committed
                            setPtr->insert(faceI);
                        }
    
    franjo's avatar
    franjo committed
                        ++severeNonOrth;
                    }
                    else
                    {
                        ++errorNonOrth;
    
    franjo's avatar
    franjo committed
                        if( setPtr )
                        {
    
    franjo's avatar
    franjo committed
                            # pragma omp critical
    
    franjo's avatar
    franjo committed
                            setPtr->insert(faceI);
                        }
                    }
                }
    
    franjo's avatar
    franjo committed
                localMinDDotS = Foam::min(dDotS, localMinDDotS);
                sumDDotS += dDotS;
            }
    
    franjo's avatar
    franjo committed
            # pragma omp critical
    
    franjo's avatar
    franjo committed
            minDDotS = Foam::min(minDDotS, localMinDDotS);
        }
    
    franjo's avatar
    franjo committed
        counter += nInternalFaces;
    
    franjo's avatar
    franjo committed
        if( Pstream::parRun() )
        {
    
            const PtrList<processorBoundaryPatch>& procBoundaries =
    
    franjo's avatar
    franjo committed
                mesh.procBoundaries();
    
            const label start = procBoundaries[0].patchStart();
    
            # ifdef USE_OMP
            # pragma omp parallel \
            reduction(+ : severeNonOrth, errorNonOrth, sumDDotS, counter)
            # endif
    
    franjo's avatar
    franjo committed
            {
    
                scalar localMinDDotS(VGREAT);
    
                # pragma omp for schedule(dynamic, 100)
    
                for(label faceI=start;faceI<own.size();++faceI)
    
    franjo's avatar
    franjo committed
                {
    
                    if( changedFacePtr && !(*changedFacePtr)[start+faceI] )
                                continue;
    
                    const scalar dDotS = faceDotProduct[faceI];
    
                    if( dDotS < severeNonorthogonalityThreshold )
                    {
                        if( dDotS > SMALL )
    
    franjo's avatar
    franjo committed
                        {
    
    franjo's avatar
    franjo committed
                            {
    
                                // Severe non-orthogonality but mesh still OK
                                # ifdef USE_OMP
                                # pragma omp critical
                                # endif
    
    franjo's avatar
    franjo committed
                                {
    
                                    const scalar angle
                                    (
                                        Foam::acos(dDotS) /
                                        M_PI * 180.0
                                    );
                                    Pout<< "Severe non-orthogonality for face "
                                        << start+faceI
                                        << ": Angle = "
                                        << angle
                                        << " deg." << endl;
    
    franjo's avatar
    franjo committed
                                }
                            }
    
                            if( setPtr )
                            {
                                # ifdef USE_OMP
                                # pragma omp critical
                                # endif
                                setPtr->insert(start+faceI);
    
    franjo's avatar
    franjo committed
                            }
    
    franjo's avatar
    franjo committed
                        }
    
                            if( setPtr )
                            {
                                # ifdef USE_OMP
                                # pragma omp critical
                                # endif
                                setPtr->insert(start+faceI);
                            }
                        }
    
    franjo's avatar
    franjo committed
                    }
    
                    localMinDDotS = Foam::min(dDotS, localMinDDotS);
                    sumDDotS += 0.5 * dDotS;
    
                    ++counter;
    
    franjo's avatar
    franjo committed
                }
    
                # ifdef USE_OMP
                # pragma omp critical
                # endif
                minDDotS = Foam::min(minDDotS, localMinDDotS);
    
    franjo's avatar
    franjo committed
            }
        }
    
    franjo's avatar
    franjo committed
        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
    
    franjo's avatar
    franjo committed
                    << " average: " <<
    
                       ::acos(sumDDotS/counter)/M_PI*180.0
    
    franjo's avatar
    franjo committed
                    << 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;
    
    
    franjo's avatar
    franjo committed
        # pragma omp parallel for schedule(guided) reduction(+ : nErrorPyrs)
    
    franjo's avatar
    franjo committed
        forAll(faces, faceI)
        {
    
            if( changedFacePtr && !(*changedFacePtr)[faceI] )
                continue;
    
    franjo's avatar
    franjo committed
            // Create the owner pyramid - it will have negative volume
            const scalar pyrVol = pyramidPointFaceRef
    
            (
                faces[faceI],
                ctrs[owner[faceI]]
            ).mag(points);
    
    franjo's avatar
    franjo committed
            bool badFace(false);
    
            if( pyrVol > -minPyrVol )
            {
                if( report )