fvMeshDistribute.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2018 OpenFOAM Foundation
    Copyright (C) 2015-2021 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 "fvMeshDistribute.H"
#include "PstreamCombineReduceOps.H"
#include "fvMeshAdder.H"
#include "faceCoupleInfo.H"
#include "processorFvPatchField.H"
#include "processorFvsPatchField.H"
#include "processorCyclicPolyPatch.H"
#include "processorCyclicFvPatchField.H"
#include "polyTopoChange.H"
#include "removeCells.H"
#include "polyModifyFace.H"
#include "polyRemovePoint.H"
#include "mapDistributePolyMesh.H"
#include "surfaceFields.H"
#include "syncTools.H"
#include "CompactListList.H"
#include "fvMeshTools.H"
#include "labelPairHashes.H"
#include "ListOps.H"
#include "globalIndex.H"
#include "cyclicACMIPolyPatch.H"
#include "mappedPatchBase.H"

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

namespace Foam
{
    defineTypeNameAndDebug(fvMeshDistribute, 0);

    //- Less function class that can be used for sorting processor patches
    class lessProcPatches
    {
        const labelList& nbrProc_;
        const labelList& referPatchID_;

    public:

        lessProcPatches(const labelList& nbrProc, const labelList& referPatchID)
        :
            nbrProc_(nbrProc),
            referPatchID_(referPatchID)
        {}

        bool operator()(const label a, const label b)
        {
            if (nbrProc_[a] < nbrProc_[b])
            {
                return true;
            }
            else if (nbrProc_[a] > nbrProc_[b])
            {
                return false;
            }
            else
            {
                // Equal neighbour processor
                return referPatchID_[a] < referPatchID_[b];
            }
        }
    };
}


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

void Foam::fvMeshDistribute::inplaceRenumberWithFlip
(
    const labelUList& oldToNew,
    const bool oldToNewHasFlip,
    const bool lstHasFlip,
    labelUList& lst
)
{
    if (!lstHasFlip && !oldToNewHasFlip)
    {
        Foam::inplaceRenumber(oldToNew, lst);
    }
    else
    {
        // Either input data or map encodes sign so result encodes sign

        forAll(lst, elemI)
        {
            // Extract old value and sign
            label val = lst[elemI];
            label sign = 1;
            if (lstHasFlip)
            {
                if (val > 0)
                {
                    val = val-1;
                }
                else if (val < 0)
                {
                    val = -val-1;
                    sign = -1;
                }
                else
                {
                    FatalErrorInFunction
                        << "Problem : zero value " << val
                        << " at index " << elemI << " out of " << lst.size()
                        << " list with flip bit" << exit(FatalError);
                }
            }


            // Lookup new value and possibly change sign
            label newVal = oldToNew[val];

            if (oldToNewHasFlip)
            {
                if (newVal > 0)
                {
                    newVal = newVal-1;
                }
                else if (newVal < 0)
                {
                    newVal = -newVal-1;
                    sign = -sign;
                }
                else
                {
                    FatalErrorInFunction
                        << "Problem : zero value " << newVal
                        << " at index " << elemI << " out of "
                        << oldToNew.size()
                        << " list with flip bit" << exit(FatalError);
                }
            }


            // Encode new value and sign
            lst[elemI] = sign*(newVal+1);
        }
    }
}


Foam::labelList Foam::fvMeshDistribute::select
(
    const bool selectEqual,
    const labelList& values,
    const label value
)
{
    label n = 0;

    forAll(values, i)
    {
        if (selectEqual == (values[i] == value))
        {
            n++;
        }
    }

    labelList indices(n);
    n = 0;

    forAll(values, i)
    {
        if (selectEqual == (values[i] == value))
        {
            indices[n++] = i;
        }
    }
    return indices;
}


Foam::wordList Foam::fvMeshDistribute::mergeWordList(const wordList& procNames)
{
    List<wordList> allNames(Pstream::nProcs());
    allNames[Pstream::myProcNo()] = procNames;
    Pstream::gatherList(allNames);

    // Assume there are few zone names to merge. Use HashSet otherwise (but
    // maintain ordering)
    DynamicList<word> mergedNames;
    if (Pstream::master())
    {
        mergedNames = procNames;
        for (const wordList& names : allNames)
        {
            for (const word& name : names)
            {
                mergedNames.appendUniq(name);
            }
        }
    }
    Pstream::scatter(mergedNames);

    return mergedNames;
}


void Foam::fvMeshDistribute::printMeshInfo(const fvMesh& mesh)
{
    Pout<< "Primitives:" << nl
        << "    points       :" << mesh.nPoints() << nl
        << "    bb           :" << boundBox(mesh.points(), false) << nl
        << "    internalFaces:" << mesh.nInternalFaces() << nl
        << "    faces        :" << mesh.nFaces() << nl
        << "    cells        :" << mesh.nCells() << nl;

    const fvBoundaryMesh& patches = mesh.boundary();

    Pout<< "Patches:" << endl;
    forAll(patches, patchi)
    {
        const polyPatch& pp = patches[patchi].patch();

        Pout<< "    " << patchi << " name:" << pp.name()
            << " size:" << pp.size()
            << " start:" << pp.start()
            << " type:" << pp.type()
            << endl;
    }

    if (mesh.pointZones().size())
    {
        Pout<< "PointZones:" << endl;
        forAll(mesh.pointZones(), zoneI)
        {
            const pointZone& pz = mesh.pointZones()[zoneI];
            Pout<< "    " << zoneI << " name:" << pz.name()
                << " size:" << pz.size()
                << endl;
        }
    }
    if (mesh.faceZones().size())
    {
        Pout<< "FaceZones:" << endl;
        forAll(mesh.faceZones(), zoneI)
        {
            const faceZone& fz = mesh.faceZones()[zoneI];
            Pout<< "    " << zoneI << " name:" << fz.name()
                << " size:" << fz.size()
                << endl;
        }
    }
    if (mesh.cellZones().size())
    {
        Pout<< "CellZones:" << endl;
        forAll(mesh.cellZones(), zoneI)
        {
            const cellZone& cz = mesh.cellZones()[zoneI];
            Pout<< "    " << zoneI << " name:" << cz.name()
                << " size:" << cz.size()
                << endl;
        }
    }
}


void Foam::fvMeshDistribute::printCoupleInfo
(
    const primitiveMesh& mesh,
    const labelList& sourceFace,
    const labelList& sourceProc,
    const labelList& sourcePatch,
    const labelList& sourceNewNbrProc
)
{
    Pout<< nl
        << "Current coupling info:"
        << endl;

    forAll(sourceFace, bFacei)
    {
        label meshFacei = mesh.nInternalFaces() + bFacei;

        Pout<< "    meshFace:" << meshFacei
            << " fc:" << mesh.faceCentres()[meshFacei]
            << " connects to proc:" << sourceProc[bFacei]
            << "/face:" << sourceFace[bFacei]
            << " which will move to proc:" << sourceNewNbrProc[bFacei]
            << endl;
    }
}


Foam::label Foam::fvMeshDistribute::findNonEmptyPatch() const
{
    // Finds (non-empty) patch that exposed internal and proc faces can be
    // put into.
    const polyBoundaryMesh& patches = mesh_.boundaryMesh();


    // Mark 'special' patches : -coupled, -duplicate faces. These probably
    // should not be used to (temporarily) store processor faces ...

    bitSet isCoupledPatch(patches.size());
    forAll(patches, patchi)
    {
        const polyPatch& pp = patches[patchi];
        const auto* cpp = isA<cyclicACMIPolyPatch>(pp);

        if (cpp)
        {
            isCoupledPatch.set(patchi);
            const label dupPatchID = cpp->nonOverlapPatchID();
            if (dupPatchID != -1)
            {
                isCoupledPatch.set(dupPatchID);
            }
        }
        else if (pp.coupled())
        {
            isCoupledPatch.set(patchi);
        }
    }

    label nonEmptyPatchi = -1;

    forAllReverse(patches, patchi)
    {
        const polyPatch& pp = patches[patchi];

        if
        (
           !isA<emptyPolyPatch>(pp)
        && !isCoupledPatch(patchi)
        && !isA<mappedPatchBase>(pp)
        )
        {
            nonEmptyPatchi = patchi;
            break;
        }
    }

    if (nonEmptyPatchi == -1)
    {
        FatalErrorInFunction
            << "Cannot find a patch which is neither of type empty nor"
            << " coupled in patches " << patches.names() << endl
            << "There has to be at least one such patch for"
            << " distribution to work" << abort(FatalError);
    }

    if (debug)
    {
        Pout<< "findNonEmptyPatch : using patch " << nonEmptyPatchi
            << " name:" << patches[nonEmptyPatchi].name()
            << " type:" << patches[nonEmptyPatchi].type()
            << " to put exposed faces into." << endl;
    }


    // Do additional test for processor patches intermingled with non-proc
    // patches.
    label procPatchi = -1;

    forAll(patches, patchi)
    {
        if (isA<processorPolyPatch>(patches[patchi]))
        {
            procPatchi = patchi;
        }
        else if (procPatchi != -1)
        {
            FatalErrorInFunction
                << "Processor patches should be at end of patch list."
                << endl
                << "Have processor patch " << procPatchi
                << " followed by non-processor patch " << patchi
                << " in patches " << patches.names()
                << abort(FatalError);
        }
    }

    return nonEmptyPatchi;
}


Foam::tmp<Foam::surfaceScalarField> Foam::fvMeshDistribute::generateTestField
(
    const fvMesh& mesh
)
{
    const vector testNormal = normalised(vector::one);

    tmp<surfaceScalarField> tfld
    (
        new surfaceScalarField
        (
            IOobject
            (
                "myFlux",
                mesh.time().timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            dimensionedScalar(dimless, Zero)
        )
    );
    surfaceScalarField& fld = tfld.ref();

    const surfaceVectorField n(mesh.Sf()/mesh.magSf());

    forAll(fld, facei)
    {
        fld[facei] = (n[facei] & testNormal);
    }

    surfaceScalarField::Boundary& fluxBf = fld.boundaryFieldRef();
    const surfaceVectorField::Boundary& nBf = n.boundaryField();

    forAll(fluxBf, patchi)
    {
        fvsPatchScalarField& fvp = fluxBf[patchi];

        scalarField newPfld(fvp.size());
        forAll(newPfld, i)
        {
            newPfld[i] = (nBf[patchi][i] & testNormal);
        }
        fvp == newPfld;
    }

    return tfld;
}


void Foam::fvMeshDistribute::testField(const surfaceScalarField& fld)
{
    const fvMesh& mesh = fld.mesh();

    const vector testNormal = normalised(vector::one);

    const surfaceVectorField n(mesh.Sf()/mesh.magSf());

    forAll(fld, facei)
    {
        scalar cos = (n[facei] & testNormal);

        if (mag(cos - fld[facei]) > 1e-6)
        {
            //FatalErrorInFunction
            WarningInFunction
                << "On internal face " << facei << " at "
                << mesh.faceCentres()[facei]
                << " the field value is " << fld[facei]
                << " whereas cos angle of " << testNormal
                << " with mesh normal " << n[facei]
                << " is " << cos
                //<< exit(FatalError);
                << endl;
        }
    }
    forAll(fld.boundaryField(), patchi)
    {
        const fvsPatchScalarField& fvp = fld.boundaryField()[patchi];
        const fvsPatchVectorField& np = n.boundaryField()[patchi];

        forAll(fvp, i)
        {
            scalar cos = (np[i] & testNormal);

            if (mag(cos - fvp[i]) > 1e-6)
            {
                label facei = fvp.patch().start()+i;
                //FatalErrorInFunction
                WarningInFunction
                    << "On face " << facei
                    << " on patch " << fvp.patch().name()
                    << " at " << mesh.faceCentres()[facei]
                    << " the field value is " << fvp[i]
                    << " whereas cos angle of " << testNormal
                    << " with mesh normal " << np[i]
                    << " is " << cos
                    //<< exit(FatalError);
                    << endl;
            }
        }
    }
}


Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::deleteProcPatches
(
    const label destinationPatch
)
{
    // Delete all processor patches. Move any processor faces into the last
    // non-processor patch.

    // New patchID per boundary faces to be repatched. Is -1 (no change)
    // or new patchID
    labelList newPatchID(mesh_.nBoundaryFaces(), -1);

    for (const polyPatch& pp : mesh_.boundaryMesh())
    {
        if (isA<processorPolyPatch>(pp))
        {
            if (debug)
            {
                Pout<< "Moving all faces of patch " << pp.name()
                    << " into patch " << destinationPatch
                    << endl;
            }

            SubList<label>
            (
                newPatchID,
                pp.size(),
                pp.offset()
            ) = destinationPatch;
        }
    }

    // Note: order of boundary faces been kept the same since the
    // destinationPatch is at the end and we have visited the patches in
    // incremental order.
    labelListList dummyFaceMaps;
    autoPtr<mapPolyMesh> map = repatch(newPatchID, dummyFaceMaps);


    // Delete (now empty) processor patches.
    {
        labelList oldToNew(identity(mesh_.boundaryMesh().size()));
        label newi = 0;
        // Non processor patches first
        forAll(mesh_.boundaryMesh(), patchi)
        {
            if (!isA<processorPolyPatch>(mesh_.boundaryMesh()[patchi]))
            {
                oldToNew[patchi] = newi++;
            }
        }
        label nNonProcPatches = newi;

        // Processor patches as last
        forAll(mesh_.boundaryMesh(), patchi)
        {
            if (isA<processorPolyPatch>(mesh_.boundaryMesh()[patchi]))
            {
                oldToNew[patchi] = newi++;
            }
        }
        fvMeshTools::reorderPatches(mesh_, oldToNew, nNonProcPatches, false);
    }

    return map;
}


Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::repatch
(
    const labelList& newPatchID,         // per boundary face -1 or new patchID
    labelListList& constructFaceMap
)
{
    polyTopoChange meshMod(mesh_);

    forAll(newPatchID, bFacei)
    {
        if (newPatchID[bFacei] != -1)
        {
            label facei = mesh_.nInternalFaces() + bFacei;

            label zoneID = mesh_.faceZones().whichZone(facei);
            bool zoneFlip = false;

            if (zoneID >= 0)
            {
                const faceZone& fZone = mesh_.faceZones()[zoneID];
                zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
            }

            meshMod.setAction
            (
                polyModifyFace
                (
                    mesh_.faces()[facei],       // modified face
                    facei,                      // label of face
                    mesh_.faceOwner()[facei],   // owner
                    -1,                         // neighbour
                    false,                      // face flip
                    newPatchID[bFacei],         // patch for face
                    false,                      // remove from zone
                    zoneID,                     // zone for face
                    zoneFlip                    // face flip in zone
                )
            );
        }
    }


    // Do mapping of fields from one patchField to the other ourselves since
    // is currently not supported by updateMesh.

    // Store boundary fields (we only do this for surfaceFields)
    PtrList<FieldField<fvsPatchField, scalar>> sFlds;
    saveBoundaryFields<scalar, surfaceMesh>(sFlds);
    PtrList<FieldField<fvsPatchField, vector>> vFlds;
    saveBoundaryFields<vector, surfaceMesh>(vFlds);
    PtrList<FieldField<fvsPatchField, sphericalTensor>> sptFlds;
    saveBoundaryFields<sphericalTensor, surfaceMesh>(sptFlds);
    PtrList<FieldField<fvsPatchField, symmTensor>> sytFlds;
    saveBoundaryFields<symmTensor, surfaceMesh>(sytFlds);
    PtrList<FieldField<fvsPatchField, tensor>> tFlds;
    saveBoundaryFields<tensor, surfaceMesh>(tFlds);

    // Change the mesh (no inflation). Note: parallel comms allowed.
    //
    // NOTE: there is one very particular problem with this ordering.
    // We first create the processor patches and use these to merge out
    // shared points (see mergeSharedPoints below). So temporarily points
    // and edges do not match!

    autoPtr<mapPolyMesh> mapPtr = meshMod.changeMesh(mesh_, false, true);
    mapPolyMesh& map = *mapPtr;

    // Update fields. No inflation, parallel sync.
    mesh_.updateMesh(map);

    // Map patch fields using stored boundary fields. Note: assumes order
    // of fields has not changed in object registry!
    mapBoundaryFields<scalar, surfaceMesh>(map, sFlds);
    mapBoundaryFields<vector, surfaceMesh>(map, vFlds);
    mapBoundaryFields<sphericalTensor, surfaceMesh>(map, sptFlds);
    mapBoundaryFields<symmTensor, surfaceMesh>(map, sytFlds);
    mapBoundaryFields<tensor, surfaceMesh>(map, tFlds);


    // Move mesh (since morphing does not do this)
    if (map.hasMotionPoints())
    {
        mesh_.movePoints(map.preMotionPoints());
    }

    // Adapt constructMaps.

    if (debug)
    {
        label index = map.reverseFaceMap().find(-1);

        if (index != -1)
        {
            FatalErrorInFunction
                << "reverseFaceMap contains -1 at index:"
                << index << endl
                << "This means that the repatch operation was not just"
                << " a shuffle?" << abort(FatalError);
        }
    }

    forAll(constructFaceMap, proci)
    {
        inplaceRenumberWithFlip
        (
            map.reverseFaceMap(),
            false,
            true,
            constructFaceMap[proci]
        );
    }

    return mapPtr;
}


// Detect shared points. Need processor patches to be present.
// Background: when adding bits of mesh one can get points which
// share the same position but are only detectable to be topologically
// the same point when doing parallel analysis. This routine will
// merge those points.
Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
(
    const labelList& pointToGlobalMaster,
    labelListList& constructPointMap
)
{
    // Find out which sets of points get merged and create a map from
    // mesh point to unique point.

    label nShared = 0;
    forAll(pointToGlobalMaster, pointi)
    {
        if (pointToGlobalMaster[pointi] != -1)
        {
            nShared++;
        }
    }

    if (debug)
    {
        Pout<< "mergeSharedPoints : found " << nShared
            << " points on processor boundaries" << nl << endl;
    }

    Map<label> globalMasterToLocalMaster(2*nShared);
    Map<label> pointToMaster(2*nShared);
    label nMatch = 0;

    forAll(pointToGlobalMaster, pointi)
    {
        label globali = pointToGlobalMaster[pointi];
        if (globali != -1)
        {
            const auto iter = globalMasterToLocalMaster.cfind(globali);

            if (iter.found())
            {
                // Matched to existing master
                pointToMaster.insert(pointi, *iter);
                nMatch++;
            }
            else
            {
                // Found first point. Designate as master
                globalMasterToLocalMaster.insert(globali, pointi);
                pointToMaster.insert(pointi, pointi);
            }
        }
    }

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

    if (debug)
    {
        Pout<< "mergeSharedPoints : found "
            << nMatch << " mergeable points" << nl << endl;
    }


    if (nMatch == 0)
    {
        return nullptr;
    }


    polyTopoChange meshMod(mesh_);

    fvMeshAdder::mergePoints(mesh_, pointToMaster, meshMod);

    // Change the mesh (no inflation). Note: parallel comms allowed.
    autoPtr<mapPolyMesh> mapPtr = meshMod.changeMesh(mesh_, false, true);
    mapPolyMesh& map = *mapPtr;

    // Update fields. No inflation, parallel sync.
    mesh_.updateMesh(map);

    // Adapt constructMaps for merged points.
    forAll(constructPointMap, proci)
    {
        labelList& constructMap = constructPointMap[proci];

        forAll(constructMap, i)
        {
            label oldPointi = constructMap[i];

            label newPointi = map.reversePointMap()[oldPointi];

            if (newPointi < -1)
            {
                constructMap[i] = -newPointi-2;
            }
            else if (newPointi >= 0)
            {
                constructMap[i] = newPointi;
            }
            else
            {
                FatalErrorInFunction
                    << "Problem. oldPointi:" << oldPointi
                    << " newPointi:" << newPointi << abort(FatalError);
            }
        }
    }

    return mapPtr;
}


void Foam::fvMeshDistribute::getCouplingData
(
    const labelList& distribution,
    labelList& sourceFace,
    labelList& sourceProc,
    labelList& sourcePatch,
    labelList& sourceNewNbrProc,
    labelList& sourcePointMaster
) const
{
    // Construct the coupling information for all (boundary) faces and
    // points

    const label nBnd = mesh_.nBoundaryFaces();
    sourceFace.setSize(nBnd);
    sourceProc.setSize(nBnd);
    sourcePatch.setSize(nBnd);
    sourceNewNbrProc.setSize(nBnd);

    const polyBoundaryMesh& patches = mesh_.boundaryMesh();

    // Get neighbouring meshFace labels and new processor of coupled boundaries.
    labelList nbrFaces(nBnd, -1);
    labelList nbrNewNbrProc(nBnd, -1);

    forAll(patches, patchi)
    {
        const polyPatch& pp = patches[patchi];

        if (pp.coupled())
        {
            label offset = pp.start() - mesh_.nInternalFaces();

            // Mesh labels of faces on this side
            forAll(pp, i)
            {
                label bndI = offset + i;
                nbrFaces[bndI] = pp.start()+i;
            }

            // Which processor they will end up on
            SubList<label>(nbrNewNbrProc, pp.size(), offset) =
                labelUIndList(distribution, pp.faceCells())();
        }
    }


    // Exchange the boundary data
    syncTools::swapBoundaryFaceList(mesh_, nbrFaces);
    syncTools::swapBoundaryFaceList(mesh_, nbrNewNbrProc);


    forAll(patches, patchi)
    {
        const polyPatch& pp = patches[patchi];
        label offset = pp.start() - mesh_.nInternalFaces();

        if (isA<processorPolyPatch>(pp))
        {
            const processorPolyPatch& procPatch =
                refCast<const processorPolyPatch>(pp);

            // Check which of the two faces we store.

            if (procPatch.owner())
            {
                // Use my local face labels
                forAll(pp, i)
                {
                    label bndI = offset + i;
                    sourceFace[bndI] = pp.start()+i;
                    sourceProc[bndI] = Pstream::myProcNo();
                    sourceNewNbrProc[bndI] = nbrNewNbrProc[bndI];
                }
            }
            else
            {
                // Use my neighbours face labels
                forAll(pp, i)
                {
                    label bndI = offset + i;
                    sourceFace[bndI] = nbrFaces[bndI];
                    sourceProc[bndI] = procPatch.neighbProcNo();
                    sourceNewNbrProc[bndI] = nbrNewNbrProc[bndI];
                }
            }


            label patchi = -1;
            if (isA<processorCyclicPolyPatch>(pp))
            {
                patchi = refCast<const processorCyclicPolyPatch>
                (
                    pp
                ).referPatchID();
            }

            forAll(pp, i)
            {
                label bndI = offset + i;
                sourcePatch[bndI] = patchi;
            }
        }
        else if (isA<cyclicPolyPatch>(pp))
        {
            const cyclicPolyPatch& cpp = refCast<const cyclicPolyPatch>(pp);

            if (cpp.owner())
            {
                forAll(pp, i)
                {
                    label bndI = offset + i;
                    sourceFace[bndI] = pp.start()+i;
                    sourceProc[bndI] = Pstream::myProcNo();
                    sourcePatch[bndI] = patchi;
                    sourceNewNbrProc[bndI] = nbrNewNbrProc[bndI];
                }
            }
            else
            {
                forAll(pp, i)
                {
                    label bndI = offset + i;
                    sourceFace[bndI] = nbrFaces[bndI];
                    sourceProc[bndI] = Pstream::myProcNo();
                    sourcePatch[bndI] = patchi;
                    sourceNewNbrProc[bndI] = nbrNewNbrProc[bndI];
                }
            }
        }
        else
        {
            // Normal (physical) boundary
            forAll(pp, i)
            {
                label bndI = offset + i;
                sourceFace[bndI] = -1;
                sourceProc[bndI] = -1;
                sourcePatch[bndI] = patchi;
                sourceNewNbrProc[bndI] = -1;
            }
        }
    }


    // Collect coupled (collocated) points
    sourcePointMaster.setSize(mesh_.nPoints());
    sourcePointMaster = -1;
    {
        // Assign global master point
        const globalIndex globalPoints(mesh_.nPoints());

        const globalMeshData& gmd = mesh_.globalData();
        const indirectPrimitivePatch& cpp = gmd.coupledPatch();
        const labelList& meshPoints = cpp.meshPoints();
        const mapDistribute& slavesMap = gmd.globalCoPointSlavesMap();
        const labelListList& slaves = gmd.globalCoPointSlaves();

        labelList elems(slavesMap.constructSize(), -1);
        forAll(meshPoints, pointi)
        {
            const labelList& slots = slaves[pointi];

            if (slots.size())
            {
                // pointi is a master. Assign a unique label.

                label globalPointi = globalPoints.toGlobal(meshPoints[pointi]);
                elems[pointi] = globalPointi;
                forAll(slots, i)
                {
                    label sloti = slots[i];
                    if (sloti >= meshPoints.size())
                    {
                        // Filter out local collocated points. We don't want
                        // to merge these
                        elems[slots[i]] = globalPointi;
                    }
                }
            }
        }

        // Push slave-slot data back to slaves
        slavesMap.reverseDistribute(elems.size(), elems, false);

        // Extract back onto mesh
        forAll(meshPoints, pointi)
        {
            sourcePointMaster[meshPoints[pointi]] = elems[pointi];
        }
    }
}