/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | cfMesh: A library for mesh generation
   \\    /   O peration     |
    \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
     \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
-------------------------------------------------------------------------------
License
    This file is part of cfMesh.

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

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

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

Description

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

#include "meshSurfaceEngine.H"
#include "demandDrivenData.H"
#include "boolList.H"
#include "helperFunctions.H"
#include "VRWGraphSMPModifier.H"
#include "labelledPoint.H"
#include "HashSet.H"

#include <map>
#include <set>

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

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

namespace Foam
{

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

void meshSurfaceEngine::calculateBoundaryFaces() const
{
    if( mesh_.boundaries().size() != 0 )
    {
        const faceListPMG& faces = mesh_.faces();
        const PtrList<boundaryPatch>& boundaries = mesh_.boundaries();

        label nBoundaryFaces(0);
        if( activePatch_ < 0 )
        {
            //- take all patches
            forAll(boundaries, patchI)
                nBoundaryFaces += boundaries[patchI].patchSize();

            boundaryFacesPtr_ =
                new faceList::subList
                (
                    faces,
                    nBoundaryFaces,
                    boundaries[0].patchStart()
                );
        }
        else if( activePatch_ < boundaries.size() )
        {
            nBoundaryFaces = boundaries[activePatch_].patchSize();

            boundaryFacesPtr_ =
                new faceList::subList
                (
                    faces,
                    nBoundaryFaces,
                    boundaries[activePatch_].patchStart()
                );
        }
        else
        {
            FatalErrorIn
            (
                "void meshSurfaceEngine::calculateBoundaryFaces() const"
            ) << "Cannot select boundary faces. Invalid patch index "
              << activePatch_ << exit(FatalError);
        }

        reduce(nBoundaryFaces, sumOp<label>());
        Info << "Found " << nBoundaryFaces << " boundary faces " << endl;
    }
    else
    {
        FatalErrorIn
        (
            "void meshSurfaceEngine::calculateBoundaryFaces() const"
        ) << "Boundary faces are not at the end of the face list!"
            << exit(FatalError);
    }
}

void meshSurfaceEngine::calculateBoundaryOwners() const
{
    const labelList& owner = mesh_.owner();

    const faceList::subList& boundaryFaces = this->boundaryFaces();

    if( !boundaryFaceOwnersPtr_ )
    boundaryFaceOwnersPtr_ = new labelList(boundaryFaces.size());

    labelList& owners = *boundaryFaceOwnersPtr_;

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

    # ifdef USE_OMP
    # pragma omp parallel for schedule(static, 1)
    # endif
    forAll(boundaryFaces, fI)
    owners[fI] = owner[start+fI];
}

void meshSurfaceEngine::calculateBoundaryNodes() const
{
    //- mark boundary points
    label pointI(0);
    if( !bppPtr_ )
        bppPtr_ = new labelList(mesh_.points().size(), -1);
    labelList& bp = *bppPtr_;

    const faceList::subList& boundaryFaces = this->boundaryFaces();

    boolList isBndPoint(bp.size(), false);

    # ifdef USE_OMP
    const label nThreads = 3 * omp_get_num_procs();
    # pragma omp parallel for num_threads(nThreads) schedule(static, 1)
    # endif
    forAll(boundaryFaces, bfI)
    {
        const face& bf = boundaryFaces[bfI];

        forAll(bf, pI)
            isBndPoint[bf[pI]] = true;
    }

    forAll(isBndPoint, pI)
    {
        if( isBndPoint[pI] )
            bp[pI] = pointI++;
    }

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

        //- exchange information with processors
        //- this is need sometimes to find all nodes at the boundary
        bool found;
        do
        {
            found = false;

            //- send bnd nodes to other processor
            forAll(procBoundaries, patchI)
            {
                const label start = procBoundaries[patchI].patchStart();
                const label end = start + procBoundaries[patchI].patchSize();

                labelLongList dts;
                labelHashSet addedPoint;
                for(label faceI=start;faceI<end;++faceI)
                {
                    const face& f = faces[faceI];
                    forAll(f, pI)
                        if( (bp[f[pI]] != -1) && !addedPoint.found(f[pI]) )
                        {
                            addedPoint.insert(f[pI]);

                            //- data is sent as follows
                            //- 1. local face label in patch
                            //- 2. local node in face
                            dts.append(faceI-start);
                            dts.append((f.size() - pI) % f.size());
                        }
                }

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

            //- receive data and update positions if needed
            forAll(procBoundaries, patchI)
            {
                const label start = procBoundaries[patchI].patchStart();
                labelList receiveData;
                IPstream fromOtherProc
                (
                    Pstream::blocking,
                    procBoundaries[patchI].neiProcNo()
                );
                fromOtherProc >> receiveData;

                label counter(0);
                while( counter < receiveData.size() )
                {
                    const label fI = receiveData[counter++];
                    const label pI = receiveData[counter++];

                    if( bp[faces[start+fI][pI]] == -1 )
                    {
                        bp[faces[start+fI][pI]] = pointI++;
                        found = true;
                    }
                }
            }

            reduce(found, maxOp<bool>());

        } while( found );
    }

    if( !boundaryPointsPtr_ )
        boundaryPointsPtr_ = new labelList();

    labelList& boundaryPoints = *boundaryPointsPtr_;
    boundaryPoints.setSize(pointI);

    //- fill the boundaryPoints list
    # ifdef USE_OMP
    # pragma omp parallel for num_threads(nThreads) schedule(static, 1)
    # endif
    forAll(bp, bpI)
    {
        if( bp[bpI] != -1 )
            boundaryPoints[bp[bpI]] = bpI;
    }
}

void meshSurfaceEngine::calculateBoundaryFacePatches() const
{
    const faceList::subList& bFaces = this->boundaryFaces();
    boundaryFacePatchPtr_ = new labelList(bFaces.size());
    labelList& facePatch = *boundaryFacePatchPtr_;

    label faceI(0);
    const PtrList<boundaryPatch>& boundaries = mesh_.boundaries();
    forAll(boundaries, patchI)
    {
        const label nFaces = boundaries[patchI].patchSize();
        for(label patchFaceI=0;patchFaceI<nFaces;++patchFaceI)
        {
            facePatch[faceI] = patchI;
            ++faceI;
        }
    }
}

void meshSurfaceEngine::calculatePointFaces() const
{
    //- fill pointFacesAddr
    if( !pointFacesPtr_ )
        pointFacesPtr_ = new VRWGraph();
    VRWGraph& pointFacesAddr = *pointFacesPtr_;

    if( !pointInFacePtr_ )
        pointInFacePtr_ = new VRWGraph();
    VRWGraph& pointInFaceAddr = *pointInFacePtr_;

    const labelList& bPoints = this->boundaryPoints();
    const faceList::subList& bFaces = this->boundaryFaces();

    //- create boundary points
    const labelList& bp = this->bp();

    labelLongList npf;

    # ifdef USE_OMP
    label nThreads = 3 * omp_get_num_procs();
    if( bPoints.size() < 1000 )
        nThreads = 1;
    # else
    const label nThreads(1);
    # endif

    label minRow(INT_MAX), maxRow(0);
    List<List<LongList<labelPair> > > dataForOtherThreads(nThreads);

    # ifdef USE_OMP
    # pragma omp parallel num_threads(nThreads)
    # endif
    {
        # ifdef USE_OMP
        const label threadI = omp_get_thread_num();
        # else
        const label threadI(0);
        # endif

        List<LongList<labelPair> >& dot = dataForOtherThreads[threadI];
        dot.setSize(nThreads);

        //- find min and max entry in the graph
        //- they are used for assigning ranges of values local for each process
        label localMinRow(minRow), localMaxRow(0);
        # ifdef USE_OMP
        # pragma omp for schedule(static)
        # endif
        forAll(bFaces, bfI)
        {
            const face& bf = bFaces[bfI];
            forAll(bf, pI)
            {
                const label bpI = bp[bf[pI]];
                localMaxRow = Foam::max(localMaxRow, bpI);
                localMinRow = Foam::min(localMinRow, bpI);
            }
        }

        ++localMaxRow;

        # ifdef USE_OMP
        # pragma omp critical
        # endif
        {
            minRow = Foam::min(minRow, localMinRow);
            minRow = Foam::max(minRow, 0);
            maxRow = Foam::max(maxRow, localMaxRow);

            npf.setSize(maxRow);
        }

        # ifdef USE_OMP
        # pragma omp barrier
        # endif

        //- initialise appearances
        # ifdef USE_OMP
        # pragma omp for schedule(static)
        # endif
        for(label i=0;i<maxRow;++i)
            npf[i] = 0;

        # ifdef USE_OMP
        # pragma omp barrier
        # endif

        const label range = (maxRow - minRow) / nThreads + 1;
        const label localMin = minRow + threadI * range;
        const label localMax = Foam::min(localMin + range, maxRow);

        //- find the number of appearances of each element in the original graph
        # ifdef USE_OMP
        # pragma omp for schedule(static)
        # endif
        forAll(bFaces, bfI)
        {
            const face& bf = bFaces[bfI];

            forAll(bf, pI)
            {
                const label bpI = bp[bf[pI]];

                const label threadNo = (bpI - minRow) / range;

                if( threadNo == threadI )
                {
                    ++npf[bpI];
                }
                else
                {
                    dot[threadNo].append(labelPair(bpI, bfI));
                }
            }
        }

        # ifdef USE_OMP
        # pragma omp barrier
        # endif

        //- count the appearances which are not local to the processor
        for(label i=0;i<nThreads;++i)
        {
            const LongList<labelPair>& data =
                dataForOtherThreads[i][threadI];

            forAll(data, j)
                ++npf[data[j].first()];
        }

        # ifdef USE_OMP
        # pragma omp barrier
        # endif

        //- allocate graph
        # ifdef USE_OMP
        # pragma omp master
        # endif
        {
            VRWGraphSMPModifier(pointFacesAddr).setSizeAndRowSize(npf);
            VRWGraphSMPModifier(pointInFaceAddr).setSizeAndRowSize(npf);
        }

        # ifdef USE_OMP
        # pragma omp barrier
        # endif

        for(label i=localMin;i<localMax;++i)
            npf[i] = 0;

        //- start filling reverse addressing graph
        //- update data from processors with smaller labels
        for(label i=0;i<threadI;++i)
        {
            const LongList<labelPair>& data =
                dataForOtherThreads[i][threadI];

            forAll(data, j)
            {
                const label bpI = data[j].first();
                const label bfI = data[j].second();

                pointFacesAddr(bpI, npf[bpI]) = bfI;
                pointInFaceAddr(bpI, npf[bpI]) =
                    bFaces[bfI].which(bPoints[bpI]);

                ++npf[bpI];
            }
        }

        //- update data local to the processor
        # ifdef USE_OMP
        # pragma omp for schedule(static)
        # endif
        forAll(bFaces, bfI)
        {
            const face& bf = bFaces[bfI];

            forAll(bf, pI)
            {
                const label bpI = bp[bf[pI]];

                if( (bpI >= localMin) && (bpI < localMax) )
                {
                    pointInFaceAddr(bpI, npf[bpI]) = pI;
                    pointFacesAddr(bpI, npf[bpI]++) = bfI;
                }
            }
        }

        //- update data from the processors with higher labels
        for(label i=threadI+1;i<nThreads;++i)
        {
            const LongList<labelPair>& data =
                dataForOtherThreads[i][threadI];

            forAll(data, j)
            {
                const label bpI = data[j].first();
                const label bfI = data[j].second();

                pointFacesAddr(bpI, npf[bpI]) = bfI;
                pointInFaceAddr(bpI, npf[bpI]) =
                    bFaces[bfI].which(bPoints[bpI]);

                ++npf[bpI];
            }
        }
    }

    pointFacesAddr.setSize(bPoints.size());
    pointInFaceAddr.setSize(bPoints.size());
}

void meshSurfaceEngine::calculatePointPatches() const
{
    if( !pointPatchesPtr_ )
        pointPatchesPtr_ = new VRWGraph();
    VRWGraph& pPatches = *pointPatchesPtr_;

    const labelList& facePatch = boundaryFacePatches();
    const VRWGraph& pFaces = pointFaces();

    # ifdef USE_OMP
    const label nThreads = 3 * omp_get_num_procs();
    # endif

    labelList npPatches(pFaces.size());

    # ifdef USE_OMP
    # pragma omp parallel num_threads(nThreads)
    # endif
    {
        # ifdef USE_OMP
        # pragma omp for schedule(static)
        # endif
        forAll(npPatches, i)
            npPatches[i] = 0;

        # ifdef USE_OMP
        # pragma omp for schedule(static)
        # endif
        forAll(pFaces, bpI)
        {
            DynList<label> pf;
            forAllRow(pFaces, bpI, pfI)
                pf.appendIfNotIn(facePatch[pFaces(bpI, pfI)]);

            npPatches[bpI] = pf.size();
        }

        # ifdef USE_OMP
        # pragma omp barrier

        # pragma omp master
        # endif
        VRWGraphSMPModifier(pPatches).setSizeAndRowSize(npPatches);

        # ifdef USE_OMP
        # pragma omp barrier

        # pragma omp for schedule(static)
        # endif
        forAll(pFaces, bpI)
        {
            DynList<label> pf;
            forAllRow(pFaces, bpI, pfI)
                pf.appendIfNotIn(facePatch[pFaces(bpI, pfI)]);

            pPatches.setRow(bpI, pf);
        }
    }

    if( Pstream::parRun() )
    {
        const labelList& globalPointLabel = globalBoundaryPointLabel();
        const VRWGraph& bpAtProcs = this->bpAtProcs();
        const Map<label>& globalToLocal = globalToLocalBndPointAddressing();

        std::map<label, labelLongList> exchangeData;

        forAllConstIter(Map<label>, globalToLocal, iter)
        {
            const label bpI = iter();

            forAllRow(bpAtProcs, bpI, procI)
            {
                const label neiProc = bpAtProcs(bpI, procI);
                if( neiProc == Pstream::myProcNo() )
                    continue;

                if( exchangeData.find(neiProc) == exchangeData.end() )
                {
                    exchangeData.insert
                    (
                        std::make_pair(neiProc, labelLongList())
                    );
                }

                labelLongList& dataToSend = exchangeData[neiProc];

                //- prepare data which will be sent
                //- data is sent as follows
                //- 1. global point label
                //- 2. number of local patches for point
                //- 3. patch labels for a given point
                dataToSend.append(globalPointLabel[bpI]);
                dataToSend.append(pPatches.sizeOfRow(bpI));
                forAllRow(pPatches, bpI, patchI)
                    dataToSend.append(pPatches(bpI, patchI));
            }
        }

        //- exchange data with other processors
        labelLongList receivedData;
        help::exchangeMap(exchangeData, receivedData);

        label counter(0);
        while( counter < receivedData.size() )
        {
            const label bpI = globalToLocal[receivedData[counter++]];
            const label nPatches = receivedData[counter++];
            for(label i=0;i<nPatches;++i)
                pPatches.appendIfNotIn(bpI, receivedData[counter++]);
        }
    }
}

void meshSurfaceEngine::calculatePointPoints() const
{
    if( !pointPointsPtr_ )
        pointPointsPtr_ = new VRWGraph();

    VRWGraph& pointPoints = *pointPointsPtr_;

    const labelList& boundaryPoints = this->boundaryPoints();
    const faceList::subList& bFaces = this->boundaryFaces();
    const VRWGraph& pFaces = this->pointFaces();
    const labelList& bp = this->bp();

    # ifdef USE_OMP
    const label nThreads = 3 * omp_get_num_procs();
    # endif

    labelList npp(boundaryPoints.size());

    # ifdef USE_OMP
    # pragma omp parallel num_threads(nThreads)
    # endif
    {
        # ifdef USE_OMP
        # pragma omp for schedule(static)
        # endif
        forAll(npp, i)
            npp[i] = 0;

        # ifdef USE_OMP
        # pragma omp for schedule(static)
        # endif
        forAll(pFaces, bpI)
        {
            DynList<label> pPoints;

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

                const label pos = bf.which(boundaryPoints[bpI]);

                pPoints.appendIfNotIn(bp[bf.nextLabel(pos)]);
                pPoints.appendIfNotIn(bp[bf.prevLabel(pos)]);
            }

            npp[bpI] = pPoints.size();
        }

        # ifdef USE_OMP
        # pragma omp barrier

        # pragma omp master
        # endif
        VRWGraphSMPModifier(pointPoints).setSizeAndRowSize(npp);

        # ifdef USE_OMP
        # pragma omp barrier

        # pragma omp for schedule(static)
        # endif
        forAll(pFaces, bpI)
        {
            DynList<label> pPoints;

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

                const label pos = bf.which(boundaryPoints[bpI]);

                pPoints.appendIfNotIn(bp[bf.nextLabel(pos)]);
                pPoints.appendIfNotIn(bp[bf.prevLabel(pos)]);
            }

            pointPoints.setRow(bpI, pPoints);
        }
    }

    if( Pstream::parRun() )
    {
        //- this is needed to make the connection matrix symmetric
        //- on all processors. In some cases the points on a given processor
        //- may not be connected because of a single layer of faces on some
        //- other processor. P0, P0 | P1 | P0 P0
        const labelList& globalPointLabel = this->globalBoundaryPointLabel();
        const Map<label>& globalToLocal =
            this->globalToLocalBndPointAddressing();
        const VRWGraph& bpAtProcs = this->bpAtProcs();

        std::map<label, labelLongList> exchangeData;
        forAllConstIter(Map<label>, globalToLocal, iter)
        {
            const label bpI = iter();

            DynList<label> neiToSend;
            forAllRow(pointPoints, bpI, j)
            {
                const label bpJ = pointPoints(bpI, j);
                if( bpAtProcs.sizeOfRow(bpJ) != 0 )
                    neiToSend.append(bpJ);
            }

            forAllRow(bpAtProcs, bpI, procI)
            {
                const label neiProc = bpAtProcs(bpI, procI);
                if( neiProc == Pstream::myProcNo() )
                    continue;

                if( exchangeData.find(neiProc) == exchangeData.end() )
                    exchangeData.insert(std::make_pair(neiProc,labelLongList()));

                if( neiToSend.size() != 0 )
                {
                    labelLongList& dts = exchangeData[neiProc];
                    dts.append(globalPointLabel[bpI]);
                    dts.append(neiToSend.size());
                    forAll(neiToSend, i)
                        dts.append(globalPointLabel[neiToSend[i]]);
                }
            }
        }

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

        label counter(0);
        while( counter < receivedData.size() )
        {
            const label bpI = globalToLocal[receivedData[counter++]];
            const label size = receivedData[counter++];
            for(label i=0;i<size;++i)
            {
                const label gpI = receivedData[counter++];
                if( globalToLocal.found(gpI) )
                    pointPoints.appendIfNotIn(bpI, globalToLocal[gpI]);
            }
        }
    }
}

void meshSurfaceEngine::calculatePointNormals() const
{
    const VRWGraph& pFaces = pointFaces();
    const vectorField& fNormals = faceNormals();

    pointNormalsPtr_ = new vectorField(pFaces.size());

    # ifdef USE_OMP
    # pragma omp parallel for if( pFaces.size() > 1000 ) schedule(dynamic, 50)
    # endif
    forAll(pFaces, pI)
    {
        vector normal(vector::zero);

        forAllRow(pFaces, pI, pfI)
            normal += fNormals[pFaces(pI, pfI)];

        const scalar d = mag(normal);
        if( d > VSMALL )
        {
            normal /= d;
        }
        else
        {
            normal = vector::zero;
        }

        (*pointNormalsPtr_)[pI] = normal;
    }

    updatePointNormalsAtProcBoundaries();
}

void meshSurfaceEngine::calculateFaceNormals() const
{
    const faceList::subList& bFaces = this->boundaryFaces();
    const pointFieldPMG& points = mesh_.points();

    faceNormalsPtr_ = new vectorField(bFaces.size());

    # ifdef USE_OMP
    # pragma omp parallel for if( bFaces.size() > 1000 )
    # endif
    forAll(bFaces, bfI)
    {
        const face& bf = bFaces[bfI];

        faceNormalsPtr_->operator[](bfI) = bf.normal(points);
    }
}

void meshSurfaceEngine::calculateFaceCentres() const
{
    const faceList::subList& bFaces = this->boundaryFaces();
    const pointFieldPMG& points = mesh_.points();

    faceCentresPtr_ = new vectorField(bFaces.size());

    # ifdef USE_OMP
    # pragma omp parallel for if( bFaces.size() > 1000 )
    # endif
    forAll(bFaces, bfI)
        faceCentresPtr_->operator[](bfI) = bFaces[bfI].centre(points);
}

void meshSurfaceEngine::updatePointNormalsAtProcBoundaries() const
{
    if( !Pstream::parRun() )
        return;

    const VRWGraph& pFaces = pointFaces();
    const vectorField& fNormals = faceNormals();
    const labelList& globalPointLabel = this->globalBoundaryPointLabel();
    const Map<label>& globalToLocal =
        this->globalToLocalBndPointAddressing();
    const VRWGraph& bpAtProcs = this->bpAtProcs();

    vectorField& pNormals = *pointNormalsPtr_;

    //- create data which will be sent to other processors
    std::map<label, LongList<labelledPoint> > exchangeData;

    forAllConstIter(Map<label>, globalToLocal, iter)
    {
        const label bpI = iter();

        vector& n = pNormals[bpI];
        n = vector::zero;

        forAllRow(pFaces, bpI, pfI)
            n += fNormals[pFaces(bpI, pfI)];

        forAllRow(bpAtProcs, bpI, procI)
        {
            const label neiProc = bpAtProcs(bpI, procI);
            if( neiProc == Pstream::myProcNo() )
                continue;
            if( exchangeData.find(neiProc) == exchangeData.end() )
            {
                exchangeData.insert
                (
                    std::make_pair(neiProc, LongList<labelledPoint>())
                );
            }

            exchangeData[neiProc].append
            (
                labelledPoint(globalPointLabel[bpI], n)
            );
        }
    }

    //- exchange data with other procs
    LongList<labelledPoint> receivedData;
    help::exchangeMap(exchangeData, receivedData);

    forAll(receivedData, i)
    {
        const label bpI = globalToLocal[receivedData[i].pointLabel()];
        pNormals[bpI] += receivedData[i].coordinates();
    }

    //- normalize vectors
    # ifdef USE_OMP
    # pragma omp parallel for if( bpAtProcs.size() > 1000 ) \
    schedule(guided)
    # endif
    forAll(bpAtProcs, bpI)
    {
        if( bpAtProcs.sizeOfRow(bpI) == 0 )
            continue;

        vector normal = pNormals[bpI];
        const scalar d = mag(normal);
        if( d > VSMALL )
        {
            normal /= d;
        }
        else
        {
            normal = vector::zero;
        }

        pNormals[bpI] = normal;
    }
}

void meshSurfaceEngine::calculateEdgesAndAddressing() const
{
    const VRWGraph& pFaces = pointFaces();
    const faceList::subList& bFaces = boundaryFaces();
    const labelList& bPoints = boundaryPoints();
    const labelList& bp = this->bp();

    edgesPtr_ = new edgeList();
    edgeList& edges = *edgesPtr_;

    bpEdgesPtr_ = new VRWGraph();
    VRWGraph& bpEdges = *bpEdgesPtr_;

    # ifdef USE_OMP
    label nThreads = 3 * omp_get_num_procs();
    if( pFaces.size() < 1000 )
        nThreads = 1;
    # else
    const label nThreads(1);
    # endif

    labelList nEdgesForThread(nThreads);

    label edgeI(0);

    # ifdef USE_OMP
    # pragma omp parallel num_threads(nThreads)
    # endif
    {
        edgeLongList edgesHelper;

        # ifdef USE_OMP
        # pragma omp for schedule(static)
        # endif
        forAll(pFaces, bpI)
        {
            std::set<std::pair<label, label> > edgesAtPoint;

            forAllRow(pFaces, bpI, pfI)
            {
                const label bfI = pFaces(bpI, pfI);
                const face& bf = bFaces[bfI];

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

                if( bp[bf.nextLabel(pos)] >= bpI )
                {
                    edgesAtPoint.insert
                    (
                        std::make_pair(bf[pos], bf.nextLabel(pos))
                    );
                }
                if( bp[bf.prevLabel(pos)] >= bpI )
                {
                    edgesAtPoint.insert
                    (
                        std::make_pair(bf[pos], bf.prevLabel(pos))
                    );
                }
            }

            std::set<std::pair<label, label> >::const_iterator it;
            for(it=edgesAtPoint.begin();it!=edgesAtPoint.end();++it)
                edgesHelper.append(edge(it->first, it->second));
        }

        //- this enables other threads to see the number of edges
        //- generated by each thread
        # ifdef USE_OMP
        const label threadI = omp_get_thread_num();
        # else
        const label threadI(0);
        # endif
        nEdgesForThread[threadI] = edgesHelper.size();

        # ifdef USE_OMP
        # pragma omp critical
        # endif
        edgeI += edgesHelper.size();

        # ifdef USE_OMP
        # pragma omp barrier

        # pragma omp master
        # endif
        edgesPtr_->setSize(edgeI);

        # ifdef USE_OMP
        # pragma omp barrier
        # endif

        //- find the starting position of the edges generated by this thread
        //- in the global list of edges
        label localStart(0);
        for(label i=0;i<threadI;++i)
            localStart += nEdgesForThread[i];

        //- store edges into the global list
        forAll(edgesHelper, i)
            edgesPtr_->operator[](localStart+i) = edgesHelper[i];
    }

    //- set the bpEdges
    VRWGraphSMPModifier(bpEdges).reverseAddressing(bp, edges);
    bpEdges.setSize(pFaces.size());

    if( !Pstream::parRun() )
        return;

    bool addEdges;
    do
    {
        addEdges = false;

        //- mark boundary edges for processors which do not contain
        //- boundary faces. This procedure is needed to identify boundary
        //- edges which are not part of any boundary face on their processor
        const faceListPMG& faces = mesh_.faces();
        const PtrList<processorBoundaryPatch>& procBoundaries =
            mesh_.procBoundaries();

        //- send boundary edges to neighbour processors
        forAll(procBoundaries, patchI)
        {
            const label start = procBoundaries[patchI].patchStart();
            const label end = start + procBoundaries[patchI].patchSize();

            labelLongList dts;
            for(label faceI=start;faceI<end;++faceI)
            {
                const face& f = faces[faceI];
                forAll(f, eI)
                {
                    const edge e = f.faceEdge(eI);
                    const label s = bp[e.start()];
                    if( s < 0 )
                        continue;

                    forAllRow(bpEdges, s, peI)
                        if( edges[bpEdges(s, peI)] == e )
                        {
                            dts.append(faceI-start);
                            dts.append((f.size()-1-eI)%f.size());
                            break;
                        }
                }
            }

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

        //- receive data from other processors. Mark edges which are not yet
        //- marked as boundary edges
        forAll(procBoundaries, patchI)
        {
            labelList receivedEdges;
            IPstream fromOtherProc
            (
                Pstream::blocking,
                procBoundaries[patchI].neiProcNo()
            );
            fromOtherProc >> receivedEdges;

            const label start = procBoundaries[patchI].patchStart();
            label nReceivedEdges(0);
            while( nReceivedEdges < receivedEdges.size() )
            {
                const face& f = faces[start+receivedEdges[nReceivedEdges++]];
                const label eI = receivedEdges[nReceivedEdges++];

                const edge e = f.faceEdge(eI);
                const label s = bp[e.start()];

                bool found(false);
                forAllRow(bpEdges, s, peI)
                    if( edges[bpEdges(s, peI)] == e )
                    {
                        found = true;
                        break;
                    }

                if( !found )
                {
                    //- create a new edge
                    addEdges = true;
                    edges.newElmt(edgeI) = e;

                    bpEdges.append(bp[e.start()], edgeI);
                    bpEdges.append(bp[e.end()], edgeI);
                    ++edgeI;
                }
            }
        }

        reduce(addEdges, maxOp<bool>());
    } while( addEdges );

    edges.setSize(edgeI);
}

void meshSurfaceEngine::calculateFaceEdgesAddressing() const
{
    const faceList::subList& bFaces = this->boundaryFaces();
    const labelList& bp = this->bp();
    const edgeList& edges = this->edges();
    const VRWGraph& bpEdges = this->boundaryPointEdges();

    faceEdgesPtr_ = new VRWGraph(bFaces.size());
    VRWGraph& faceEdges = *faceEdgesPtr_;

    labelList nfe(bFaces.size());

    # ifdef USE_OMP
    const label nThreads = 3 * omp_get_num_procs();

    # pragma omp parallel num_threads(nThreads)
    # endif
    {
        # ifdef USE_OMP
        # pragma omp for schedule(static)
        # endif
        forAll(bFaces, bfI)
            nfe[bfI] = bFaces[bfI].size();

        # ifdef USE_OMP
        # pragma omp barrier

        # pragma omp master
        {
        # endif

        VRWGraphSMPModifier(faceEdges).setSizeAndRowSize(nfe);

        # ifdef USE_OMP
        }

        # pragma omp barrier

        # pragma omp for schedule(dynamic, 100)
        # endif
        forAll(faceEdges, bfI)
        {
            const face& bf = bFaces[bfI];

            forAll(bf, eI)
            {
                const edge e = bf.faceEdge(eI);

                const label bps = bp[e.start()];

                forAllRow(bpEdges, bps, peI)
                {
                    const label beI = bpEdges(bps, peI);
                    const edge& ee = edges[beI];

                    if( e == ee )
                    {
                        faceEdges(bfI, eI) = beI;
                        break;
                    }
                }
            }
        }
    }
}

void meshSurfaceEngine::calculateEdgeFacesAddressing() const
{
    const faceList::subList& bFaces = this->boundaryFaces();
    const VRWGraph& pointFaces = this->pointFaces();
    const edgeList& edges = this->edges();
    const labelList& bp = this->bp();

    edgeFacesPtr_ = new VRWGraph();
    VRWGraph& edgeFaces = *edgeFacesPtr_;

    labelList nef(edges.size());

    # ifdef USE_OMP
    const label nThreads = 3 * omp_get_num_procs();

    # pragma omp parallel num_threads(nThreads)
    # endif
    {
        # ifdef USE_OMP
        # pragma omp for schedule(static)
        # endif
        forAll(nef, edgeI)
            nef[edgeI] = 0;

        # ifdef USE_OMP
        # pragma omp for schedule(static)
        # endif
        forAll(edges, edgeI)
        {
            const edge& ee = edges[edgeI];
            const label bpI = bp[ee.start()];

            forAllRow(pointFaces, bpI, pfI)
            {
                const label bfI = pointFaces(bpI, pfI);

                const face& bf = bFaces[bfI];

                forAll(bf, eI)
                {
                    if( bf.faceEdge(eI) == ee )
                    {
                        ++nef[edgeI];
                        break;
                    }
                }
            }
        }

        # ifdef USE_OMP
        # pragma omp barrier

        # pragma omp master
        # endif
        VRWGraphSMPModifier(edgeFaces).setSizeAndRowSize(nef);

        # ifdef USE_OMP
        # pragma omp barrier

        # pragma omp for schedule(static)
        # endif
        forAll(edges, edgeI)
        {
            const edge& ee = edges[edgeI];
            const label bpI = bp[ee.start()];

            //- find boundary faces attached to this edge
            DynList<label> eFaces;
            forAllRow(pointFaces, bpI, pfI)
            {
                const label bfI = pointFaces(bpI, pfI);

                const face& bf = bFaces[bfI];

                forAll(bf, eI)
                {
                    if( bf.faceEdge(eI) == ee )
                    {
                        eFaces.append(bfI);
                        break;
                    }
                }
            }

            //- the face that owns the edge shall be the first one in the list
            // TODO: find out whether this will be necessary
            if( eFaces.size() == 2 )
            {
                const face& bf = bFaces[eFaces[1]];

                const label pos = bf.which(ee.start());

                if( bf.nextLabel(pos) == ee.end() )
                {
                    //- this face shall be the first one in the list
                    const label helper = eFaces[0];
                    eFaces[0] = eFaces[1];
                    eFaces[1] = helper;
                }
            }

            edgeFaces.setRow(edgeI, eFaces);
        }
    }
}

void meshSurfaceEngine::calculateEdgePatchesAddressing() const
{
    edgePatchesPtr_ = new VRWGraph();
    VRWGraph& edgePatches = *edgePatchesPtr_;

    const VRWGraph& edgeFaces = this->edgeFaces();
    const labelList& facePatch = this->boundaryFacePatches();

    edgePatches.setSize(edgeFaces.size());

    forAll(edgeFaces, eI)
    {
        DynList<label> ePatches;

        forAllRow(edgeFaces, eI, i)
        {
            const label patchI = facePatch[edgeFaces(eI, i)];

            ePatches.appendIfNotIn(patchI);
        }

        edgePatches.setRow(eI, ePatches);
    }

    if( Pstream::parRun() )
    {
        const Map<label>& globalToLocal = globalToLocalBndEdgeAddressing();
        const Map<label>& otherPatch = this->otherEdgeFacePatch();

        forAllConstIter(Map<label>, globalToLocal, it)
        {
            const label beI = it();

            edgePatches.appendIfNotIn(beI, otherPatch[beI]);
        }
    }
}

void meshSurfaceEngine::calculateFaceFacesAddressing() const
{
    const VRWGraph& edgeFaces = this->edgeFaces();

    const faceList::subList& bFaces = boundaryFaces();
    faceFacesPtr_ = new VRWGraph(bFaces.size());
    VRWGraph& faceFaces = *faceFacesPtr_;

    forAll(bFaces, bfI)
        faceFaces.setRowSize(bfI, bFaces[bfI].size());

    labelList nAppearances(bFaces.size(), 0);

    forAll(edgeFaces, efI)
    {
        if( edgeFaces.sizeOfRow(efI) == 2 )
        {
            const label f0 = edgeFaces(efI, 0);
            const label f1 = edgeFaces(efI, 1);

            faceFaces(f0, nAppearances[f0]++) = f1;
            faceFaces(f1, nAppearances[f1]++) = f0;
        }
        else if( Pstream::parRun() && (edgeFaces.sizeOfRow(efI) == 1) )
        {
            const label f0 = edgeFaces(efI, 0);
            faceFaces(f0, nAppearances[f0]++) = -1;
        }
        else if( Pstream::parRun() && (edgeFaces.sizeOfRow(efI) != 0 ) )
        {
            FatalErrorIn
            (
                "void meshSurfaceEngine::calculateFaceFacesAddressing() const"
            ) << "The surface of the mesh is invalid!"
                << " The number of faces containing edge " << efI
                << " is " << edgeFaces.sizeOfRow(efI)
                << " Cannot continue" << exit(FatalError);
        }
    }
}

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

} // End namespace Foam

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