lduPrimitiveMesh.C

Go to the documentation of this file.

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

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

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

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

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

#include "lduPrimitiveMesh.H"
#include "processorLduInterface.H"
#include "processorCyclicGAMGInterface.H"
#include "edgeHashes.H"
#include "labelPair.H"
#include "processorGAMGInterface.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(lduPrimitiveMesh, 0);
}


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

void Foam::lduPrimitiveMesh::checkUpperTriangular
(
    const label size,
    const labelUList& l,
    const labelUList& u
)
{
    forAll(l, facei)
    {
        if (u[facei] < l[facei])
        {
            FatalErrorInFunction
                << "Reversed face. Problem at face " << facei
                << " l:" << l[facei] << " u:" << u[facei]
                << abort(FatalError);
        }
        if (l[facei] < 0 || u[facei] < 0 || u[facei] >= size)
        {
            FatalErrorInFunction
                << "Illegal cell label. Problem at face " << facei
                << " l:" << l[facei] << " u:" << u[facei]
                << abort(FatalError);
        }
    }

    for (label facei=1; facei < l.size(); ++facei)
    {
        if (l[facei-1] > l[facei])
        {
            FatalErrorInFunction
                << "Lower not in incremental cell order."
                << " Problem at face " << facei
                << " l:" << l[facei] << " u:" << u[facei]
                << " previous l:" << l[facei-1]
                << abort(FatalError);
        }
        else if (l[facei-1] == l[facei])
        {
            // Same cell.
            if (u[facei-1] > u[facei])
            {
                FatalErrorInFunction
                    << "Upper not in incremental cell order."
                    << " Problem at face " << facei
                    << " l:" << l[facei] << " u:" << u[facei]
                    << " previous u:" << u[facei-1]
                    << abort(FatalError);
            }
        }
    }
}


Foam::labelListList Foam::lduPrimitiveMesh::globalCellCells
(
    const lduMesh& mesh,
    const globalIndex& globalNumbering
)
{
    const lduAddressing& addr = mesh.lduAddr();
    lduInterfacePtrsList interfaces = mesh.interfaces();

    const label myProci = UPstream::myProcNo(mesh.comm());

    const labelList globalIndices
    (
        Foam::identity
        (
            addr.size(),
            globalNumbering.localStart(myProci)
        )
    );

    // Get the interface cells
    PtrList<labelList> nbrGlobalCells(interfaces.size());
    {
        const label startOfRequests = UPstream::nRequests();

        // Initialise transfer of restrict addressing on the interface
        forAll(interfaces, inti)
        {
            if (interfaces.set(inti))
            {
                interfaces[inti].initInternalFieldTransfer
                (
                    Pstream::commsTypes::nonBlocking,
                    globalIndices
                );
            }
        }

        UPstream::waitRequests(startOfRequests);

        forAll(interfaces, inti)
        {
            if (interfaces.set(inti))
            {
                nbrGlobalCells.set
                (
                    inti,
                    new labelList
                    (
                        interfaces[inti].internalFieldTransfer
                        (
                            Pstream::commsTypes::nonBlocking,
                            globalIndices
                        )
                    )
                );
            }
        }
    }


    // Scan the neighbour list to find out how many times the cell
    // appears as a neighbour of the face. Done this way to avoid guessing
    // and resizing list
    labelList nNbrs(addr.size(), Zero);

    const labelUList& nbr = addr.upperAddr();
    const labelUList& own = addr.lowerAddr();

    {
        forAll(nbr, facei)
        {
            nNbrs[nbr[facei]]++;
            nNbrs[own[facei]]++;
        }

        forAll(interfaces, inti)
        {
            if (interfaces.set(inti))
            {
                for (const label celli : interfaces[inti].faceCells())
                {
                    nNbrs[celli]++;
                }
            }
        }
    }


    // Create cell-cells addressing
    labelListList cellCells(addr.size());

    forAll(cellCells, celli)
    {
        cellCells[celli].setSize(nNbrs[celli], -1);
    }

    // Reset the list of number of neighbours to zero
    nNbrs = 0;

    // Scatter the neighbour faces
    forAll(nbr, facei)
    {
        const label c0 = own[facei];
        const label c1 = nbr[facei];

        cellCells[c0][nNbrs[c0]++] = globalIndices[c1];
        cellCells[c1][nNbrs[c1]++] = globalIndices[c0];
    }
    forAll(interfaces, inti)
    {
        if (interfaces.set(inti))
        {
            const labelUList& faceCells = interfaces[inti].faceCells();

            forAll(faceCells, facei)
            {
                const label c0 = faceCells[facei];
                cellCells[c0][nNbrs[c0]++] = nbrGlobalCells[inti][facei];
            }
        }
    }

    return cellCells;
}


Foam::label Foam::lduPrimitiveMesh::totalSize
(
    const PtrList<lduPrimitiveMesh>& meshes
)
{
    label size = 0;

    for (const lduPrimitiveMesh& msh : meshes)
    {
        size += msh.lduAddr().size();
    }
    return size;
}


Foam::labelList Foam::lduPrimitiveMesh::upperTriOrder
(
    const label nCells,
    const labelUList& lower,
    const labelUList& upper
)
{
    labelList nNbrs(nCells, Zero);

    // Count number of upper neighbours
    forAll(lower, facei)
    {
        if (upper[facei] < lower[facei])
        {
            FatalErrorInFunction
                << "Problem at face:" << facei
                << " lower:" << lower[facei]
                << " upper:" << upper[facei]
                << exit(FatalError);
        }
        nNbrs[lower[facei]]++;
    }

    // Construct cell-upper cell addressing
    labelList offsets(nCells+1);
    offsets[0] = 0;
    forAll(nNbrs, celli)
    {
        offsets[celli+1] = offsets[celli]+nNbrs[celli];
    }

    nNbrs = offsets;

    labelList cellToFaces(offsets.last());
    forAll(upper, facei)
    {
        const label celli = lower[facei];
        cellToFaces[nNbrs[celli]++] = facei;
    }

    // Sort

    labelList oldToNew(lower.size());

    DynamicList<label> order;
    DynamicList<label> nbr;

    label newFacei = 0;

    for (label celli = 0; celli < nCells; ++celli)
    {
        const label startOfCell = offsets[celli];
        const label nNbr = offsets[celli+1] - startOfCell;

        nbr.resize(nNbr);
        order.resize(nNbr);

        forAll(nbr, i)
        {
            nbr[i] = upper[cellToFaces[offsets[celli]+i]];
        }
        sortedOrder(nbr, order);

        for (const label index : order)
        {
            oldToNew[cellToFaces[startOfCell + index]] = newFacei++;
        }
    }

    return oldToNew;
}


Foam::label Foam::lduPrimitiveMesh::findConnectedInterface
(
    const lduMesh& myMesh,
    const PtrList<lduPrimitiveMesh>& otherMeshes,
    const labelPairList& procAndInterfaces,

    const label nbrProci,
    const label myRank
) const
{
    // Find mesh, interfacei in procAndInterfaces
    label nbrInti = -1;
    for (const auto& procAndInterface : procAndInterfaces)
    {
        const label proci = procAndInterface[0];

        if (proci == nbrProci)
        {
            const label interfacei = procAndInterface[1];
            const lduInterfacePtrsList interfaces =
                mesh
                (
                    myMesh,
                    otherMeshes,
                    proci
                ).interfaces();
            const processorLduInterface& pldui =
                refCast<const processorLduInterface>
                (
                    interfaces[interfacei]
                );

            if (pldui.neighbProcNo() == myRank)
            {
                nbrInti = procAndInterface[1];
                break;
            }
        }
    }


    if (nbrInti == -1)
    {
        FatalErrorInFunction
            << "procAndInterfaces:" << procAndInterfaces << abort(FatalError);
    }
    return nbrInti;
}


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

Foam::lduPrimitiveMesh::lduPrimitiveMesh
(
    const label nCells,
    labelList& l,
    labelList& u,
    const label comm,
    bool reuse
)
:
    lduAddressing(nCells),
    lowerAddr_(l, reuse),
    upperAddr_(u, reuse),
    comm_(comm)
{
    if (debug && lowerAddr_.size())
    {
        if (max(lowerAddr_) >= nCells || min(lowerAddr_) < 0)
        {
            FatalErrorInFunction << "Illegal lower addressing."
                << " nCells:" << nCells
                << " max(lower):" << max(lowerAddr_)
                << " min(lower):" << min(lowerAddr_)
                << exit(FatalError);
        }
    }
    if (debug && upperAddr_.size())
    {
        if (max(upperAddr_) >= nCells || min(upperAddr_) < 0)
        {
            FatalErrorInFunction << "Illegal upper addressing."
                << " nCells:" << nCells
                << " max(upper):" << max(upperAddr_)
                << " min(upper):" << min(upperAddr_)
                << exit(FatalError);
        }
    }
}


void Foam::lduPrimitiveMesh::addInterfaces
(
    lduInterfacePtrsList& interfaces,
    const lduSchedule& ps
)
{
    interfaces_ = interfaces;
    patchSchedule_ = ps;

    // Create interfaces
    primitiveInterfaces_.setSize(interfaces_.size());
    forAll(interfaces_, i)
    {
        if (interfaces_.set(i))
        {
            primitiveInterfaces_.set(i, &interfaces_[i]);
        }
    }
}


Foam::lduPrimitiveMesh::lduPrimitiveMesh
(
    const label nCells
)
:
    lduAddressing(nCells),
    comm_(0)
{}


Foam::lduPrimitiveMesh::lduPrimitiveMesh
(
    const label nCells,
    labelList& l,
    labelList& u,
    PtrList<const lduInterface>& primitiveInterfaces,
    const lduSchedule& ps,
    const label comm
)
:
    lduAddressing(nCells),
    lowerAddr_(l, true),
    upperAddr_(u, true),
    primitiveInterfaces_(),
    patchSchedule_(ps),
    comm_(comm)
{
    if (debug && lowerAddr_.size())
    {
        if (max(lowerAddr_) >= nCells || min(lowerAddr_) < 0)
        {
            FatalErrorInFunction << "Illegal lower addressing."
                << " nCells:" << nCells
                << " max(lower):" << max(lowerAddr_)
                << " min(lower):" << min(lowerAddr_)
                << exit(FatalError);
        }
    }
    if (debug && upperAddr_.size())
    {
        if (max(upperAddr_) >= nCells || min(upperAddr_) < 0)
        {
            FatalErrorInFunction << "Illegal upper addressing."
                << " nCells:" << nCells
                << " max(upper):" << max(upperAddr_)
                << " min(upper):" << min(upperAddr_)
                << exit(FatalError);
        }
    }


    primitiveInterfaces_.transfer(primitiveInterfaces);

    // Create interfaces
    interfaces_.setSize(primitiveInterfaces_.size());
    forAll(primitiveInterfaces_, i)
    {
        if (primitiveInterfaces_.set(i))
        {
            interfaces_.set(i, &primitiveInterfaces_[i]);
        }
    }
}


Foam::lduPrimitiveMesh::lduPrimitiveMesh
(
    const label comm,
    const labelList& procAgglomMap,

    const labelList& procIDs,
    const lduMesh& myMesh,
    const PtrList<lduPrimitiveMesh>& otherMeshes,

    labelList& cellOffsets,
    labelList& faceOffsets,
    labelListList& faceMap,
    labelListList& boundaryMap,
    labelListListList& boundaryFaceMap
)
:
    lduAddressing(myMesh.lduAddr().size() + totalSize(otherMeshes)),
    comm_(comm)
{
    const label currentComm = myMesh.comm();

    forAll(otherMeshes, i)
    {
        if (otherMeshes[i].comm() != currentComm)
        {
            WarningInFunction
                << "Communicator " << otherMeshes[i].comm()
                << " at index " << i
                << " differs from that of predecessor "
                << currentComm
                << endl;    //exit(FatalError);
        }
    }

    const label nMeshes = otherMeshes.size()+1;

    const label myAgglom = procAgglomMap[UPstream::myProcNo(currentComm)];
    if (lduPrimitiveMesh::debug)
    {
        Pout<< "I am " << UPstream::myProcNo(currentComm)
            << " agglomerating into " << myAgglom
            << " as are " << findIndices(procAgglomMap, myAgglom)
            << endl;
    }

    const lduInterfacePtrsList myInterfaces = myMesh.interfaces();

    forAll(procIDs, i)
    {
        if (procAgglomMap[procIDs[i]] != procAgglomMap[procIDs[0]])
        {
            FatalErrorInFunction
                << "Processor " << procIDs[i]
                << " agglomerates to " << procAgglomMap[procIDs[i]]
                << " whereas other processors " << procIDs
                << " agglomerate to "
                << labelUIndList(procAgglomMap, procIDs)
                << exit(FatalError);
        }
    }


    // Cells get added in order.
    cellOffsets.setSize(nMeshes+1);
    cellOffsets[0] = 0;
    for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
    {
        const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);

        cellOffsets[procMeshI+1] =
            cellOffsets[procMeshI]
          + procMesh.lduAddr().size();
    }


    // Faces initially get added in order, sorted later
    labelList internalFaceOffsets(nMeshes+1);
    internalFaceOffsets[0] = 0;
    for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
    {
        const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);

        internalFaceOffsets[procMeshI+1] =
            internalFaceOffsets[procMeshI]
          + procMesh.lduAddr().lowerAddr().size();
    }

    // Count how faces get added. Proc interfaces inbetween get merged.

    // Merged proc interfaces: map from two coarse processors back to
    // - procMeshes
    // - interface in procMesh
    // (estimate size from number of patches of mesh0)
    EdgeMap<labelPairList> mergedMap(2*myInterfaces.size());

    // Unmerged proc interfaces: map from two coarse processors back to
    // - procMeshes
    // - interface in procMesh
    EdgeMap<labelPairList> unmergedMap(mergedMap.size());

    // (unmerged) global interfaces. These are present on all processors
    // in the same order (and keep the order in the merged mesh)
    List<DynamicList<label>> procToGlobal(nMeshes);


    boundaryMap.setSize(nMeshes);
    boundaryFaceMap.setSize(nMeshes);


    labelList nCoupledFaces(nMeshes, Zero);

    for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
    {
        const lduInterfacePtrsList interfaces =
            mesh(myMesh, otherMeshes, procMeshI).interfaces();

        // Initialise all boundaries as merged
        boundaryMap[procMeshI].setSize(interfaces.size(), -1);
        boundaryFaceMap[procMeshI].setSize(interfaces.size());

        // Get sorted order of processors
        forAll(interfaces, intI)
        {
            if (interfaces.set(intI))
            {
                const lduInterface& ldui = interfaces[intI];

                if (isA<processorLduInterface>(ldui))
                {
                    if (isA<processorCyclicGAMGInterface>(ldui))
                    {
                        FatalErrorInFunction
                            << "At mesh from processor " << procIDs[procMeshI]
                            << " have interface " << intI
                            << " of unhandled type " << ldui.type()
                            << ". Adapt decomposition to avoid these"
                            << exit(FatalError);
                    }

                    const processorLduInterface& pldui =
                        refCast<const processorLduInterface>(ldui);

                    label agglom0 = procAgglomMap[pldui.myProcNo()];
                    label agglom1 = procAgglomMap[pldui.neighbProcNo()];

                    const edge procEdge(agglom0, agglom1);

                    if (agglom0 != myAgglom && agglom1 != myAgglom)
                    {
                        FatalErrorInFunction
                            << "At mesh from processor " << procIDs[procMeshI]
                            << " have interface " << intI
                            << " with myProcNo:" << pldui.myProcNo()
                            << " with neighbProcNo:" << pldui.neighbProcNo()
                            << exit(FatalError);
                    }
                    else if (agglom0 == myAgglom && agglom1 == myAgglom)
                    {
                        // Merged interface
                        if (debug)
                        {
                            Pout<< "merged proc interface: myProcNo:"
                                << pldui.myProcNo()
                                << " nbr:" << pldui.neighbProcNo()
                                << " size:" << ldui.faceCells().size()
                                << endl;
                        }

                        const label nbrProcMeshI =
                            procIDs.find(pldui.neighbProcNo());

                        if (procMeshI < nbrProcMeshI)
                        {
                            // I am 'master' since get lowest numbered cells
                            nCoupledFaces[procMeshI] += ldui.faceCells().size();
                        }

                        mergedMap(procEdge).append
                        (
                            labelPair(procMeshI, intI)
                        );
                    }
                    else
                    {
                        if (debug)
                        {
                            Pout<< "external proc interface: myProcNo:"
                                << pldui.myProcNo()
                                << " nbr:" << pldui.neighbProcNo()
                                << " size:" << ldui.faceCells().size()
                                << endl;
                        }

                        unmergedMap(procEdge).append
                        (
                            labelPair(procMeshI, intI)
                        );
                    }
                }
                else
                {
                    // Still external (non proc) interface
                    procToGlobal[procMeshI].append(intI);
                }
            }
        }
    }


    // Check that all processors have any global patches in the same order
    if (debug)
    {
        const auto& global0 = procToGlobal[0];
        for (label procMeshI = 1; procMeshI < nMeshes; ++procMeshI)
        {
            const auto& global = procToGlobal[procMeshI];
            if (global != global0)
            {
                FatalErrorInFunction
                    << "At mesh from processor " << procIDs[procMeshI]
                    << " have global interfaces " << global
                    << " which differ from those on processor "
                    << procIDs[procMeshI]
                    << " : " << global0
                    << exit(FatalError);
            }
        }
    }


    if (debug)
    {
        Pout<< "Global interfaces:" << endl;
        const auto& global0 = procToGlobal[0];
        for (const label intI : global0)
        {
            Pout<< "        interfacei:" << intI
                << " type:" << myInterfaces[intI].type()
                << endl;

        }
        Pout<< endl;

        Pout<< "Remaining interfaces:" << endl;
        for (const auto& iter : unmergedMap.csorted())
        {
            Pout<< "    agglom procEdge:" << iter.key() << endl;
            const labelPairList& elems = iter.val();
            forAll(elems, i)
            {
                label procMeshI = elems[i][0];
                label interfacei = elems[i][1];
                const lduInterfacePtrsList interfaces =
                    mesh(myMesh, otherMeshes, procMeshI).interfaces();

                const processorLduInterface& pldui =
                    refCast<const processorLduInterface>
                    (
                        interfaces[interfacei]
                    );

                Pout<< "        proc:" << procIDs[procMeshI]
                    << " interfacei:" << interfacei
                    << " between:" << pldui.myProcNo()
                    << " and:" << pldui.neighbProcNo()
                    << " localsize:"
                    << interfaces[interfacei].faceCells().size()
                    << endl;
            }
            Pout<< endl;
        }
        Pout<< endl;

        Pout<< "Merged interfaces:" << endl;
        //forAllConstIters(mergedMap, iter)
        for (const auto& iter : mergedMap.csorted())
        {
            Pout<< "    agglom procEdge:" << iter.key() << endl;
            const labelPairList& elems = iter.val();

            forAll(elems, i)
            {
                label procMeshI = elems[i][0];
                label interfacei = elems[i][1];
                const lduInterfacePtrsList interfaces =
                    mesh(myMesh, otherMeshes, procMeshI).interfaces();
                const processorLduInterface& pldui =
                    refCast<const processorLduInterface>
                    (
                        interfaces[interfacei]
                    );

                Pout<< "        proc:" << procIDs[procMeshI]
                    << " interfacei:" << interfacei
                    << " between:" << pldui.myProcNo()
                    << " and:" << pldui.neighbProcNo()
                    << " localsize:"
                    << interfaces[interfacei].faceCells().size()
                    << endl;
            }
            Pout<< endl;
        }
        Pout<< endl;
    }


    // Adapt faceOffsets for internal interfaces
    faceOffsets.setSize(nMeshes+1);
    faceOffsets[0] = 0;
    faceMap.setSize(nMeshes);
    for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
    {
        const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
        label nInternal = procMesh.lduAddr().lowerAddr().size();

        faceOffsets[procMeshI+1] =
            faceOffsets[procMeshI]
          + nInternal
          + nCoupledFaces[procMeshI];

        labelList& map = faceMap[procMeshI];
        map.setSize(nInternal);
        forAll(map, i)
        {
            map[i] = faceOffsets[procMeshI] + i;
        }
    }


    // Combine upper and lower
    lowerAddr_.setSize(faceOffsets.last(), -1);
    upperAddr_.setSize(lowerAddr_.size(), -1);


    // Old internal faces and resolved coupled interfaces
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
    {
        const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);

        const labelUList& l = procMesh.lduAddr().lowerAddr();
        const labelUList& u = procMesh.lduAddr().upperAddr();

        // Add internal faces
        label allFacei = faceOffsets[procMeshI];

        forAll(l, facei)
        {
            lowerAddr_[allFacei] = cellOffsets[procMeshI]+l[facei];
            upperAddr_[allFacei] = cellOffsets[procMeshI]+u[facei];
            allFacei++;
        }


        // Add merged interfaces
        const lduInterfacePtrsList interfaces = procMesh.interfaces();

        forAll(interfaces, intI)
        {
            if (interfaces.set(intI))
            {
                const lduInterface& ldui = interfaces[intI];

                if (isA<processorLduInterface>(ldui))
                {
                    const processorLduInterface& pldui =
                        refCast<const processorLduInterface>(ldui);

                    // Look up corresponding interfaces
                    label myP = pldui.myProcNo();
                    label nbrP = pldui.neighbProcNo();
                    label nbrProcMeshI = procIDs.find(nbrP);

                    if (procMeshI < nbrProcMeshI)
                    {
                        // I am 'master' since my cell numbers will be lower
                        // since cells get added in procMeshI order.

                        const label agglom0 = procAgglomMap[myP];
                        const label agglom1 = procAgglomMap[nbrP];

                        const auto fnd =
                            mergedMap.cfind(edge(agglom0, agglom1));

                        if (fnd.good())
                        {
                            const labelPairList& elems = fnd();

                            const label nbrIntI = findConnectedInterface
                            (
                                myMesh,
                                otherMeshes,
                                elems,
                                nbrProcMeshI,
                                procIDs[procMeshI]
                            );

                            const lduInterfacePtrsList nbrInterfaces = mesh
                            (
                                myMesh,
                                otherMeshes,
                                nbrProcMeshI
                            ).interfaces();


                            const labelUList& faceCells =
                                ldui.faceCells();
                            const labelUList& nbrFaceCells =
                                nbrInterfaces[nbrIntI].faceCells();

                            if (faceCells.size() != nbrFaceCells.size())
                            {
                                FatalErrorInFunction
                                    << "faceCells:" << faceCells
                                    << " nbrFaceCells:" << nbrFaceCells
                                    << abort(FatalError);
                            }


                            labelList& bfMap =
                                boundaryFaceMap[procMeshI][intI];
                            labelList& nbrBfMap =
                                boundaryFaceMap[nbrProcMeshI][nbrIntI];

                            bfMap.setSize(faceCells.size());
                            nbrBfMap.setSize(faceCells.size());

                            forAll(faceCells, pfI)
                            {
                                lowerAddr_[allFacei] =
                                    cellOffsets[procMeshI]+faceCells[pfI];
                                bfMap[pfI] = allFacei;
                                upperAddr_[allFacei] =
                                    cellOffsets[nbrProcMeshI]+nbrFaceCells[pfI];
                                nbrBfMap[pfI] = (-allFacei-1);
                                allFacei++;
                            }
                        }
                    }
                }
            }
        }
    }


    // Sort upper-tri order
    {
        labelList oldToNew
        (
            upperTriOrder
            (
                cellOffsets.last(), //nCells
                lowerAddr_,
                upperAddr_
            )
        );

        forAll(faceMap, procMeshI)
        {
            labelList& map = faceMap[procMeshI];
            forAll(map, i)
            {
                if (map[i] >= 0)
                {
                    map[i] = oldToNew[map[i]];
                }
                else
                {
                    label allFacei = -map[i]-1;
                    map[i] = -oldToNew[allFacei]-1;
                }
            }
        }


        inplaceReorder(oldToNew, lowerAddr_);
        inplaceReorder(oldToNew, upperAddr_);

        forAll(boundaryFaceMap, proci)
        {
            const labelList& bMap = boundaryMap[proci];
            forAll(bMap, intI)
            {
                if (bMap[intI] == -1)
                {
                    // Merged interface
                    labelList& bfMap = boundaryFaceMap[proci][intI];

                    forAll(bfMap, i)
                    {
                        if (bfMap[i] >= 0)
                        {
                            bfMap[i] = oldToNew[bfMap[i]];
                        }
                        else
                        {
                            label allFacei = -bfMap[i]-1;
                            bfMap[i] = (-oldToNew[allFacei]-1);
                        }
                    }
                }
            }
        }
    }


    // Kept interfaces
    // ~~~~~~~~~~~~~~~

    interfaces_.setSize(unmergedMap.size() + procToGlobal[0].size());
    primitiveInterfaces_.setSize(interfaces_.size());

    label allInterfacei = 0;


    // Global interfaces
    // ~~~~~~~~~~~~~~~~~
    // (e.g. cyclicAMI)

    {
        const auto& global0 = procToGlobal[0];

        for (const label interfacei : global0)
        {
            //Pout<< "        interfacei:" << interfacei
            //    << " type:" << interfaces[interfacei].type()
            //    << endl;

            // Just add all individual face-cells in processor order

            // Count
            label n = 0;
            for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
            {
                const auto& procMesh = mesh(myMesh, otherMeshes, procMeshI);
                n += procMesh.interfaces()[interfacei].faceCells().size();
            }

            // Size
            labelField allFaceCells(n);
            labelField allFaceRestrictAddressing(n);
            labelList faceOffsets(nMeshes+1, 0);
            lduInterfacePtrsList allProcInterfaces(nMeshes);
            n = 0;

            // Fill
            for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
            {
                faceOffsets[procMeshI] = n;

                const auto& procMesh = mesh(myMesh, otherMeshes, procMeshI);
                const lduInterfacePtrsList interfaces = procMesh.interfaces();

                allProcInterfaces.set(procMeshI, &interfaces[interfacei]);
                boundaryMap[procMeshI][interfacei] = allInterfacei;
                labelList& bfMap = boundaryFaceMap[procMeshI][interfacei];

                const labelUList& l = interfaces[interfacei].faceCells();
                bfMap.setSize(l.size());

                forAll(l, facei)
                {
                    allFaceCells[n] = cellOffsets[procMeshI]+l[facei];
                    allFaceRestrictAddressing[n] = n;
                    bfMap[facei] = n;
                    n++;
                }
            }

            // For convenience populate last element with size
            faceOffsets.last() = n;


            const auto& myFineInterface =
                refCast<const GAMGInterface>(myInterfaces[interfacei]);

            autoPtr<GAMGInterface> ppPtr
            (
                myFineInterface.clone
                (
                    allInterfacei,
                    interfaces_,
                    global0,
                    allFaceCells,
                    allFaceRestrictAddressing,
                    faceOffsets,
                    allProcInterfaces,
                    comm_,
                    myAgglom,
                    procAgglomMap
                )
            );

            primitiveInterfaces_.set
            (
                allInterfacei,
                ppPtr.ptr()
            );

            interfaces_.set(allInterfacei, &primitiveInterfaces_[allInterfacei]);

            if (debug)
            {
                Pout<< "Created " << interfaces_[allInterfacei].type()
                    << " interface at " << allInterfacei
                    << " comm:" << comm_
                    << " myProcNo:" << myAgglom
                    << " nFaces:" << allFaceCells.size()
                    << endl;
            }

            allInterfacei++;
        }
    }


    // Unmerged processor interfaces
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    //forAllConstIters(unmergedMap, iter)
    for (const auto& iter : unmergedMap.csorted())
    {
        const labelPairList& elems = iter.val();

        // Sort processors in increasing order so both sides walk through in
        // same order.
        labelPairList procPairs(elems.size());
        forAll(elems, i)
        {
            const labelPair& elem = elems[i];
            label procMeshI = elem[0];
            label interfacei = elem[1];
            const lduInterfacePtrsList interfaces = mesh
            (
                myMesh,
                otherMeshes,
                procMeshI
            ).interfaces();

            const processorLduInterface& pldui =
                refCast<const processorLduInterface>
                (
                    interfaces[interfacei]
                );
            label myProcNo = pldui.myProcNo();
            label nbrProcNo = pldui.neighbProcNo();
            procPairs[i] = labelPair
            (
                min(myProcNo, nbrProcNo),
                max(myProcNo, nbrProcNo)
            );
        }

        labelList order(sortedOrder(procPairs));

        // Count
        label n = 0;

        forAll(order, i)
        {
            const labelPair& elem = elems[order[i]];
            label procMeshI = elem[0];
            label interfacei = elem[1];
            const lduInterfacePtrsList interfaces = mesh
            (
                myMesh,
                otherMeshes,
                procMeshI
            ).interfaces();

            n += interfaces[interfacei].faceCells().size();
        }

        // Size
        labelField allFaceCells(n);
        labelField allFaceRestrictAddressing(n);
        n = 0;

        // Fill
        forAll(order, i)
        {
            const labelPair& elem = elems[order[i]];
            label procMeshI = elem[0];
            label interfacei = elem[1];
            const lduInterfacePtrsList interfaces = mesh
            (
                myMesh,
                otherMeshes,
                procMeshI
            ).interfaces();

            boundaryMap[procMeshI][interfacei] = allInterfacei;
            labelList& bfMap = boundaryFaceMap[procMeshI][interfacei];

            const labelUList& l = interfaces[interfacei].faceCells();
            bfMap.setSize(l.size());

            forAll(l, facei)
            {
                allFaceCells[n] = cellOffsets[procMeshI]+l[facei];
                allFaceRestrictAddressing[n] = n;
                bfMap[facei] = n;
                n++;
            }
        }


        // Find out local and remote processor in new communicator

        label neighbProcNo = -1;

        // See what the two processors map onto

        if (iter.key()[0] == myAgglom)
        {
            if (iter.key()[1] == myAgglom)
            {
                FatalErrorInFunction
                    << "problem procEdge:" << iter.key()
                    << exit(FatalError);
            }

            neighbProcNo = iter.key()[1];
        }
        else
        {
            if (iter.key()[1] != myAgglom)
            {
                FatalErrorInFunction
                    << "problem procEdge:" << iter.key()
                    << exit(FatalError);
            }

            neighbProcNo = iter.key()[0];
        }

        primitiveInterfaces_.set
        (
            allInterfacei,
            new processorGAMGInterface
            (
                allInterfacei,
                interfaces_,
                allFaceCells,
                allFaceRestrictAddressing,
                comm_,
                myAgglom,
                neighbProcNo,
                tensorField(),          // forwardT
                Pstream::msgType()      // tag
            )
        );
        interfaces_.set(allInterfacei, &primitiveInterfaces_[allInterfacei]);

        if (debug)
        {
            Pout<< "Created " << interfaces_[allInterfacei].type()
                << " interface at " << allInterfacei
                << " comm:" << comm_
                << " myProcNo:" << myAgglom
                << " neighbProcNo:" << neighbProcNo
                << " nFaces:" << allFaceCells.size()
                << endl;
        }


        allInterfacei++;
    }


    if (allInterfacei != interfaces_.size())
    {
        FatalErrorInFunction << "allInterfacei:" << allInterfacei
            << " interfaces_:" << interfaces_.size() << exit(FatalError);
    }


    if (debug)
    {
        Pout<< endl
            << "Created new lduPrimitiveMesh:" << nl
            << "    cells:" << this->lduAddr().size() << nl
            << "    internal face lower:"
            << this->lduAddr().lowerAddr().size() << nl
            << "    internal faces upper:"
            << this->lduAddr().upperAddr().size()
            << endl;
        forAll(interfaces_, i)
        {
            if (interfaces_.set(i))
            {
                Pout<< "    interface:" << i << " type:" << interfaces_[i].type()
                    << nl
                    << "    faceCells:" << flatOutput(interfaces_[i].faceCells())
                    << endl;
            }
        }


        Pout<< "Original input meshes:" << endl;
        forAll(boundaryMap, procMeshI)
        {
            const auto& procMesh = mesh(myMesh, otherMeshes, procMeshI);
            const lduInterfacePtrsList interfaces = procMesh.interfaces();

            Pout<< "    proc:" << procMeshI
                << " interfaces:" << interfaces.size() << endl;
            forAll(interfaces, inti)
            {
                if (interfaces.set(inti))
                {
                    Pout<< "        int:" << inti
                        << " type:" << interfaces[inti].type()
                        << " size:" << interfaces[inti].faceCells().size()
                        << " maps to:" << boundaryMap[procMeshI][inti]
                        << endl;
                }
            }
        }
    }


    patchSchedule_ = nonBlockingSchedule<processorGAMGInterface>(interfaces_);

    if (debug)
    {
        checkUpperTriangular(cellOffsets.last(), lowerAddr_, upperAddr_);
    }
}


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

const Foam::lduMesh& Foam::lduPrimitiveMesh::mesh
(
    const lduMesh& myMesh,
    const PtrList<lduPrimitiveMesh>& otherMeshes,
    const label meshI
)
{
    return (meshI == 0 ? myMesh : otherMeshes[meshI-1]);
}


void Foam::lduPrimitiveMesh::gather
(
    const label comm,
    const lduMesh& mesh,
    PtrList<lduPrimitiveMesh>& otherMeshes
)
{
    // Force calculation of schedule (since does parallel comms)
    (void)mesh.lduAddr().patchSchedule();

    // Use PstreamBuffers
    PstreamBuffers pBufs(comm);

    // Send to master
    if (!Pstream::master(comm))
    {
        const lduAddressing& addressing = mesh.lduAddr();

        lduInterfacePtrsList interfaces(mesh.interfaces());
        boolList validInterface(interfaces.size());
        forAll(interfaces, intI)
        {
            validInterface[intI] = interfaces.set(intI);
        }

        UOPstream toMaster(UPstream::masterNo(), pBufs);

        toMaster
            << addressing.size()
            << addressing.lowerAddr()
            << addressing.upperAddr()
            << validInterface;

        forAll(interfaces, intI)
        {
            if (interfaces.set(intI))
            {
                const GAMGInterface& interface = refCast<const GAMGInterface>
                (
                    interfaces[intI]
                );

                toMaster << interface.type();
                interface.write(toMaster);
            }
        }
    }

    // Wait for finish
    pBufs.finishedGathers();

    // Consume
    if (Pstream::master(comm))
    {
        const label nProcs = UPstream::nProcs(comm);

        // Master.
        otherMeshes.setSize(nProcs-1);

        for (const int proci : UPstream::subProcs(comm))
        {
            UIPstream fromProc(proci, pBufs);

            const label nCells = readLabel(fromProc);
            labelList lowerAddr(fromProc);
            labelList upperAddr(fromProc);
            const boolList validInterface(fromProc);


            // Construct mesh without interfaces
            otherMeshes.set
            (
                proci-1,
                new lduPrimitiveMesh
                (
                    nCells,
                    lowerAddr,
                    upperAddr,
                    mesh.comm(),
                    true    // reuse
                )
            );

            // Construct GAMGInterfaces
            lduInterfacePtrsList newInterfaces(validInterface.size());
            forAll(validInterface, intI)
            {
                if (validInterface[intI])
                {
                    word coupleType(fromProc);

                    newInterfaces.set
                    (
                        intI,
                        GAMGInterface::New
                        (
                            coupleType,
                            intI,
                            otherMeshes[proci-1].rawInterfaces(),
                            fromProc
                        ).ptr()
                    );
                }
            }

            otherMeshes[proci-1].addInterfaces
            (
                newInterfaces,
                nonBlockingSchedule<processorGAMGInterface>
                (
                    newInterfaces
                )
            );
        }
    }
}


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