/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | cfMesh: A library for mesh generation
   \\    /   O peration     |
    \\  /    A nd           | www.cfmesh.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2014-2017 Creative Fields, Ltd.
-------------------------------------------------------------------------------
Author
     Franjo Juretic (franjo.juretic@c-fields.com)

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 "error.H"
#include "polyMeshGenModifier.H"
#include "edgeExtractor.H"
#include "meshSurfaceEngine.H"
#include "meshSurfaceEngineModifier.H"
#include "meshSurfaceOptimizer.H"
#include "meshOctree.H"
#include "triSurf.H"
#include "triSurfModifier.H"
#include "helperFunctions.H"
#include "DynList.H"
#include "labelPair.H"
#include "labelledScalar.H"
#include "labelledPoint.H"
#include "refLabelledPoint.H"
#include "HashSet.H"
#include "triSurfacePartitioner.H"
#include "triSurfaceClassifyEdges.H"
#include "demandDrivenData.H"

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

//#define DEBUGEdgeExtractor

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

void Foam::Module::edgeExtractor::faceEvaluator::calculateNeiPatchesParallel()
{
    otherFacePatch_.clear();

    const labelList& fPatches = extractor_.facePatch_;

    if (Pstream::parRun())
    {
        const meshSurfaceEngine& mse = extractor_.surfaceEngine();
        const VRWGraph& edgeFaces = mse.edgeFaces();
        const Map<label>& otherProc = mse.otherEdgeFaceAtProc();
        const Map<label>& globalToLocal =
            mse.globalToLocalBndEdgeAddressing();

        // create communication matrix
        std::map<label, labelLongList> exchangeData;
        const DynList<label>& neiProcs = mse.beNeiProcs();
        forAll(neiProcs, procI)
            exchangeData.insert
            (
                std::make_pair(neiProcs[procI], labelLongList())
            );

        forAllConstIters(globalToLocal, it)
        {
            const label beI = it();

            if (edgeFaces.sizeOfRow(beI) == 1)
            {
                labelLongList& dts = exchangeData[otherProc[beI]];
                // send data as follows:
                // 1. global edge label
                // 2. patch of the attached boundary face
                dts.append(it.key());
                dts.append(fPatches[edgeFaces(beI, 0)]);
            }
        }

        labelLongList receivedData;
        help::exchangeMap(exchangeData, receivedData);

        label counter(0);
        while (counter < receivedData.size())
        {
            const label beI = globalToLocal[receivedData[counter++]];
            const label fPatch = receivedData[counter++];

            otherFacePatch_.insert(beI, fPatch);
        }
    }
}


void Foam::Module::edgeExtractor::faceEvaluator::
calculateNeiPatchesParallelNewPatches()
{
    if (newOtherFacePatchPtr_)
        return;

    if (!newBoundaryPatchesPtr_)
        FatalErrorInFunction
            << "newBoundaryPatchesPtr_ are not set" << exit(FatalError);

    newOtherFacePatchPtr_ = new Map<label>();
    Map<label>& otherFacePatch = *newOtherFacePatchPtr_;

    const labelList& fPatches = *newBoundaryPatchesPtr_;

    if (Pstream::parRun())
    {
        const meshSurfaceEngine& mse = extractor_.surfaceEngine();
        const VRWGraph& edgeFaces = mse.edgeFaces();
        const Map<label>& otherProc = mse.otherEdgeFaceAtProc();
        const Map<label>& globalToLocal =
            mse.globalToLocalBndEdgeAddressing();

        // create communication matrix
        std::map<label, labelLongList> exchangeData;
        const DynList<label>& neiProcs = mse.beNeiProcs();
        forAll(neiProcs, procI)
            exchangeData.insert
            (
                std::make_pair(neiProcs[procI], labelLongList())
            );

        forAllConstIters(globalToLocal, it)
        {
            const label beI = it();

            if (edgeFaces.sizeOfRow(beI) == 1)
            {
                labelLongList& dts = exchangeData[otherProc[beI]];
                // send data as follows:
                // 1. global edge label
                // 2. patch of the attached boundary face
                dts.append(it.key());
                dts.append(fPatches[edgeFaces(beI, 0)]);
            }
        }

        labelLongList receivedData;
        help::exchangeMap(exchangeData, receivedData);

        label counter(0);
        while (counter < receivedData.size())
        {
            const label beI = globalToLocal[receivedData[counter++]];
            const label fPatch = receivedData[counter++];

            otherFacePatch.insert(beI, fPatch);
        }
    }
}


void Foam::Module::edgeExtractor::faceEvaluator::neiFacesOverEdges
(
    const label bfI,
    DynList<label>& neiFaces
) const
{
    const meshSurfaceEngine& mse = extractor_.surfaceEngine();

    const VRWGraph& faceEdges = mse.faceEdges();
    const VRWGraph& edgeFaces = mse.edgeFaces();

    neiFaces.setSize(faceEdges.sizeOfRow(bfI));

    forAllRow(faceEdges, bfI, feI)
    {
        const label beI = faceEdges(bfI, feI);

        if (edgeFaces.sizeOfRow(beI) == 2)
        {
            neiFaces[feI] = edgeFaces(beI, 0);
            if (neiFaces[feI] == bfI)
                neiFaces[feI] = edgeFaces(beI, 1);
        }
        else
        {
            neiFaces[feI] = -1;
        }
    }
}


void Foam::Module::edgeExtractor::faceEvaluator::neiFacesProcs
(
    const label bfI,
    DynList<label>& neiProcs
) const
{
    const meshSurfaceEngine& mse = extractor_.surfaceEngine();

    const VRWGraph& faceEdges = mse.faceEdges();

    neiProcs.setSize(faceEdges.sizeOfRow(bfI));
    neiProcs = Pstream::myProcNo();

    if (Pstream::parRun())
    {
        const Map<label>& otherFaceProc = mse.otherEdgeFaceAtProc();

        forAllRow(faceEdges, bfI, feI)
        {
            const label beI = faceEdges(bfI, feI);

            const Map<label>::const_iterator it = otherFaceProc.find(beI);

            if (it != otherFaceProc.end())
                neiProcs[feI] = it();
        }
    }
}


void Foam::Module::edgeExtractor::faceEvaluator::neiPatchesOverEdges
(
    const label bfI,
    const labelList& fPatches,
    const Map<label>& otherFacePatch,
    DynList<label> &neiPatches
) const
{
    const meshSurfaceEngine& mse = extractor_.surfaceEngine();

    const VRWGraph& faceEdges = mse.faceEdges();
    const VRWGraph& edgeFaces = mse.edgeFaces();

    neiPatches.setSize(faceEdges.sizeOfRow(bfI));

    forAllRow(faceEdges, bfI, feI)
    {
        const label beI = faceEdges(bfI, feI);

        if (edgeFaces.sizeOfRow(beI) == 2)
        {
            label nei = edgeFaces(beI, 0);
            if (nei == bfI)
                nei = edgeFaces(beI, 1);

            neiPatches[feI] = fPatches[nei];
        }
        else if (Pstream::parRun() && (edgeFaces.sizeOfRow(beI) == 1))
        {
            neiPatches[feI] = otherFacePatch[beI];
        }
    }
}


Foam::label Foam::Module::edgeExtractor::faceEvaluator::bestPatchTopological
(
    const DynList<label>& neiPatches,
    const label currentPatch
)
{
    // find indices of all neighbour patches
    DynList<label> allNeiPatches;
    forAll(neiPatches, i)
        allNeiPatches.appendUniq(neiPatches[i]);

    if ((allNeiPatches.size() == 1) && (allNeiPatches[0] == currentPatch))
        return currentPatch;

    // counter the number of neighbours in a patch
    Map<label> nNeiInPatch(allNeiPatches.size());
    forAll(allNeiPatches, i)
        nNeiInPatch.insert(allNeiPatches[i], 0);
    forAll(neiPatches, eI)
        ++nNeiInPatch[neiPatches[eI]];

    label newPatch = -1;
    label nNeiEdges(0);
    forAllConstIters(nNeiInPatch, it)
    {
        if (it() > nNeiEdges)
        {
            newPatch = it.key();
            nNeiEdges = it();
        }
        else if
        (
            (it() == nNeiEdges) && (it.key() == currentPatch)
        )
        {
            newPatch = it.key();
        }
    }

    // do not swap if the situation allows for more than one edge
    // shared with faces in other patches than the dominant one
    if (nNeiEdges < (neiPatches.size() - 1))
        return currentPatch;

    // do not swap if the best face is in the current patch
    if ((newPatch < 0) || (newPatch == currentPatch))
        return currentPatch;

    // check whether the edges shared ith the neighbour patch form
    // a singly linked chain
    DynList<bool> sharedEdge;
    sharedEdge.setSize(neiPatches.size());
    sharedEdge = false;

    forAll(neiPatches, eI)
        if (neiPatches[eI] == newPatch)
            sharedEdge[eI] = true;

    if (help::areElementsInChain(sharedEdge))
        return newPatch;

    return currentPatch;
}


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

Foam::Module::edgeExtractor::faceEvaluator::faceEvaluator
(
    const edgeExtractor& ee
)
:
    extractor_(ee),
    otherFacePatch_(),
    newBoundaryPatchesPtr_(nullptr),
    newOtherFacePatchPtr_(nullptr)
{
    if (Pstream::parRun())
        calculateNeiPatchesParallel();
}


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

Foam::Module::edgeExtractor::faceEvaluator::~faceEvaluator()
{
    deleteDemandDrivenData(newOtherFacePatchPtr_);
}


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

void Foam::Module::edgeExtractor::faceEvaluator::setNewBoundaryPatches
(
    const labelList& newBoudaryPatches
)
{
    newBoundaryPatchesPtr_ = &newBoudaryPatches;

    if (Pstream::parRun())
        calculateNeiPatchesParallelNewPatches();
}


void Foam::Module::edgeExtractor::faceEvaluator::neiPatchesOverEdges
(
    const label bfI,
    DynList<label>& neiPatches
) const
{
    neiPatchesOverEdges
    (
        bfI,
        extractor_.facePatch_,
        otherFacePatch_,
        neiPatches
    );
}


Foam::label Foam::Module::edgeExtractor::faceEvaluator::bestPatchTopological
(
    const label bfI
) const
{
    // get neighbour patches over edges
    DynList<label> neiPatches;
    neiPatchesOverEdges
    (
        bfI,
        extractor_.facePatch_,
        otherFacePatch_,
        neiPatches
    );

    return bestPatchTopological(neiPatches, extractor_.facePatch_[bfI]);
}


Foam::label
Foam::Module::edgeExtractor::faceEvaluator::bestPatchAfterModification
(
    const label bfI
) const
{
    const label patchI = newBoundaryPatchesPtr_->operator[](bfI);

    if (patchI != extractor_.facePatch_[bfI])
    {
        DynList<label> newNeiPatches, oldNeiPatches;
        neiPatchesOverEdges
        (
            bfI,
            *newBoundaryPatchesPtr_,
            *newOtherFacePatchPtr_,
            newNeiPatches
        );

        neiPatchesOverEdges
        (
            bfI,
            extractor_.facePatch_,
            otherFacePatch_,
            oldNeiPatches
        );

        DynList<label> neiFaces, neiProcs;
        neiFacesOverEdges(bfI, neiFaces);
        neiFacesProcs(bfI, neiProcs);

        forAll(neiFaces, eI)
        {
            const label origPatchI = oldNeiPatches[eI];
            const label newPatchI = newNeiPatches[eI];

            if (neiFaces[eI] > bfI)
            {
                if (newPatchI != origPatchI)
                    newNeiPatches[eI] = origPatchI;
            }
            else if (neiFaces[eI] < 0)
            {
                if (neiProcs[eI] > Pstream::myProcNo())
                    newNeiPatches[eI] = origPatchI;
            }
        }

        return bestPatchTopological(newNeiPatches, patchI);
    }

    return patchI;
}


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

void Foam::Module::edgeExtractor::cornerEvaluator::createParallelAddressing()
{
    const labelHashSet& corners = partitioner_.corners();

    const labelList& facePatch = extractor_.facePatch_;

    typedef Map<label> mapType;

    const meshSurfaceEngine& mse = extractor_.surfaceEngine();
    const faceList::subList& bFaces = mse.boundaryFaces();
    const pointFieldPMG& points = mse.points();
    const labelList& bp = mse.bp();
    const VRWGraph& pointFaces = mse.pointFaces();
    const labelList& globalLabel = mse.globalBoundaryPointLabel();
    const mapType& globalToLocal = mse.globalToLocalBndPointAddressing();
    const VRWGraph& bpAtProcs = mse.bpAtProcs();
    const DynList<label>& neiProcs = mse.bpNeiProcs();

    typedef std::map<label, LongList<labelledPoint>> exchangeMapType;
    exchangeMapType exchangeData;
    forAll(neiProcs, i)
        exchangeData[neiProcs[i]].clear();

    faceMap_.clear();
    faceAtProc_.clear();
    facePatches_.clear();

    forAllConstIters(globalToLocal, it)
    {
        const label bpI = it();

        if (corners.found(bpI))
        {
            forAllRow(bpAtProcs, bpI, i)
            {
                const label neiProc = bpAtProcs(bpI, i);

                if (neiProc == Pstream::myProcNo())
                    continue;

                LongList<labelledPoint>& dts = exchangeData[neiProc];

                // data is send in the ollowing form
                // 1. globsl label of the current point
                // 2. number of faces attached to the point
                // 3. for each face send its patch, number of points and point
                dts.append(labelledPoint(it.key(), vector::zero));
                dts.append
                (
                    labelledPoint(pointFaces.sizeOfRow(bpI), vector::zero)
                );
                forAllRow(pointFaces, bpI, i)
                {
                    const face& bf = bFaces[pointFaces(bpI, i)];
                    dts.append
                    (
                        labelledPoint
                        (
                            facePatch[pointFaces(bpI, i)],
                            vector::zero
                        )
                    );
                    dts.append(labelledPoint(bf.size(), vector::zero));

                    DynList<labelledPoint> lpf;
                    forAll(bf, pI)
                    {
                        const labelledPoint lp
                        (
                            globalLabel[bp[bf[bpI]]],
                            points[bf[pI]]
                        );

                        dts.append(lp);
                        lpf.append(lp);
                    }

                    faceMap_[bpI].append(lpf);
                    faceAtProc_[bpI].append(Pstream::myProcNo());
                    facePatches_[bpI].append(facePatch[pointFaces(bpI, i)]);
                }
            }
        }
    }

    std::map<label, List<labelledPoint>> receivedDataMap;
    help::exchangeMap(exchangeData, receivedDataMap);

    forAllConstIters(receivedDataMap, it)
    {
        const label procI = it->first;
        const List<labelledPoint>& receivedData = it->second;

        for (label i = 0; i < receivedData.size(); /*nil*/ )
        {
            const label bpI = globalToLocal[receivedData[i++].pointLabel()];

            const label nFaces = receivedData[i++].pointLabel();

            for (label fI = 0; fI < nFaces; ++fI)
            {
                const label patchI = receivedData[i++].pointLabel();
                const label size = receivedData[i++].pointLabel();

                DynList<labelledPoint> lpf(size);
                forAll(lpf, pI)
                {
                    lpf[pI] = receivedData[i++];
                }

                faceMap_[bpI].append(lpf);
                faceAtProc_[bpI].append(procI);
                facePatches_[bpI].append(patchI);
            }
        }
    }

    // sort the faces in the counter-clockwise order
    forAllIters(faceAtProc_, it)
    {
        const label bpI = it->first;

        DynList<label>& faceProc = it->second;
        DynList<label>& fPatches = facePatches_[bpI];
        DynList<DynList<labelledPoint, 6>>& pFaces = faceMap_[bpI];

        forAll(pFaces, i)
        {
            const DynList<labelledPoint, 6>& lpf = pFaces[i];

            label pos(-1);
            forAll(lpf, pI)
            {
                if (lpf[i].pointLabel() == globalLabel[bpI])
                {
                    pos = pI;
                    break;
                }
            }

            for (label j = i + 2; j < pFaces.size(); ++j)
            {
                const DynList<labelledPoint, 6>& olpf = pFaces[j];

                if (olpf.found(lpf[lpf.rcIndex(pos)]))
                {
                    label helper = faceProc[j];
                    faceProc[j] = faceProc[i + 1];
                    faceProc[i + 1] = helper;

                    helper = fPatches[j];
                    fPatches[j] = fPatches[i + 1];
                    fPatches[i + 1] = helper;

                    DynList<labelledPoint, 6> lpfHelper;
                    lpfHelper = pFaces[j];
                    pFaces[j] = pFaces[i + 1];
                    pFaces[i + 1] = lpfHelper;
                }
            }
        }
    }
}


void Foam::Module::edgeExtractor::cornerEvaluator::sortedFacesAtPoint
(
    const label bpI,
    DynList<label>& pFaces
) const
{
    const meshSurfaceEngine& mse = extractor_.surfaceEngine();
    const faceList::subList& bFaces = mse.boundaryFaces();
    const VRWGraph& pointFaces = mse.pointFaces();
    const VRWGraph& pointInFace = mse.pointInFaces();

    pFaces = pointFaces[bpI];

    forAll(pFaces, i)
    {
        const face& bf = bFaces[pFaces[i]];
        const label pos = pointFaces.find(bpI, pFaces[i]);

        const edge e = bf.faceEdge(bf.rcIndex(pointInFace(bpI, pos)));

        for (label j = i + 2; j < pFaces.size(); ++j)
        {
            if (bFaces[pFaces[j]].which(e.start()) >= 0)
            {
                label bfJ = pFaces[i + 1];
                pFaces[i + 1] = pFaces[j];
                pFaces[j] = bfJ;
            }
        }
    }
}


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

Foam::Module::edgeExtractor::cornerEvaluator::cornerEvaluator
(
    const edgeExtractor& ee,
    const meshSurfacePartitioner& mPart
)
:
    extractor_(ee),
    partitioner_(mPart),
    faceMap_(),
    facePatches_(),
    faceAtProc_()
{
    if (Pstream::parRun())
        createParallelAddressing();
}


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

void Foam::Module::edgeExtractor::cornerEvaluator::improveCorners
(
    labelList& newBoundaryPatches
)
{
    const meshOctree& mo = extractor_.meshOctree_;

    const labelHashSet& corners = partitioner_.corners();

    const meshSurfaceEngine& mse = extractor_.surfaceEngine();
    const labelList& bPoints = mse.boundaryPoints();
    const pointFieldPMG& points = mse.points();
    const faceList::subList& bFaces = mse.boundaryFaces();

    for (const label corneri : corners)
    {
        const point& p = points[bPoints[corneri]];

        if (faceMap_.find(corneri) == faceMap_.end())
        {
            const labelList& facePatch = extractor_.facePatch_;
            Info<< "Checking corner " << corneri << endl;

            // get faces attached to this point
            // sorted in the counter-clocwise order
            DynList<label> pFaces;
            sortedFacesAtPoint(corneri, pFaces);

            Info<< "Faces at corner " << pFaces << endl;

            // find feature edges attached to the point
            DynList<labelPair> edgePatches;
            DynList<label> edgeIndex;
            forAll(pFaces, i)
            {
                const label patch0 = facePatch[pFaces[i]];
                const label patch1 = facePatch[pFaces[pFaces.fcIndex(i)]];

                if (patch0 != patch1)
                {
                    edgeIndex.append(i);
                    edgePatches.append(labelPair(patch0, patch1));
                }
            }

            Info<< "Edge patches " << edgePatches << endl;
            Info<< "Edge indices " << edgeIndex << endl;

            // find the best aligned feature edge to the feature edge
            DynList<label> bestEdgeIndices(edgePatches.size());
            forAll(edgePatches, i)
            {
                const label patch0 = edgePatches[i].first();
                const label patch1 = edgePatches[i].second();
                DynList<label> patches(2);
                patches[0] = patch0;
                patches[1] = patch1;

                // find the proection of the first point onto the feature edge
                point mps;
                scalar dSqS;
                mo.findNearestPointToPatches(mps, dSqS, p, patches);

                scalar bestAlignment(0.0);
                label bestEdgeIndex(-1);

                forAll(pFaces, pfI)
                {
                    const face& bf = bFaces[pFaces[pfI]];

                    const label pos = bf.which(bPoints[corneri]);

                    const point& pe = points[bf[bf.rcIndex(pos)]];

                    // calculate the vector of current edge
                    vector ev = pe - p;
                    ev /= (mag(ev) + VSMALL);

                    // project the endpoint onto the feature edge
                    point mpe;
                    scalar dSqE;
                    mo.findNearestPointToPatches(mpe, dSqE, pe, patches);

                    vector fv = mpe - mps;
                    fv /= (mag(fv) + VSMALL);

                    // calculate alignment metrix between the current edge
                    // and the feature edge
                    const scalar align = 0.5*(1.0 + (ev & fv));

                    if (align > bestAlignment)
                    {
                        bestAlignment = align;
                        bestEdgeIndex = pfI;
                    }
                }

                bestEdgeIndices[i] = bestEdgeIndex;
            }

            Info<< "\nPatches of edges at corner " << corneri
                << " are " << edgePatches << endl;
            Info<< "Orig edge indices " << edgeIndex << endl;
            Info<< "New edge indices " << bestEdgeIndices << endl;

            DynList<label> newFacePatches(pFaces.size());
            forAll(bestEdgeIndices, i)
            {
                const label patch = edgePatches[i].first();
                label index = bestEdgeIndices[i];

                while (index != bestEdgeIndices[bestEdgeIndices.fcIndex(i)])
                {
                    newFacePatches[index] = patch;
                    index = (index + 1) % pFaces.size();
                }
            }

            bool allPatchesRemain(true);
            forAll(edgePatches, i)
            {
                if (!newFacePatches.found(edgePatches[i].first()))
                {
                    allPatchesRemain = false;
                    break;
                }
            }

            Info<< "New patches at point" << newFacePatches << endl;
            if (allPatchesRemain)
            {
                forAll(pFaces, i)
                    newBoundaryPatches[pFaces[i]] = newFacePatches[i];
            }
        }
        else
        {
        }
    }
}


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

bool Foam::Module::edgeExtractor::findCornerCandidates()
{
    bool changed(false);

    const triSurf& surf = meshOctree_.surface();
    const pointField& sp = surf.points();

    // create the surface partitioner and find the corners on the surface mesh
    triSurfacePartitioner surfPartitioner(surf);

    const labelList& corners = surfPartitioner.corners();

    Map<label> surfacePointToCorner;
    forAll(corners, cornerI)
        surfacePointToCorner.insert(corners[cornerI], cornerI);

    List<labelledScalar> nearestPoint
    (
        corners.size(),
        labelledScalar(0, VGREAT)
    );

    // calculate the search range
    const meshSurfaceEngine& mse = this->surfaceEngine();
    const pointFieldPMG& points = mesh_.points();
    const labelList& bPoints = mse.boundaryPoints();
    const labelList& bp = mse.bp();
    const faceList::subList& bFaces = mse.boundaryFaces();

    scalarList searchRange(bPoints.size(), 0.0);
    forAll(bFaces, bfI)
    {
        const face& bf = bFaces[bfI];
        const point c = bf.centre(points);

        forAll(bf, pI)
        {
            const scalar d = 2.0*Foam::mag(c - points[bf[pI]]);
            const label bpI = bp[bf[pI]];

            searchRange[bpI] = Foam::max(d, searchRange[bpI]);
        }
    }

    if (Pstream::parRun())
    {
        const VRWGraph& bpAtProcs = mse.beAtProcs();
        const DynList<label>& bpNeiProcs = mse.bpNeiProcs();
        const Map<label>& globalToLocal =
            mse.globalToLocalBndPointAddressing();

        std::map<label, LongList<labelledScalar>> exchangeData;
        forAll(bpNeiProcs, i)
            exchangeData.insert
            (
                std::make_pair(bpNeiProcs[i], LongList<labelledScalar>())
            );

        forAllConstIters(globalToLocal, iter)
        {
            const label bpI = iter();

            forAllRow(bpAtProcs, bpI, i)
            {
                const label neiProc = bpAtProcs(bpI, i);

                if (neiProc == Pstream::myProcNo())
                    continue;

                exchangeData[neiProc].append
                (
                    labelledScalar(iter.key(), searchRange[bpI])
                );
            }
        }

        LongList<labelledScalar> receivedData;
        help::exchangeMap(exchangeData, receivedData);

        forAll(receivedData, i)
        {
            const label bpI = globalToLocal[receivedData[i].scalarLabel()];

            searchRange[bpI] =
                Foam::max(searchRange[bpI], receivedData[i].value());
        }
    }

    DynList<label> containedTriangles;
    # ifdef USE_OMP
    # pragma omp parallel for schedule(dynamic, 40) private(containedTriangles)
    # endif
    forAll(bPoints, bpI)
    {
        const point& p = points[bPoints[bpI]];
        const scalar s = searchRange[bpI];

        boundBox bb(p - vector(s, s, s), p + vector(s, s, s));
        meshOctree_.findTrianglesInBox(bb, containedTriangles);

        // find the nearest corner on the surface mesh
        forAll(containedTriangles, i)
        {
            const labelledTri& tri = surf[containedTriangles[i]];

            forAll(tri, pI)
            {
                const label spI = tri[pI];
                if (!surfacePointToCorner.found(spI))
                    continue;

                const label cornerI = surfacePointToCorner[spI];

                const scalar d = Foam::magSqr(sp[spI] - points[bPoints[bpI]]);

                if (nearestPoint[cornerI].value() > d)
                {
                    // update the nearest point to the found corner
                    # ifdef USE_OMP
                    # pragma omp critical
                    # endif
                    {
                        nearestPoint[cornerI] = labelledScalar(bpI, d);
                    }
                }
            }
        }
    }

    # ifdef DEBUGEdgeExtractor
    Info<< "Nearest points to corners " << nearestPoint << endl;
    # endif

    return changed;
}


bool Foam::Module::edgeExtractor::checkCorners()
{
    bool changed(false);

    # ifdef DEBUGEdgeExtractor
    const triSurf* surfPtr = surfaceWithPatches();
    surfPtr->writeSurface("checkCornersBefore.fms");
    deleteDemandDrivenData(surfPtr);
    # endif

    const pointFieldPMG& points = mesh_.points();
    const meshSurfaceEngine& mse = this->surfaceEngine();
    const labelList& bPoints = mse.boundaryPoints();
    const faceList::subList& bFaces = mse.boundaryFaces();
    const labelList& bp = mse.bp();
    const edgeList& edges = mse.edges();
    const VRWGraph& pointFaces = mse.pointFaces();
    const VRWGraph& faceEdges = mse.faceEdges();
    const VRWGraph& edgeFaces = mse.edgeFaces();

    const triSurfacePartitioner& sPartitioner = this->partitioner();
    const List<labelHashSet>& patchPatches = sPartitioner.patchPatches();

    // allocate a copy of boundary patches
    labelList newBoundaryPatches(facePatch_.size());

    label nCorrected;
    Map<label> otherProcNewPatch;

    label nIteration(0);

    do
    {
        nCorrected = 0;
        newBoundaryPatches = facePatch_;

        // check whether there exist situations where a boundary face
        // is surrounded by more faces in different patches than the
        // faces in the current patch
        if (Pstream::parRun())
        {
            findOtherFacePatchesParallel
            (
                otherProcNewPatch,
                &newBoundaryPatches
            );
        }

        // update the information which edges are assigned as feature edges
        meshSurfacePartitioner mPart(mse, newBoundaryPatches);

        // find corners in the current constelation
        const labelHashSet& corners = mPart.corners();
        const VRWGraph& pPatches = mPart.pointPatches();

        typedef std::map<label, label> mapType;
        typedef std::map<label, std::pair<point, scalar>> distMapType;

        mapType cornerIndex;
        distMapType nearestToCorner;

        # ifdef DEBUGEdgeExtractor
        Info<< "Finding corners " << endl;
        # endif

        // find nearest corners in the surface mesh
        for (const label bpI : corners)
        {
            const DynList<label> neiPatches = pPatches[bpI];

            const point& p = points[bPoints[bpI]];
            point pMap;
            scalar dSq;
            label nsp;

            if (!meshOctree_.findNearestCorner(pMap, dSq, nsp, p, neiPatches))
            {
                nsp = -1;
                meshOctree_.findNearestPointToPatches(pMap, dSq, p, neiPatches);
            }

            nearestToCorner[bpI] = std::make_pair(pMap, dSq);
            cornerIndex[bpI] = nsp;
        }

        std::map<label, DynList<label>> bpMappedAtSurfaceCorner;
        forAllConstIters(cornerIndex, mIt)
        {
            if (mIt->second < 0)
            {
                // the corner does not exist in the surface mesh
                // this may be a situation where parts of the surface are in
                // close proximity and are not topologically connected
                Warning << "Should not get in here. Not implemented" << endl;
            }
            else
            {
                bpMappedAtSurfaceCorner[mIt->second].append(mIt->first);
            }
        }

        # ifdef DEBUGEdgeExtractor
        forAllConstIters(bpMappedAtSurfaceCorner, mIt)
        {
            if (mIt->second.size() > 1)
            {
                Info<< "Surface corner " << mIt->first
                    << " mapped mesh points " << mIt->second << endl;
            }
        }
        # endif

        // for all edge nodes find their nearest counterparts in the surface
        # ifdef DEBUGEdgeExtractor
        Info<< "Finding nearest edges" << endl;
        # endif

        const labelHashSet& featureEdge = mPart.featureEdges();

        distMapType nearestToEdgePoint;
        mapType edgePointIndex;
        std::map<std::pair<label, label>, labelLongList> edgesSharedByPatches;

        for (const label beI : featureEdge)
        {
            // get the patches bounded by this feature edge
            DynList<label> patches(2);

            if (edgeFaces.sizeOfRow(beI) == 2)
            {
                patches[0] = newBoundaryPatches[edgeFaces(beI, 0)];
                patches[1] = newBoundaryPatches[edgeFaces(beI, 1)];
            }
            else if (edgeFaces.sizeOfRow(beI) == 1)
            {
                patches[0] = newBoundaryPatches[edgeFaces(beI, 0)];
                patches[1] = otherProcNewPatch[beI];
            }

            // store the edge into the right group
            std::pair<label, label> patchPair
            (
                Foam::min(patches[0], patches[1]),
                Foam::max(patches[0], patches[1])
            );
            edgesSharedByPatches[patchPair].append(beI);

            // check if some points have already been checked
            const edge& e = edges[beI];
            const label bps = bp[e.start()];
            const label bpe = bp[e.end()];

            DynList<label> checkPoints;
            if
            (
                corners.found(bps) &&
                (nearestToEdgePoint.find(bps) == nearestToEdgePoint.end())
            )
                checkPoints.append(bps);

            if
            (
                corners.found(bpe) &&
                (nearestToEdgePoint.find(bpe) == nearestToEdgePoint.end())
            )
                checkPoints.append(bpe);

            forAll(checkPoints, i)
            {
                const point& p = points[bPoints[checkPoints[i]]];
                point pMap;
                scalar dSq;
                label nse;

                if (!meshOctree_.findNearestEdgePoint(pMap, dSq, nse, p, patches))
                {
                    nse = -1;
                    meshOctree_.findNearestPointToPatches(pMap, dSq, p, patches);
                }

                nearestToEdgePoint[checkPoints[i]] = std::make_pair(pMap, dSq);
                edgePointIndex[checkPoints[i]] = nse;
            }
        }

        labelHashSet invalidFeatureEdges;
        forAllConstIters(edgesSharedByPatches, mIt)
        {
            const bool validConnection =
                patchPatches[mIt->first.first].found(mIt->first.second);

            # ifdef DEBUGEdgeExtractor
            Info<< "Patches " << mIt->first.first
                 << " and " << mIt->first.second
                 << " share " << mIt->second
                 << " edges " << validConnection << endl;
            # endif

            if (!validConnection)
            {
                // find surface facets in the vicinity of the edge
                // and check whether there exist various disconnected surface
                // parts in the vicinity of this group of edge

                const auto& invalidEdges = mIt->second;

                forAll(invalidEdges, i)
                    invalidFeatureEdges.insert(invalidEdges[i]);
            }
        }

        # ifdef DEBUGEdgeExtractor
        Info<< "Invalid feature edges " << invalidFeatureEdges << endl;
        ::exit(0);
        # endif

        // check the vicinity of the corner point and check whether
        // it shall be replaced by some other point in the vicinity
        std::map<label, DynList<labelPair>> facesAtCornerNewPatches;
        forAllConstIters(nearestToCorner, it)
        {
            const label bpI = it->first;

            // find all points connected to the corner point via a face
            labelHashSet nearPoints;
            std::map<label, DynList<label>> facesContainingPoint;
            forAllRow(pointFaces, bpI, pfI)
            {
                const label bfI = pointFaces(bpI, pfI);
                const face& bf = bFaces[bfI];

                forAll(bf, pI)
                {
                    nearPoints.insert(bf[pI]);
                    facesContainingPoint[bf[pI]].append(bfI);
                }
            }

            # ifdef DEBUGEdgeExtractor
            for (const label pointi : nearPoints)
            {
                Info<< "Near point " << pointi << endl;
                Info<< "Faces containing near point are "
                     << facesContainingPoint[pointi] << endl;
            }
            # endif

            // find the nearest point to the location where the corner
            // shall be mapped onto the surface mesh
            label bestPoint(-1);
            scalar minDistSq(VGREAT);
            for (const label pointi : nearPoints)
            {
                const point& p = points[pointi];

                const scalar dSq = magSqr(p - nearestToCorner[bpI].first);

                if (dSq < minDistSq)
                {
                    minDistSq = dSq;
                    bestPoint = pointi;
                }
            }

            if (bestPoint == bPoints[bpI])
            {
                # ifdef DEBUGEdgeExtractor
                Info<< "Current corner is nearest" << endl;
                # endif

                continue;
            }

            # ifdef DEBUGEdgeExtractor
            surfPtr = surfaceWithPatches(bpI);
            surfPtr->writeSurface
            (
                "corner_" + Foam::name(bpI)
              + "_iter_" + Foam::name(nIteration) +".fms"
            );
            deleteDemandDrivenData(surfPtr);
            Info<< "Best candidate " << bestPoint << endl;
            # endif

            // sort faces and edges at the corner in counter-clockwise order
            DynList<label> pFaces, pEdges;

            DynList<label> front;
            front.append(0);

            while (front.size())
            {
                const label fI = front.remove();
                const label bfI = pointFaces(bpI, fI);

                pFaces.append(bfI);

                const face& bf = bFaces[bfI];
                const label pos = bf.which(bPoints[bpI]);

                const label beI = faceEdges(bfI, bf.rcIndex(pos));

                pEdges.append(beI);

                label nei = edgeFaces(beI, 0);
                if (nei == bfI)
                    nei = edgeFaces(beI, 1);

                if (pFaces.found(nei) || (nei < 0))
                    continue;

                front.append(pointFaces.find(bpI, nei));
            }

            # ifdef DEBUGEdgeExtractor
            Info<< "Boundary point " << bpI
                 << " pFaces " << pFaces
                 << " pEdges " << pEdges << endl;
            forAll(pFaces, i)
            {
                Info<< "Face " << pFaces[i] << " has nodes " << bFaces[pFaces[i]] << endl;
                Info<< "Face " << pFaces[i] << " is in patch " << facePatch_[pFaces[i]] << endl;
                Info<< "Edge " << pEdges[i] << " has nodes " << edges[pEdges[i]] << endl;
            }
            # endif

            // find the best fitting edge to move towards the corner
            label bestEdge(-1);
            scalar bestAlignment(0.0);

            vector dv = it->second.first - points[bPoints[bpI]];
            dv /= (mag(dv) + VSMALL);

            if (facesContainingPoint[bestPoint].size() == 1)
            {
                const label bfI = facesContainingPoint[bestPoint][0];

                // only the edges of this face can be candidates
                forAll(pEdges, i)
                {
                    const label beI = pEdges[i];
                    const edge& e = edges[pEdges[i]];

                    if (!(edgeFaces(beI, 0) == bfI || edgeFaces(beI, 1) == bfI))
                        continue;

                    vector ev
                    (
                        points[e.otherVertex(bPoints[bpI])] - points[bPoints[bpI]]
                    );
                    ev /= (mag(ev) + VSMALL);

                    const scalar metric = 0.5*(1.0 + (dv & ev));

                    if (metric > bestAlignment)
                    {
                        bestAlignment = metric;
                        bestEdge = pEdges[i];
                    }
                }
            }
            else
            {
                // the edge containing the best point is the candidate
                forAll(pEdges, i)
                {
                    const edge& e = edges[pEdges[i]];

                    if (e.otherVertex(bestPoint) != -1)
                    {
                        // select the edge which contains best fitting point
                        bestEdge = pEdges[i];
                        break;
                    }
                }
            }

            # ifdef DEBUGEdgeExtractor
            Info<< "Best edge " << bestEdge
                << " has alignment " << bestAlignment << endl;
            # endif

            // find groups of sorted faces at this corner which are bounded by
            // existing feature edges and the new candidate
            DynList<label> faceInGroup(pFaces.size(), -1);
            label groupI(0);
            forAll(pFaces, i)
            {
                if (faceInGroup[i] != -1)
                    continue;

                DynList<label> front;
                front.append(i);
                faceInGroup[i] = groupI;

                while (front.size() != 0)
                {
                    const label currIndex = front.remove();

                    const label nbeI = pEdges[currIndex];
                    const label pbeI = pEdges[pFaces.rcIndex(currIndex)];
                    if ((nbeI != bestEdge) && !featureEdge.found(nbeI))
                    {
                        const label nfI = pFaces.fcIndex(currIndex);

                        if (faceInGroup[nfI] == -1)
                        {
                            faceInGroup[nfI] = groupI;
                            front.append(nfI);
                        }
                    }
                    else if ((pbeI != bestEdge) && !featureEdge.found(pbeI))
                    {
                        const label pfI = pFaces.rcIndex(currIndex);

                        if (faceInGroup[pfI] == -1)
                        {
                            faceInGroup[pfI] = groupI;
                            front.append(pfI);
                        }
                    }
                }

                ++groupI;
            }

            # ifdef DEBUGEdgeExtractor
            Info<< "faceInGroup " << faceInGroup << endl;
            # endif

            // check which group of faces shall change patch in order to make
            // the best fitting edge a feature edge
            DynList<labelPair> groupPairs, groupPatches;
            labelPair groupsForChanging(-1, -1);
            forAll(pEdges, i)
            {
                const label beI = pEdges[i];

                if (beI == bestEdge)
                {
                    label pos = pFaces.find(edgeFaces(beI, 0));
                    groupsForChanging.first() = faceInGroup[pos];
                    pos = pFaces.find(edgeFaces(beI, 1));
                    groupsForChanging.second() = faceInGroup[pos];
                }
                else if (featureEdge.found(beI))
                {
                    labelPair lp, lpp;
                    label pos = pFaces.find(edgeFaces(beI, 0));
                    lpp.first() = facePatch_[pFaces[pos]];
                    lp.first() = faceInGroup[pos];
                    pos = pFaces.find(edgeFaces(beI, 1));
                    lpp.second() = facePatch_[pFaces[pos]];
                    lp.second() = faceInGroup[pos];

                    groupPairs.append(lp);
                    groupPatches.append(lpp);
                }
            }

            # ifdef DEBUGEdgeExtractor
            Info<< "group pairs " << groupPairs << endl;
            Info<< "Group patches " << groupPatches << endl;
            Info<< "Groups for changing " << groupsForChanging << endl;
            # endif

            // check which groups shall change their patch
            // desired patches at the best edge are determined by finding
            // the best alignment between the bestEdge and the feature edges
            scalar maxAlignment(0.0);
            label bestGroupsOfPatches(-1);
            forAll(groupPatches, i)
            {
                const labelPair& pp = groupPatches[i];

                const point& bes = points[edges[bestEdge].start()];
                const point& bee = points[edges[bestEdge].end()];

                DynList<label> patches(2);
                patches[0] = pp.first();
                patches[1] = pp.second();

                point mps, mpe;
                scalar dSqS, dSqE;
                label nse;

                meshOctree_.findNearestEdgePoint(mps, dSqS, nse, bes, patches);
                meshOctree_.findNearestEdgePoint(mpe, dSqE, nse, bee, patches);

                const scalar align = 0.5*(1.0 + ((bee - bes) & (mpe - mps)));

                if (align > maxAlignment)
                {
                    maxAlignment = align;
                    bestGroupsOfPatches = i;
                }
            }

            # ifdef DEBUGEdgeExtractor
            Info<< "Best groups of patches " << bestGroupsOfPatches << endl;
            # endif

            // assign new patches
            FixedList<label, 2> newPatchForGroup(-1);
            DynList<labelPair>& newPatchForFace = facesAtCornerNewPatches[bpI];
            forAll(groupPatches, i)
            {
                if (i == bestGroupsOfPatches)
                    continue;

                const labelPair& gp = groupPairs[i];
                const labelPair& gpp = groupPatches[i];

                if
                (
                    gp.first() == groupsForChanging.first() ||
                    gp.first() == groupsForChanging.second()
                )
                {
                    const label otherPatch = gpp.second();

                    scalar Eold(0.0), Enew(0.0);
                    forAll(faceInGroup, j)
                    {
                        if (faceInGroup[j] != gp.first())
                            continue;

                        const label bfI = pFaces[j];

                        forAllRow(faceEdges, bfI, feI)
                        {
                            const label beI = faceEdges(bfI, feI);

                            # ifdef DEBUGEdgeExtractor
                            Info<< "2. beI " << beI << endl;
                            # endif

                            const point& ps = points[edges[beI].start()];
                            const point& pe = points[edges[beI].end()];
                            const scalar magE = edges[beI].mag(points) + VSMALL;

                            vector ev = pe - ps;
                            ev /= magE;

                            label nei = edgeFaces(beI, 0);
                            if (nei == bfI)
                                nei = edgeFaces(beI, 1);

                            const label posNei = pFaces.find(nei);
                            if
                            (
                                (posNei < 0) ||
                                (faceInGroup[posNei] != faceInGroup[j])
                            )
                            {
                                if (facePatch_[bfI] != facePatch_[nei])
                                {
                                    DynList<label> patches(2);
                                    patches[0] = facePatch_[bfI];
                                    patches[1] = facePatch_[nei];
                                    // calculate deformation energy
                                    // of the old state
                                    point mps, mpe;
                                    scalar dSqS, dSqE;

                                    meshOctree_.findNearestPointToPatches
                                    (
                                        mps,
                                        dSqS,
                                        ps,
                                        patches
                                    );
                                    meshOctree_.findNearestPointToPatches
                                    (
                                        mpe,
                                        dSqE,
                                        pe,
                                        patches
                                    );

                                    vector fv = mpe - mps;
                                    fv /= (mag(fv) + VSMALL);

                                    scalar c = min(fv & ev, 1.0);
                                    c = max(-1.0, c);
                                    const scalar angle = acos(c);

                                    Eold +=
                                        1.0/magE*(dSqS + dSqE) + magE*angle;
                                }
                                if (otherPatch != facePatch_[nei])
                                {
                                    // calculate deformation energy
                                    // of the new state
                                    DynList<label> patches(2);
                                    patches[0] = otherPatch;
                                    patches[1] = facePatch_[nei];

                                    point mps, mpe;
                                    scalar dSqS, dSqE;

                                    meshOctree_.findNearestPointToPatches
                                    (
                                        mps,
                                        dSqS,
                                        ps,
                                        patches
                                    );
                                    meshOctree_.findNearestPointToPatches
                                    (
                                        mpe,
                                        dSqE,
                                        pe,
                                        patches
                                    );

                                    vector fv = mpe - mps;
                                    fv /= (mag(fv) + VSMALL);

                                    scalar c = min(fv & ev, 1.0);
                                    c = max(-1.0, c);
                                    const scalar angle = acos(c);

                                    Enew +=
                                        1.0/magE*(dSqS + dSqE) + magE*angle;
                                }
                            }
                        }
                    }

                    # ifdef DEBUGEdgeExtractor
                    Info<< "1. Eold " << Eold << " Enew " << Enew << endl;
                    # endif

                    if (Enew <= Eold)
                    {
                        newPatchForGroup[0] = otherPatch;

                        forAll(faceInGroup, j)
                        {
                            if (faceInGroup[j] != gp.first())
                                continue;

                            newPatchForFace.append
                            (
                                labelPair(pFaces[j], otherPatch)
                            );
                        }
                    }
                }
                else if
                (
                    gp.second() == groupsForChanging.first() ||
                    gp.second() == groupsForChanging.second()
                )
                {
                    const label otherPatch = gpp.first();

                    scalar Eold(0.0), Enew(0.0);
                    forAll(faceInGroup, j)
                    {
                        if (faceInGroup[j] != gp.second())
                            continue;

                        const label bfI = pFaces[j];

                        forAllRow(faceEdges, bfI, feI)
                        {
                            const label beI = faceEdges(bfI, feI);

                            # ifdef DEBUGEdgeExtractor
                            Info<< "2. beI " << beI << endl;
                            # endif

                            const point& ps = points[edges[beI].start()];
                            const point& pe = points[edges[beI].end()];
                            const scalar magE = edges[beI].mag(points) + VSMALL;

                            vector ev = pe - ps;
                            ev /= magE;

                            label nei = edgeFaces(beI, 0);
                            if (nei == bfI)
                                nei = edgeFaces(beI, 1);

                            const label posNei = pFaces.find(nei);
                            if
                            (
                                (posNei < 0) ||
                                (faceInGroup[posNei] != faceInGroup[j])
                            )
                            {
                                if (facePatch_[bfI] != facePatch_[nei])
                                {
                                    DynList<label> patches(2);
                                    patches[0] = facePatch_[bfI];
                                    patches[1] = facePatch_[nei];

                                    // calculate deformation energy
                                    // of the old state
                                    point mps, mpe;
                                    scalar dSqS, dSqE;

                                    meshOctree_.findNearestPointToPatches
                                    (
                                        mps,
                                        dSqS,
                                        ps,
                                        patches
                                    );
                                    meshOctree_.findNearestPointToPatches
                                    (
                                        mpe,
                                        dSqE,
                                        pe,
                                        patches
                                    );

                                    vector fv = mpe - mps;
                                    fv /= (mag(fv) + VSMALL);

                                    scalar c = min(fv & ev, 1.0);
                                    c = max(-1.0, c);
                                    const scalar angle = acos(c);

                                    Eold +=
                                        1.0/magE*(dSqS + dSqE) + magE*angle;
                                }
                                if (otherPatch != facePatch_[nei])
                                {
                                    // calculate deformation energy
                                    // of the new state
                                    DynList<label> patches(2);
                                    patches[0] = otherPatch;
                                    patches[1] = facePatch_[nei];

                                    point mps, mpe;
                                    scalar dSqS, dSqE;

                                    meshOctree_.findNearestPointToPatches
                                    (
                                        mps,
                                        dSqS,
                                        ps,
                                        patches
                                    );
                                    meshOctree_.findNearestPointToPatches
                                    (
                                        mpe,
                                        dSqE,
                                        pe,
                                        patches
                                    );

                                    vector fv = mpe - mps;
                                    fv /= (mag(fv) + VSMALL);

                                    scalar c = min(fv & ev, 1.0);
                                    c = max(-1.0, c);
                                    const scalar angle = acos(c);

                                    Enew +=
                                        1.0/magE*(dSqS + dSqE) + magE*angle;
                                }
                            }
                        }
                    }

                    # ifdef DEBUGEdgeExtractor
                    Info<< "2. Eold " << Eold << " Enew " << Enew << endl;
                    # endif

                    if (Enew <= Eold)
                    {
                        newPatchForGroup[1] = otherPatch;

                        forAll(faceInGroup, j)
                        {
                            if (faceInGroup[j] != gp.second())
                                continue;

                            newPatchForFace.append
                            (
                                labelPair(pFaces[j], otherPatch)
                            );
                        }
                    }
                }
            }

            # ifdef DEBUGEdgeExtractor
            Info<< "New patches for boundary faces "
                 << newPatchForFace << endl;
            # endif
        }

        labelHashSet changedPatch;
        forAllConstIters(facesAtCornerNewPatches, it)
        {
            const DynList<labelPair>& lp = it->second;
            forAll(lp, i)
            {
                if (changedPatch.found(lp[i].first()))
                    FatalError << "Face " << lp[i].first()
                               << " is already modified" << abort(FatalError);

                changedPatch.insert(lp[i].first());

                newBoundaryPatches[lp[i].first()] = lp[i].second();
                ++nCorrected;
            }
        }

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

        if (nCorrected)
        {
            changed = true;
            facePatch_ = newBoundaryPatches;
        }

        if (nIteration++ > 5)
        {
            Info<< "Too many iterations" << endl;
            break;
        }

        # ifdef DEBUGEdgeExtractor
        surfPtr = surfaceWithPatches();
        surfPtr->writeSurface
        (
            "checkCornersAfterIteration_" + Foam::name(nIteration) + ".fms"
        );
        deleteDemandDrivenData(surfPtr);
        # endif

    } while (nCorrected != 0);

    return changed;
}


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