Skip to content
Snippets Groups Projects
primitiveMeshCheck.C 42.9 KiB
Newer Older
  • Learn to ignore specific revisions
  • /*---------------------------------------------------------------------------*\
      =========                 |
      \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       \\    /   O peration     |
    
    OpenFOAM bot's avatar
    OpenFOAM bot committed
        \\  /    A nd           | www.openfoam.com
    
         \\/     M anipulation  |
    
    -------------------------------------------------------------------------------
    
    OpenFOAM bot's avatar
    OpenFOAM bot committed
        Copyright (C) 2011-2016 OpenFOAM Foundation
    
        Copyright (C) 2019-2020 OpenCFD Ltd.
    
    -------------------------------------------------------------------------------
    License
        This file is part of OpenFOAM.
    
    
    mattijs's avatar
    mattijs committed
        OpenFOAM is free software: you can redistribute it and/or modify it
        under the terms of the GNU General Public License as published by
        the Free Software Foundation, either version 3 of the License, or
        (at your option) any later version.
    
    
        OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
        ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
        FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
        for more details.
    
        You should have received a copy of the GNU General Public License
    
    mattijs's avatar
    mattijs committed
        along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
    
    
    \*---------------------------------------------------------------------------*/
    
    #include "primitiveMesh.H"
    #include "pyramidPointFaceRef.H"
    #include "ListOps.H"
    
    #include "unitConversion.H"
    
    #include "SortableList.H"
    
    #include "edgeHashes.H"
    
    #include "primitiveMeshTools.H"
    
    
    // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
    
    
    Foam::scalar Foam::primitiveMesh::closedThreshold_  = 1.0e-6;
    Foam::scalar Foam::primitiveMesh::aspectThreshold_  = 1000;
    Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70;    // deg
    Foam::scalar Foam::primitiveMesh::skewThreshold_    = 4;
    
    mattijs's avatar
    mattijs committed
    Foam::scalar Foam::primitiveMesh::planarCosAngle_   = 1.0e-6;
    
    
    
    // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
    
    
    bool Foam::primitiveMesh::checkClosedBoundary
    (
        const vectorField& areas,
        const bool report,
    
        const bitSet& internalOrCoupledFaces
    
        DebugInFunction << "Checking if boundary is closed" << endl;
    
    
        // Loop through all boundary faces and sum up the face area vectors.
        // For a closed boundary, this should be zero in all vector components
    
    
        Vector<solveScalar> sumClosed(Zero);
        solveScalar sumMagClosedBoundary = 0;
    
        for (label facei = nInternalFaces(); facei < areas.size(); facei++)
    
            if (!internalOrCoupledFaces.size() || !internalOrCoupledFaces[facei])
    
                sumClosed += Vector<solveScalar>(areas[facei]);
    
                sumMagClosedBoundary += mag(areas[facei]);
    
        reduce(sumClosed, sumOp<Vector<solveScalar>>());
        reduce(sumMagClosedBoundary, sumOp<solveScalar>());
    
        Vector<solveScalar> openness = sumClosed/(sumMagClosedBoundary + VSMALL);
    
    
        if (cmptMax(cmptMag(openness)) > closedThreshold_)
        {
            if (debug || report)
            {
                Info<< " ***Boundary openness " << openness
                    << " possible hole in boundary description."
                    << endl;
            }
    
            return true;
        }
    
    
        if (debug || report)
        {
            Info<< "    Boundary openness " << openness << " OK."
                << endl;
    
    bool Foam::primitiveMesh::checkClosedCells
    
        const vectorField& faceAreas,
        const scalarField& cellVolumes,
    
        const bool report,
        labelHashSet* setPtr,
    
        labelHashSet* aspectSetPtr,
        const Vector<label>& meshD
    
        DebugInFunction << "Checking if cells are closed" << endl;
    
    
        // Check that all cells labels are valid
        const cellList& c = cells();
    
        label nErrorClosed = 0;
    
    
        {
            const cell& curCell = c[cI];
    
            if (min(curCell) < 0 || max(curCell) > nFaces())
            {
                if (setPtr)
                {
                    setPtr->insert(cI);
                }
    
                nErrorClosed++;
            }
        }
    
        if (nErrorClosed > 0)
        {
            if (debug || report)
            {
                Info<< " ***Cells with invalid face labels found, number of cells "
                    << nErrorClosed << endl;
            }
    
            return true;
        }
    
    
    
        scalarField openness;
        scalarField aspectRatio;
        primitiveMeshTools::cellClosedness
        (
            *this,
            meshD,
            faceAreas,
            cellVolumes,
            openness,
            aspectRatio
        );
    
        scalar maxOpennessCell = max(openness);
    
        label nAspect = 0;
    
        scalar maxAspectRatio = max(aspectRatio);
    
        forAll(openness, celli)
    
            if (openness[celli] > closedThreshold_)
    
                    setPtr->insert(celli);
    
            if (aspectRatio[celli] > aspectThreshold_)
    
                    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
                    << endl;
            }
    
            return true;
        }
    
        if (nAspect > 0)
        {
            if (debug || report)
            {
                Info<< " ***High aspect ratio cells found, Max aspect ratio: "
                    << maxAspectRatio
                    << ", number of cells " << nAspect
                    << endl;
            }
    
            return true;
        }
    
        if (debug || report)
        {
            Info<< "    Max cell openness = " << maxOpennessCell << " OK." << nl
                << "    Max aspect ratio = " << maxAspectRatio << " OK."
                << endl;
        }
    
        return false;
    }
    
    
    
    bool Foam::primitiveMesh::checkFaceAreas
    
        const vectorField& faceAreas,
    
        const bool report,
    
        const bool detailedReport,
    
        labelHashSet* setPtr
    ) const
    {
    
        DebugInFunction << "Checking face area magnitudes" << endl;
    
        const scalarField magFaceAreas(mag(faceAreas));
    
    
        scalar minArea = GREAT;
        scalar maxArea = -GREAT;
    
    
        forAll(magFaceAreas, facei)
    
            if (magFaceAreas[facei] < VSMALL)
    
                    setPtr->insert(facei);
    
                    if (isInternalFace(facei))
    
                    {
                        Pout<< "Zero or negative face area detected for "
    
                            << "internal face "<< facei << " between cells "
                            << faceOwner()[facei] << " and "
                            << faceNeighbour()[facei]
                            << ".  Face area magnitude = " << magFaceAreas[facei]
    
                            << endl;
                    }
                    else
                    {
                        Pout<< "Zero or negative face area detected for "
    
                            << "boundary face " << facei << " next to cell "
                            << faceOwner()[facei] << ".  Face area magnitude = "
                            << magFaceAreas[facei] << endl;
    
            minArea = min(minArea, magFaceAreas[facei]);
            maxArea = max(maxArea, magFaceAreas[facei]);
    
        }
    
        reduce(minArea, minOp<scalar>());
        reduce(maxArea, maxOp<scalar>());
    
        if (minArea < VSMALL)
        {
            if (debug || report)
            {
                Info<< " ***Zero or negative face area detected.  "
                    "Minimum area: " << minArea << endl;
            }
    
            return true;
        }
    
    
        if (debug || report)
        {
            Info<< "    Minimum face area = " << minArea
                << ". Maximum face area = " << maxArea
                << ".  Face area magnitudes OK." << endl;
    
    bool Foam::primitiveMesh::checkCellVolumes
    
        const scalarField& vols,
    
        const bool report,
    
        const bool detailedReport,
    
        labelHashSet* setPtr
    ) const
    {
    
        DebugInFunction << "Checking cell volumes" << endl;
    
    
        scalar minVolume = GREAT;
        scalar maxVolume = -GREAT;
    
        label nNegVolCells = 0;
    
    
            if (vols[celli] < VSMALL)
    
                    setPtr->insert(celli);
    
                if (detailedReport)
                {
                    Pout<< "Zero or negative cell volume detected for cell "
    
                        << celli << ".  Volume = " << vols[celli] << endl;
    
            minVolume = min(minVolume, vols[celli]);
            maxVolume = max(maxVolume, vols[celli]);
    
        }
    
        reduce(minVolume, minOp<scalar>());
        reduce(maxVolume, maxOp<scalar>());
        reduce(nNegVolCells, sumOp<label>());
    
        if (minVolume < VSMALL)
        {
            if (debug || report)
            {
                Info<< " ***Zero or negative cell volume detected.  "
                    << "Minimum negative volume: " << minVolume
                    << ", Number of negative volume cells: " << nNegVolCells
                    << endl;
            }
    
            return true;
        }
    
    
        if (debug || report)
        {
            Info<< "    Min volume = " << minVolume
                << ". Max volume = " << maxVolume
                << ".  Total volume = " << gSum(vols)
                << ".  Cell volumes OK." << endl;
    
    bool Foam::primitiveMesh::checkFaceOrthogonality
    
        const vectorField& fAreas,
        const vectorField& cellCtrs,
    
        const bool report,
        labelHashSet* setPtr
    ) const
    {
    
        DebugInFunction << "Checking mesh non-orthogonality" << endl;
    
        tmp<scalarField> tortho = primitiveMeshTools::faceOrthogonality
        (
            *this,
            fAreas,
            cellCtrs
        );
        const scalarField& ortho = tortho();
    
    
        // Severe nonorthogonality threshold
        const scalar severeNonorthogonalityThreshold =
    
            ::cos(degToRad(nonOrthThreshold_));
    
        scalar minDDotS = min(ortho);
    
        scalar sumDDotS = sum(ortho);
    
    
        label severeNonOrth = 0;
    
        label errorNonOrth = 0;
    
    
    
        forAll(ortho, facei)
    
            if (ortho[facei] < severeNonorthogonalityThreshold)
    
                if (ortho[facei] > SMALL)
    
                        setPtr->insert(facei);
    
                        setPtr->insert(facei);
    
                    }
    
                    errorNonOrth++;
                }
            }
        }
    
        reduce(minDDotS, minOp<scalar>());
        reduce(sumDDotS, sumOp<scalar>());
        reduce(severeNonOrth, sumOp<label>());
        reduce(errorNonOrth, sumOp<label>());
    
        if (debug || report)
        {
    
            label neiSize = ortho.size();
    
            reduce(neiSize, sumOp<label>());
    
            if (neiSize > 0)
            {
                if (debug || report)
                {
                    Info<< "    Mesh non-orthogonality Max: "
    
                        << radToDeg(::acos(minDDotS))
                        << " average: " << radToDeg(::acos(sumDDotS/neiSize))
    
                        << endl;
                }
            }
    
            if (severeNonOrth > 0)
            {
                Info<< "   *Number of severely non-orthogonal faces: "
                    << severeNonOrth << "." << endl;
            }
        }
    
        if (errorNonOrth > 0)
        {
            if (debug || report)
            {
                Info<< " ***Number of non-orthogonality errors: "
                    << errorNonOrth << "." << endl;
            }
    
            return true;
        }
    
    
        if (debug || report)
        {
            Info<< "    Non-orthogonality check OK." << endl;
    
    bool Foam::primitiveMesh::checkFacePyramids
    
        const pointField& points,
        const vectorField& ctrs,
    
        const bool report,
    
        const bool detailedReport,
    
        const scalar minPyrVol,
        labelHashSet* setPtr
    ) const
    {
    
        DebugInFunction << "Checking face orientation" << endl;
    
    
        const labelList& own = faceOwner();
        const labelList& nei = faceNeighbour();
        const faceList& f = faces();
    
    
    
        scalarField ownPyrVol;
        scalarField neiPyrVol;
        primitiveMeshTools::facePyramidVolume
        (
            *this,
            points,
            ctrs,
            ownPyrVol,
            neiPyrVol
        );
    
    
        forAll(ownPyrVol, facei)
    
            if (ownPyrVol[facei] < minPyrVol)
    
                    setPtr->insert(facei);
    
                    Pout<< "Negative pyramid volume: " << ownPyrVol[facei]
                        << " for face " << facei << " " << f[facei]
                        << "  and owner cell: " << own[facei] << endl
    
                        << "Owner cell vertex labels: "
    
                        << cells()[own[facei]].labels(faces())
    
            if (isInternalFace(facei))
    
                if (neiPyrVol[facei] < minPyrVol)
    
                        setPtr->insert(facei);
    
                        Pout<< "Negative pyramid volume: " << neiPyrVol[facei]
                            << " for face " << facei << " " << f[facei]
                            << "  and neighbour cell: " << nei[facei] << nl
    
                            << "Neighbour cell vertex labels: "
    
                            << cells()[nei[facei]].labels(faces())
    
                    nErrorPyrs++;
                }
            }
        }
    
        reduce(nErrorPyrs, sumOp<label>());
    
        if (nErrorPyrs > 0)
        {
            if (debug || report)
            {
                Info<< " ***Error in face pyramids: "
                    << nErrorPyrs << " faces are incorrectly oriented."
                    << endl;
            }
    
            return true;
        }
    
    
        if (debug || report)
        {
            Info<< "    Face pyramids OK." << endl;
    
    bool Foam::primitiveMesh::checkFaceSkewness
    
        const pointField& points,
        const vectorField& fCtrs,
        const vectorField& fAreas,
        const vectorField& cellCtrs,
    
        const bool report,
        labelHashSet* setPtr
    ) const
    {
    
        DebugInFunction << "Checking face skewness" << endl;
    
    
        // Warn if the skew correction vector is more than skewWarning times
        // larger than the face area vector
    
    
        tmp<scalarField> tskewness = primitiveMeshTools::faceSkewness
        (
            *this,
            points,
            fCtrs,
            fAreas,
            cellCtrs
        );
        const scalarField& skewness = tskewness();
    
        scalar maxSkew = max(skewness);
    
        forAll(skewness, facei)
    
        {
            // Check if the skewness vector is greater than the PN vector.
            // This does not cause trouble but is a good indication of a poor mesh.
    
            if (skewness[facei] > skewThreshold_)
    
                    setPtr->insert(facei);
    
                }
    
                nWarnSkew++;
            }
        }
    
        reduce(maxSkew, maxOp<scalar>());
        reduce(nWarnSkew, sumOp<label>());
    
        if (nWarnSkew > 0)
        {
            if (debug || report)
            {
                Info<< " ***Max skewness = " << maxSkew
                    << ", " << nWarnSkew << " highly skew faces detected"
                       " which may impair the quality of the results"
                    << endl;
            }
    
            return true;
        }
    
    
        if (debug || report)
        {
            Info<< "    Max skewness = " << maxSkew << " OK." << endl;
    
    bool Foam::primitiveMesh::checkFaceAngles
    
        const pointField& points,
        const vectorField& faceAreas,
    
        const bool report,
        const scalar maxDeg,
        labelHashSet* setPtr
    ) const
    {
    
        DebugInFunction << "Checking face angles" << endl;
    
    
        if (maxDeg < -SMALL || maxDeg > 180+SMALL)
        {
    
            FatalErrorInFunction
                << "maxDeg should be [0..180] but is now " << maxDeg
    
        const scalar maxSin = Foam::sin(degToRad(maxDeg));
    
        tmp<scalarField> tfaceAngles = primitiveMeshTools::faceConcavity
        (
            maxSin,
            *this,
            points,
            faceAreas
        );
        const scalarField& faceAngles = tfaceAngles();
    
        scalar maxEdgeSin = max(faceAngles);
    
        forAll(faceAngles, facei)
    
            if (faceAngles[facei] > SMALL)
    
                    setPtr->insert(facei);
    
                }
            }
        }
    
        reduce(nConcave, sumOp<label>());
        reduce(maxEdgeSin, maxOp<scalar>());
    
        if (nConcave > 0)
        {
            scalar maxConcaveDegr =
    
                radToDeg(Foam::asin(Foam::min(1.0, maxEdgeSin)));
    
    
            if (debug || report)
            {
                Info<< "   *There are " << nConcave
                    << " faces with concave angles between consecutive"
                    << " edges. Max concave angle = " << maxConcaveDegr
                    << " degrees." << endl;
            }
    
            return true;
        }
    
    
        if (debug || report)
        {
            Info<< "    All angles in faces OK." << endl;
    
    bool Foam::primitiveMesh::checkFaceFlatness
    
        const pointField& points,
        const vectorField& faceCentres,
        const vectorField& faceAreas,
    
        const bool report,
        const scalar warnFlatness,
        labelHashSet* setPtr
    ) const
    {
    
        DebugInFunction << "Checking face flatness" << endl;
    
    
        if (warnFlatness < 0 || warnFlatness > 1)
        {
    
                << "warnFlatness should be [0..1] but is " << warnFlatness
    
                << exit(FatalError);
        }
    
        const faceList& fcs = faces();
    
    
        tmp<scalarField> tfaceFlatness = primitiveMeshTools::faceFlatness
        (
            *this,
            points,
            faceCentres,
            faceAreas
        );
        const scalarField& faceFlatness = tfaceFlatness();
    
        scalarField magAreas(mag(faceAreas));
    
    
        scalar minFlatness = GREAT;
        scalar sumFlatness = 0;
        label nSummed = 0;
    
        forAll(faceFlatness, facei)
    
            if (fcs[facei].size() > 3 && magAreas[facei] > VSMALL)
    
                sumFlatness += faceFlatness[facei];
    
                minFlatness = min(minFlatness, faceFlatness[facei]);
    
                if (faceFlatness[facei] < warnFlatness)
    
                        setPtr->insert(facei);
    
                    }
                }
            }
        }
    
    
        reduce(nWarped, sumOp<label>());
        reduce(minFlatness, minOp<scalar>());
    
        reduce(nSummed, sumOp<label>());
        reduce(sumFlatness, sumOp<scalar>());
    
        if (debug || report)
        {
            if (nSummed > 0)
            {
    
                Info<< "    Face flatness (1 = flat, 0 = butterfly) : min = "
                    << minFlatness << "  average = " << sumFlatness / nSummed
                    << endl;
    
        {
            if (debug || report)
            {
                Info<< "   *There are " << nWarped
                    << " faces with ratio between projected and actual area < "
                    << warnFlatness << endl;
    
                Info<< "    Minimum ratio (minimum flatness, maximum warpage) = "
                    << minFlatness << endl;
            }
    
            return true;
        }
    
    
        if (debug || report)
        {
            Info<< "    All face flatness OK." << endl;
    
    bool Foam::primitiveMesh::checkConcaveCells
    
        const vectorField& fAreas,
        const pointField& fCentres,
    
        const bool report,
        labelHashSet* setPtr
    ) const
    {
    
        DebugInFunction << "Checking for concave cells" << endl;
    
        const cellList& c = cells();
        const labelList& fOwner = faceOwner();
    
        label nConcaveCells = 0;
    
            const cell& cFaces = c[celli];
    
                const point& fC = fCentres[fI];
    
                vector fN = fAreas[fI];
    
                fN /= max(mag(fN), VSMALL);
    
                // Flip normal if required so that it is always pointing out of
                // the cell
    
                if (fOwner[fI] != celli)
    
                // Is the centre of any other face of the cell on the
                // wrong side of the plane of this face?
    
                forAll(cFaces, j)
                {
                    if (j != i)
                    {
                        label fJ = cFaces[j];
    
                        const point& pt = fCentres[fJ];
    
                        // If the cell is concave, the point will be on the
                        // positive normal side of the plane of f, defined by
                        // its centre and normal, and the angle between (pt -
                        // fC) and fN will be less than 90 degrees, so the dot
                        // product will be positive.
    
                        vector pC = (pt - fC);
    
                        pC /= max(mag(pC), VSMALL);
    
                        if ((pC & fN) > -planarCosAngle_)
    
                            // Concave or planar face
    
                            concave = true;
    
                            if (setPtr)
    
                                setPtr->insert(celli);
    
        reduce(nConcaveCells, sumOp<label>());
    
        {
            if (debug || report)
            {
    
                Info<< " ***Concave cells (using face planes) found,"
                    << " number of cells: " << nConcaveCells << endl;
    
        if (debug || report)
        {
            Info<< "    Concave cell check OK." << endl;
    
    bool Foam::primitiveMesh::checkUpperTriangular
    
    (
        const bool report,
        labelHashSet* setPtr
    ) const
    {
    
        DebugInFunction << "Checking face ordering" << endl;
    
    
        // Check whether internal faces are ordered in the upper triangular order
        const labelList& own = faceOwner();
        const labelList& nei = faceNeighbour();
    
        const cellList& c = cells();
    
        label internal = nInternalFaces();
    
        // Has error occurred?
        bool error = false;
        // Have multiple faces been detected?
        label nMultipleCells = false;
    
        // Loop through faceCells once more and make sure that for internal cell
        // the first label is smaller
    
        for (label facei = 0; facei < internal; facei++)
    
            if (own[facei] >= nei[facei])
    
                    setPtr->insert(facei);
    
                }
            }
        }
    
        // Loop through all cells. For each cell, find the face that is internal
        // and add it to the check list (upper triangular order).
        // Once the list is completed, check it against the faceCell list
    
    
            const labelList& curFaces = c[celli];
    
    
            // Neighbouring cells
            SortableList<label> nbr(curFaces.size());
    
            forAll(curFaces, i)
            {
    
                label facei = curFaces[i];
    
                if (facei >= nInternalFaces())
    
                {
                    // Sort last
                    nbr[i] = labelMax;
                }
                else
                {
    
                    label nbrCelli = nei[facei];
    
                    }
                    else
                    {
                        // nbrCell is master. Let it handle this face.
                        nbr[i] = labelMax;
                    }
                }
            }
    
            nbr.sort();
    
            // Now nbr holds the cellCells in incremental order. Check:
            // - neighbouring cells appear only once. Since nbr is sorted this
            //   is simple check on consecutive elements
            // - faces indexed in same order as nbr are incrementing as well.
    
            label prevCell = nbr[0];