polyMeshFilter.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2012-2016 OpenFOAM Foundation
    Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

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

#include "polyMeshFilter.H"
#include "polyMesh.H"
#include "fvMesh.H"
#include "unitConversion.H"
#include "edgeCollapser.H"
#include "syncTools.H"
#include "polyTopoChange.H"
#include "globalIndex.H"
#include "bitSet.H"
#include "pointSet.H"
#include "faceSet.H"
#include "cellSet.H"

// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //

namespace Foam
{
defineTypeNameAndDebug(polyMeshFilter, 0);
}


void Foam::polyMeshFilter::updateSets(const mapPolyMesh& map)
{
    updateSets<pointSet>(map);
    updateSets<faceSet>(map);
    updateSets<cellSet>(map);
}


void Foam::polyMeshFilter::copySets
(
    const polyMesh& oldMesh,
    const polyMesh& newMesh
)
{
    copySets<pointSet>(oldMesh, newMesh);
    copySets<faceSet>(oldMesh, newMesh);
    copySets<cellSet>(oldMesh, newMesh);
}


Foam::autoPtr<Foam::fvMesh> Foam::polyMeshFilter::copyMesh(const fvMesh& mesh)
{
    polyTopoChange originalMeshToNewMesh(mesh);

    autoPtr<fvMesh> meshCopy;
    autoPtr<mapPolyMesh> mapPtr = originalMeshToNewMesh.makeMesh
    (
        meshCopy,
        IOobject
        (
            mesh.name(),
            mesh.polyMesh::instance(),
            mesh.time(),
            IOobject::READ_IF_PRESENT,  // read fv* if present
            IOobject::NO_WRITE,
            IOobject::NO_REGISTER
        ),
        mesh,
        true // parallel sync
    );

    const mapPolyMesh& map = mapPtr();

    // Update fields
    meshCopy().updateMesh(map);
    if (map.hasMotionPoints())
    {
        meshCopy().movePoints(map.preMotionPoints());
    }

    copySets(mesh, meshCopy());

    return meshCopy;
}


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

Foam::label Foam::polyMeshFilter::filterFacesLoop(const label nOriginalBadFaces)
{
    label nBadFaces = labelMax;
    label nOuterIterations = 0;

    // Maintain the number of times a point has been part of a bad face
    labelList pointErrorCount(mesh_.nPoints(), Zero);

    bitSet newErrorPoint(mesh_.nPoints());
    edgeCollapser::checkMeshQuality
    (
        mesh_,
        meshQualityCoeffDict(),
        newErrorPoint
    );

    bool newBadFaces = true;

    // Main loop
    // ~~~~~~~~~
    // It tries and do some collapses, checks the resulting mesh and
    // 'freezes' some edges (by marking in minEdgeLen) and tries again.
    // This will iterate ultimately to the situation where every edge is
    // frozen and nothing gets collapsed.
    while
    (
        nOuterIterations < maxIterations()
     //&& nBadFaces > nOriginalBadFaces
     && newBadFaces
    )
    {
        Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
            << endl;

        printScalarFieldStats("Edge Filter Factor", minEdgeLen_);
        printScalarFieldStats("Face Filter Factor", faceFilterFactor_);

        // Reset the new mesh to the old mesh
        newMeshPtr_ = copyMesh(mesh_);
        fvMesh& newMesh = newMeshPtr_();

        scalarField newMeshFaceFilterFactor = faceFilterFactor_;
        pointPriority_.reset(new labelList(originalPointPriority_));

        labelList origToCurrentPointMap(identity(newMesh.nPoints()));

        {
            label nInnerIterations = 0;
            label nPrevLocalCollapse = labelMax;

            Info<< incrIndent;

            while (true)
            {
                Info<< nl << indent << "Inner iteration = "
                    << nInnerIterations++ << nl << incrIndent << endl;

                label nLocalCollapse = filterFaces
                (
                    newMesh,
                    newMeshFaceFilterFactor,
                    origToCurrentPointMap
                );
                Info<< decrIndent;

                if
                (
                    nLocalCollapse >= nPrevLocalCollapse
                 || nLocalCollapse == 0
                )
                {
                    Info<< decrIndent;
                    break;
                }
                else
                {
                    nPrevLocalCollapse = nLocalCollapse;
                }
            }
        }

        scalarField newMeshMinEdgeLen = minEdgeLen_;

        {
            label nInnerIterations = 0;
            label nPrevLocalCollapse = labelMax;

            Info<< incrIndent;

            while (true)
            {
                Info<< nl << indent << "Inner iteration = "
                    << nInnerIterations++ << nl << incrIndent << endl;

                label nLocalCollapse = filterEdges
                (
                    newMesh,
                    newMeshMinEdgeLen,
                    origToCurrentPointMap
                );
                Info<< decrIndent;

                if
                (
                    nLocalCollapse >= nPrevLocalCollapse
                 || nLocalCollapse == 0
                )
                {
                    Info<< decrIndent;
                    break;
                }
                else
                {
                    nPrevLocalCollapse = nLocalCollapse;
                }
            }
        }


        // Mesh check
        // ~~~~~~~~~~~~~~~~~~
        // Do not allow collapses in regions of error.
        // Updates minEdgeLen, nRelaxedEdges

        if (controlMeshQuality())
        {
            bitSet isErrorPoint(newMesh.nPoints());
            nBadFaces = edgeCollapser::checkMeshQuality
            (
                newMesh,
                meshQualityCoeffDict(),
                isErrorPoint
            );

            Info<< nl << "    Number of bad faces     : " << nBadFaces << nl
                << "    Number of marked points : "
                << returnReduce(isErrorPoint.count(), sumOp<unsigned int>())
                << endl;

            updatePointErrorCount
            (
                isErrorPoint,
                origToCurrentPointMap,
                pointErrorCount
            );

            checkMeshEdgesAndRelaxEdges
            (
                newMesh,
                origToCurrentPointMap,
                isErrorPoint,
                pointErrorCount
            );

            checkMeshFacesAndRelaxEdges
            (
                newMesh,
                origToCurrentPointMap,
                isErrorPoint,
                pointErrorCount
            );

            newBadFaces = false;
            forAll(mesh_.points(), pI)
            {
                if
                (
                    origToCurrentPointMap[pI] >= 0
                 && isErrorPoint[origToCurrentPointMap[pI]]
                )
                {
                    if (!newErrorPoint[pI])
                    {
                        newBadFaces = true;
                        break;
                    }
                }
            }

            Pstream::reduceOr(newBadFaces);
        }
        else
        {
            return -1;
        }
    }

    return nBadFaces;
}


Foam::label Foam::polyMeshFilter::filterFaces
(
    polyMesh& newMesh,
    scalarField& newMeshFaceFilterFactor,
    labelList& origToCurrentPointMap
)
{
    // Per edge collapse status
    bitSet collapseEdge(newMesh.nEdges());

    Map<point> collapsePointToLocation(newMesh.nPoints());

    edgeCollapser collapser(newMesh, collapseFacesCoeffDict());

    {
        // Collapse faces
        labelPair nCollapsedPtEdge = collapser.markSmallSliverFaces
        (
            newMeshFaceFilterFactor,
            pointPriority_(),
            collapseEdge,
            collapsePointToLocation
        );

        label nCollapsed = 0;
        forAll(nCollapsedPtEdge, collapseTypeI)
        {
            nCollapsed += nCollapsedPtEdge[collapseTypeI];
        }

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

        label nToPoint = returnReduce(nCollapsedPtEdge.first(), sumOp<label>());
        label nToEdge = returnReduce(nCollapsedPtEdge.second(), sumOp<label>());

        Info<< indent
            << "Collapsing " << nCollapsed << " faces "
            << "(to point = " << nToPoint << ", to edge = " << nToEdge << ")"
            << endl;

        if (nCollapsed == 0)
        {
            return 0;
        }
    }

    // Merge edge collapses into consistent collapse-network.
    // Make sure no cells get collapsed.
    List<pointEdgeCollapse> allPointInfo;
    const globalIndex globalPoints(newMesh.nPoints());

    collapser.consistentCollapse
    (
        globalPoints,
        pointPriority_(),
        collapsePointToLocation,
        collapseEdge,
        allPointInfo
    );

    label nLocalCollapse = collapseEdge.count();

    reduce(nLocalCollapse, sumOp<label>());
    Info<< nl << indent << "Collapsing " << nLocalCollapse
        << " edges after synchronisation and PointEdgeWave" << endl;

    if (nLocalCollapse == 0)
    {
        return 0;
    }

    {
        // Apply collapses to current mesh
        polyTopoChange newMeshMod(newMesh);

        // Insert mesh refinement into polyTopoChange.
        collapser.setRefinement(allPointInfo, newMeshMod);

        Info<< indent << "Apply changes to the current mesh" << endl;

        // Apply changes to current mesh
        autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
        (
            newMesh,
            false
        );
        const mapPolyMesh& newMap = newMapPtr();

        // Update fields
        newMesh.updateMesh(newMap);
        if (newMap.hasMotionPoints())
        {
            newMesh.movePoints(newMap.preMotionPoints());
        }
        updateSets(newMap);

        updatePointPriorities(newMesh, newMap.pointMap());

        mapOldMeshFaceFieldToNewMesh
        (
            newMesh,
            newMap.faceMap(),
            newMeshFaceFilterFactor
        );

        updateOldToNewPointMap
        (
            newMap.reversePointMap(),
            origToCurrentPointMap
        );
    }

    return nLocalCollapse;
}


Foam::label Foam::polyMeshFilter::filterEdges
(
    polyMesh& newMesh,
    scalarField& newMeshMinEdgeLen,
    labelList& origToCurrentPointMap
)
{
    // Per edge collapse status
    bitSet collapseEdge(newMesh.nEdges());

    Map<point> collapsePointToLocation(newMesh.nPoints());

    edgeCollapser collapser(newMesh, collapseFacesCoeffDict());

    // Work out which edges to collapse
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // This is by looking at minEdgeLen (to avoid frozen edges)
    // and marking in collapseEdge.
    label nSmallCollapsed = collapser.markSmallEdges
    (
        newMeshMinEdgeLen,
        pointPriority_(),
        collapseEdge,
        collapsePointToLocation
    );

    reduce(nSmallCollapsed, sumOp<label>());
    Info<< indent << "Collapsing " << nSmallCollapsed
        << " small edges" << endl;

    // Merge inline edges
    label nMerged = collapser.markMergeEdges
    (
        maxCos(),
        pointPriority_(),
        collapseEdge,
        collapsePointToLocation
    );

    reduce(nMerged, sumOp<label>());
    Info<< indent << "Collapsing " << nMerged << " in line edges"
        << endl;

    if (nMerged + nSmallCollapsed == 0)
    {
        return 0;
    }

    // Merge edge collapses into consistent collapse-network.
    // Make sure no cells get collapsed.
    List<pointEdgeCollapse> allPointInfo;
    const globalIndex globalPoints(newMesh.nPoints());

    collapser.consistentCollapse
    (
        globalPoints,
        pointPriority_(),
        collapsePointToLocation,
        collapseEdge,
        allPointInfo
    );

    label nLocalCollapse = collapseEdge.count();

    reduce(nLocalCollapse, sumOp<label>());
    Info<< nl << indent << "Collapsing " << nLocalCollapse
        << " edges after synchronisation and PointEdgeWave" << endl;

    if (nLocalCollapse == 0)
    {
        return 0;
    }

    // Apply collapses to current mesh
    polyTopoChange newMeshMod(newMesh);

    // Insert mesh refinement into polyTopoChange.
    collapser.setRefinement(allPointInfo, newMeshMod);

    Info<< indent << "Apply changes to the current mesh" << endl;

    // Apply changes to current mesh
    autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
    (
        newMesh,
        false
    );
    const mapPolyMesh& newMap = newMapPtr();

    // Update fields
    newMesh.updateMesh(newMap);
    if (newMap.hasMotionPoints())
    {
        newMesh.movePoints(newMap.preMotionPoints());
    }
    updateSets(newMap);

    // Synchronise the factors
    mapOldMeshEdgeFieldToNewMesh
    (
        newMesh,
        newMap.pointMap(),
        newMeshMinEdgeLen
    );

    updateOldToNewPointMap
    (
        newMap.reversePointMap(),
        origToCurrentPointMap
    );

    updatePointPriorities(newMesh, newMap.pointMap());

    return nLocalCollapse;
}


void Foam::polyMeshFilter::updatePointErrorCount
(
    const bitSet& isErrorPoint,
    const labelList& oldToNewMesh,
    labelList& pointErrorCount
) const
{
    forAll(mesh_.points(), pI)
    {
        if (isErrorPoint[oldToNewMesh[pI]])
        {
            pointErrorCount[pI]++;
        }
    }
}


void Foam::polyMeshFilter::checkMeshEdgesAndRelaxEdges
(
    const polyMesh& newMesh,
    const labelList& oldToNewMesh,
    const bitSet& isErrorPoint,
    const labelList& pointErrorCount
)
{
    const edgeList& edges = mesh_.edges();

    forAll(edges, edgeI)
    {
        const edge& e = edges[edgeI];
        label newStart = oldToNewMesh[e[0]];
        label newEnd = oldToNewMesh[e[1]];

        if
        (
            pointErrorCount[e[0]] >= maxPointErrorCount()
         || pointErrorCount[e[1]] >= maxPointErrorCount()
        )
        {
            minEdgeLen_[edgeI] = -1;
        }

        if
        (
            (newStart >= 0 && isErrorPoint[newStart])
         || (newEnd >= 0 && isErrorPoint[newEnd])
        )
        {
            minEdgeLen_[edgeI] *= edgeReductionFactor();
        }
    }

    syncTools::syncEdgeList(mesh_, minEdgeLen_, minEqOp<scalar>(), scalar(0));

    for (label smoothIter = 0; smoothIter < maxSmoothIters(); ++smoothIter)
    {
        // Smooth minEdgeLen
        forAll(mesh_.edges(), edgeI)
        {
            const edge& e = mesh_.edges()[edgeI];

            scalar sumMinEdgeLen = 0;
            label nEdges = 0;

            forAll(e, pointi)
            {
                const labelList& pEdges = mesh_.pointEdges()[e[pointi]];

                forAll(pEdges, pEdgeI)
                {
                    const label pEdge = pEdges[pEdgeI];
                    sumMinEdgeLen += minEdgeLen_[pEdge];
                    nEdges++;
                }
            }

            minEdgeLen_[edgeI] = min
            (
                minEdgeLen_[edgeI],
                sumMinEdgeLen/nEdges
            );
        }

        syncTools::syncEdgeList
        (
            mesh_,
            minEdgeLen_,
            minEqOp<scalar>(),
            scalar(0)
        );
    }
}


void Foam::polyMeshFilter::checkMeshFacesAndRelaxEdges
(
    const polyMesh& newMesh,
    const labelList& oldToNewMesh,
    const bitSet& isErrorPoint,
    const labelList& pointErrorCount
)
{
    const faceList& faces = mesh_.faces();

    forAll(faces, facei)
    {
        const face& f = faces[facei];

        forAll(f, fpI)
        {
            const label ptIndex = oldToNewMesh[f[fpI]];

            if (pointErrorCount[f[fpI]] >= maxPointErrorCount())
            {
                faceFilterFactor_[facei] = -1;
            }

            if (isErrorPoint[ptIndex])
            {
                faceFilterFactor_[facei] *= faceReductionFactor();

                break;
            }
        }
    }

    syncTools::syncFaceList(mesh_, faceFilterFactor_, minEqOp<scalar>());

    for (label smoothIter = 0; smoothIter < maxSmoothIters(); ++smoothIter)
    {
        // Smooth faceFilterFactor
        forAll(faces, facei)
        {
            const labelList& fEdges = mesh_.faceEdges()[facei];

            scalar sumFaceFilterFactors = 0;
            label nFaces = 0;

            // This is important: Only smooth around faces that share an
            // edge with a bad face
            bool skipFace = true;

            forAll(fEdges, fEdgeI)
            {
                const labelList& eFaces = mesh_.edgeFaces()[fEdges[fEdgeI]];

                forAll(eFaces, eFacei)
                {
                    const label eFace = eFaces[eFacei];

                    const face& f = faces[eFace];

                    forAll(f, fpI)
                    {
                        const label ptIndex = oldToNewMesh[f[fpI]];

                        if (isErrorPoint[ptIndex])
                        {
                            skipFace = false;
                            break;
                        }
                    }

                    if (eFace != facei)
                    {
                        sumFaceFilterFactors += faceFilterFactor_[eFace];
                        nFaces++;
                    }
                }
            }

            if (skipFace)
            {
                continue;
            }

            faceFilterFactor_[facei] = min
            (
                faceFilterFactor_[facei],
                sumFaceFilterFactors/nFaces
            );
        }

        // Face filter factor needs to be synchronised!
        syncTools::syncFaceList(mesh_, faceFilterFactor_, minEqOp<scalar>());
    }
}


void Foam::polyMeshFilter::updatePointPriorities
(
    const polyMesh& newMesh,
    const labelList& pointMap
)
{
    labelList newPointPriority(newMesh.nPoints(), labelMin);
    const labelList& currPointPriority = pointPriority_();

    forAll(newPointPriority, ptI)
    {
        const label newPointToOldPoint = pointMap[ptI];
        const label origPointPriority = currPointPriority[newPointToOldPoint];

        newPointPriority[ptI] = max(origPointPriority, newPointPriority[ptI]);
    }

    syncTools::syncPointList
    (
        newMesh,
        newPointPriority,
        maxEqOp<label>(),
        labelMin
    );

    pointPriority_.reset(new labelList(newPointPriority));
}


void Foam::polyMeshFilter::printScalarFieldStats
(
    const string desc,
    const scalarField& fld
) const
{
    scalar sum = 0;
    scalar validElements = 0;
    scalar min = GREAT;
    scalar max = -GREAT;

    forAll(fld, i)
    {
        const scalar fldElement = fld[i];

        if (fldElement >= 0)
        {
            sum += fldElement;

            if (fldElement < min)
            {
                min = fldElement;
            }

            if (fldElement > max)
            {
                max = fldElement;
            }

            validElements++;
        }
    }

    reduce(sum, sumOp<scalar>());
    reduce(min, minOp<scalar>());
    reduce(max, maxOp<scalar>());
    reduce(validElements, sumOp<label>());
    const label totFieldSize = returnReduce(fld.size(), sumOp<label>());

    Info<< incrIndent << indent << desc
        << ": min = " << min
        << " av = " << sum/(validElements + SMALL)
        << " max = " << max << nl
        << indent
        << "    " << validElements << " / " << totFieldSize << " elements used"
        << decrIndent << endl;
}


void Foam::polyMeshFilter::mapOldMeshEdgeFieldToNewMesh
(
    const polyMesh& newMesh,
    const labelList& pointMap,
    scalarField& newMeshMinEdgeLen
) const
{
    scalarField tmp(newMesh.nEdges());

    const edgeList& newEdges = newMesh.edges();

    forAll(newEdges, newEdgeI)
    {
        const edge& newEdge = newEdges[newEdgeI];
        const label pStart = newEdge.start();
        const label pEnd = newEdge.end();

        tmp[newEdgeI] = min
        (
            newMeshMinEdgeLen[pointMap[pStart]],
            newMeshMinEdgeLen[pointMap[pEnd]]
        );
    }

    newMeshMinEdgeLen.transfer(tmp);

    syncTools::syncEdgeList
    (
        newMesh,
        newMeshMinEdgeLen,
        maxEqOp<scalar>(),
        scalar(0)
    );
}


void Foam::polyMeshFilter::mapOldMeshFaceFieldToNewMesh
(
    const polyMesh& newMesh,
    const labelList& faceMap,
    scalarField& newMeshFaceFilterFactor
) const
{
    scalarField tmp(newMesh.nFaces());

    forAll(faceMap, newFacei)
    {
        const label oldFacei = faceMap[newFacei];

        tmp[newFacei] = newMeshFaceFilterFactor[oldFacei];
    }

    newMeshFaceFilterFactor.transfer(tmp);

    syncTools::syncFaceList
    (
        newMesh,
        newMeshFaceFilterFactor,
        maxEqOp<scalar>()
    );
}


void Foam::polyMeshFilter::updateOldToNewPointMap
(
    const labelList& currToNew,
    labelList& origToCurrentPointMap
) const
{
    forAll(origToCurrentPointMap, origPointi)
    {
        label oldPointi = origToCurrentPointMap[origPointi];

        if (oldPointi != -1)
        {
            label newPointi = currToNew[oldPointi];

            if (newPointi >= 0)
            {
                origToCurrentPointMap[origPointi] = newPointi;
            }
            else if (newPointi == -1)
            {
                origToCurrentPointMap[origPointi] = -1;
            }
            else
            {
                origToCurrentPointMap[origPointi] = -newPointi-2;
            }
        }
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::polyMeshFilter::polyMeshFilter(const fvMesh& mesh)
:
    polyMeshFilterSettings
    (
        IOdictionary
        (
            IOobject
            (
                "collapseDict",
                mesh.time().system(),
                mesh.time(),
                IOobject::MUST_READ,
                IOobject::NO_WRITE
            )
        )
    ),
    mesh_(mesh),
    newMeshPtr_(),
    originalPointPriority_(mesh.nPoints(), labelMin),
    pointPriority_(),
    minEdgeLen_(),
    faceFilterFactor_()
{
    writeSettings(Info);
}


Foam::polyMeshFilter::polyMeshFilter
(
    const fvMesh& mesh,
    const labelList& pointPriority
)
:
    polyMeshFilterSettings
    (
        IOdictionary
        (
            IOobject
            (
                "collapseDict",
                mesh.time().system(),
                mesh.time(),
                IOobject::MUST_READ,
                IOobject::NO_WRITE
            )
        )
    ),
    mesh_(mesh),
    newMeshPtr_(),
    originalPointPriority_(pointPriority),
    pointPriority_(),
    minEdgeLen_(),
    faceFilterFactor_()
{
    writeSettings(Info);
}

Foam::polyMeshFilter::polyMeshFilter
(
    const fvMesh& mesh,
    const labelList& pointPriority,
    const dictionary& dict
)
:
    polyMeshFilterSettings(dict),
    mesh_(mesh),
    newMeshPtr_(),
    originalPointPriority_(pointPriority),
    pointPriority_(),
    minEdgeLen_(),
    faceFilterFactor_()
{
    writeSettings(Info);
}


// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //

Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
{
    minEdgeLen_.resize(mesh_.nEdges(), minLen());
    faceFilterFactor_.resize(mesh_.nFaces(), initialFaceLengthFactor());

    return filterFacesLoop(nOriginalBadFaces);
}


Foam::label Foam::polyMeshFilter::filter(const faceSet& fSet)
{
    minEdgeLen_.resize(mesh_.nEdges(), minLen());
    faceFilterFactor_.resize(mesh_.nFaces(), initialFaceLengthFactor());

    forAll(faceFilterFactor_, fI)
    {
        if (!fSet.found(fI))
        {
            faceFilterFactor_[fI] = -1;
        }
    }

    return filterFacesLoop(0);
}


Foam::label Foam::polyMeshFilter::filterEdges
(
    const label nOriginalBadFaces
)
{
    label nBadFaces = labelMax/2;
    label nPreviousBadFaces = labelMax;
    label nOuterIterations = 0;

    minEdgeLen_.resize(mesh_.nEdges(), minLen());
    faceFilterFactor_.clear();

    labelList pointErrorCount(mesh_.nPoints(), Zero);

    // Main loop
    // ~~~~~~~~~
    // It tries and do some collapses, checks the resulting mesh and
    // 'freezes' some edges (by marking in minEdgeLen) and tries again.
    // This will iterate ultimately to the situation where every edge is
    // frozen and nothing gets collapsed.
    while
    (
        nOuterIterations < maxIterations()
     && nBadFaces > nOriginalBadFaces
     && nBadFaces < nPreviousBadFaces
    )
    {
        Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
            << endl;

        printScalarFieldStats("Edge Filter Factor", minEdgeLen_);

        nPreviousBadFaces = nBadFaces;

        // Reset the new mesh to the old mesh
        newMeshPtr_ = copyMesh(mesh_);
        fvMesh& newMesh = newMeshPtr_();

        scalarField newMeshMinEdgeLen = minEdgeLen_;
        pointPriority_.reset(new labelList(originalPointPriority_));

        labelList origToCurrentPointMap(identity(newMesh.nPoints()));

        Info<< incrIndent;

        label nInnerIterations = 0;
        label nPrevLocalCollapse = labelMax;

        while (true)
        {
            Info<< nl
                << indent << "Inner iteration = " << nInnerIterations++ << nl
                << incrIndent << endl;

            label nLocalCollapse = filterEdges
            (
                newMesh,
                newMeshMinEdgeLen,
                origToCurrentPointMap
            );

            Info<< decrIndent;

            if
            (
                nLocalCollapse >= nPrevLocalCollapse
             || nLocalCollapse == 0
            )
            {
                Info<< decrIndent;
                break;
            }
            else
            {
                nPrevLocalCollapse = nLocalCollapse;
            }
        }

        // Mesh check
        // ~~~~~~~~~~~~~~~~~~
        // Do not allow collapses in regions of error.
        // Updates minEdgeLen, nRelaxedEdges

        if (controlMeshQuality())
        {
            bitSet isErrorPoint(newMesh.nPoints());
            nBadFaces = edgeCollapser::checkMeshQuality
            (
                newMesh,
                meshQualityCoeffDict(),
                isErrorPoint
            );

            Info<< nl << "    Number of bad faces     : " << nBadFaces << nl
                << "    Number of marked points : "
                << returnReduce(isErrorPoint.count(), sumOp<unsigned int>())
                << endl;

            updatePointErrorCount
            (
                isErrorPoint,
                origToCurrentPointMap,
                pointErrorCount
            );

            checkMeshEdgesAndRelaxEdges
            (
                newMesh,
                origToCurrentPointMap,
                isErrorPoint,
                pointErrorCount
            );
        }
        else
        {
            return -1;
        }
    }

    return nBadFaces;
}


const Foam::autoPtr<Foam::fvMesh>& Foam::polyMeshFilter::filteredMesh() const
{
    return newMeshPtr_;
}


const Foam::autoPtr<Foam::labelList>&
Foam::polyMeshFilter::pointPriority() const
{
    return pointPriority_;
}


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