polyMeshAdder.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2019-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 "polyMeshAdder.H"
#include "mapAddedPolyMesh.H"
#include "IOobject.H"
#include "faceCoupleInfo.H"
#include "processorPolyPatch.H"
#include "SortableList.H"
#include "Time.H"
#include "globalMeshData.H"
#include "mergePoints.H"
#include "polyModifyFace.H"
#include "polyRemovePoint.H"
#include "polyTopoChange.H"
#include "globalIndex.H"

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

// Get index of patch in new set of patchnames/types
Foam::label Foam::polyMeshAdder::patchIndex
(
    const polyPatch& p,
    DynamicList<word>& allPatchNames,
    DynamicList<word>& allPatchTypes
)
{
    // Find the patch name on the list.  If the patch is already there
    // and patch types match, return index
    const word& pType = p.type();
    const word& pName = p.name();

    const label patchi = allPatchNames.find(pName);

    if (patchi == -1)
    {
        // Patch not found. Append to the list
        allPatchNames.append(pName);
        allPatchTypes.append(pType);

        return allPatchNames.size() - 1;
    }
    else if (allPatchTypes[patchi] == pType)
    {
        // Found name and types match
        return patchi;
    }
    else
    {
        // Found the name, but type is different

        // Duplicate name is not allowed.  Create a composite name from the
        // patch name and case name
        const word& caseName = p.boundaryMesh().mesh().time().caseName();

        allPatchNames.append(pName + "_" + caseName);
        allPatchTypes.append(pType);

        Pout<< "label patchIndex(const polyPatch& p) : "
            << "Patch " << p.index() << " named "
            << pName << " in mesh " << caseName
            << " already exists, but patch types"
            << " do not match.\nCreating a composite name as "
            << allPatchNames.last() << endl;

        return allPatchNames.size() - 1;
    }
}


// Get index of zone in new set of zone names
Foam::label Foam::polyMeshAdder::zoneIndex
(
    const word& curName,
    DynamicList<word>& names
)
{
    const label zoneI = names.find(curName);

    if (zoneI != -1)
    {
        return zoneI;
    }
    else
    {
        // Not found.  Add new name to the list
        names.append(curName);

        return names.size() - 1;
    }
}


void Foam::polyMeshAdder::mergePatchNames
(
    const polyBoundaryMesh& patches0,
    const polyBoundaryMesh& patches1,

    DynamicList<word>& allPatchNames,
    DynamicList<word>& allPatchTypes,

    labelList& from1ToAllPatches,
    labelList& fromAllTo1Patches
)
{
    // Insert the mesh0 patches and zones
    allPatchNames.append(patches0.names());
    allPatchTypes.append(patches0.types());


    // Patches
    // ~~~~~~~
    // Patches from 0 are taken over as is; those from 1 get either merged
    // (if they share name and type) or appended.
    // Empty patches are filtered out much much later on.

    // Add mesh1 patches and build map both ways.
    from1ToAllPatches.setSize(patches1.size());

    forAll(patches1, patchi)
    {
        from1ToAllPatches[patchi] = patchIndex
        (
            patches1[patchi],
            allPatchNames,
            allPatchTypes
        );
    }
    allPatchTypes.shrink();
    allPatchNames.shrink();

    // Invert 1 to all patch map
    fromAllTo1Patches.setSize(allPatchNames.size());
    fromAllTo1Patches = -1;

    forAll(from1ToAllPatches, i)
    {
        fromAllTo1Patches[from1ToAllPatches[i]] = i;
    }
}


Foam::labelList Foam::polyMeshAdder::getPatchStarts
(
    const polyBoundaryMesh& patches
)
{
    labelList patchStarts(patches.size());
    forAll(patches, patchi)
    {
        patchStarts[patchi] = patches[patchi].start();
    }
    return patchStarts;
}


Foam::labelList Foam::polyMeshAdder::getPatchSizes
(
    const polyBoundaryMesh& patches
)
{
    labelList patchSizes(patches.size());
    forAll(patches, patchi)
    {
        patchSizes[patchi] = patches[patchi].size();
    }
    return patchSizes;
}


Foam::List<Foam::polyPatch*> Foam::polyMeshAdder::combinePatches
(
    const polyMesh& mesh0,
    const polyMesh& mesh1,
    const polyBoundaryMesh& allBoundaryMesh,
    const label nAllPatches,
    const labelList& fromAllTo1Patches,

    const label nInternalFaces,
    const labelList& nFaces,

    labelList& from0ToAllPatches,
    labelList& from1ToAllPatches
)
{
    const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
    const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();


    // Compacted new patch list.
    DynamicList<polyPatch*> allPatches(nAllPatches);


    // Map from 0 to all patches (since gets compacted)
    from0ToAllPatches.setSize(patches0.size());
    from0ToAllPatches = -1;

    label startFacei = nInternalFaces;

    // Copy patches0 with new sizes. First patches always come from
    // mesh0 and will always be present.
    forAll(patches0, patchi)
    {
        // Originates from mesh0. Clone with new size & filter out empty
        // patch.
        label filteredPatchi;

        if (nFaces[patchi] == 0 && isA<processorPolyPatch>(patches0[patchi]))
        {
            //Pout<< "Removing zero sized mesh0 patch "
            //    << patches0[patchi].name() << endl;
            filteredPatchi = -1;
        }
        else
        {
            filteredPatchi = allPatches.size();

            allPatches.append
            (
                patches0[patchi].clone
                (
                    allBoundaryMesh,
                    filteredPatchi,
                    nFaces[patchi],
                    startFacei
                ).ptr()
            );
            startFacei += nFaces[patchi];
        }

        // Record new index in allPatches
        from0ToAllPatches[patchi] = filteredPatchi;

        // Check if patch was also in mesh1 and update its addressing if so.
        if (fromAllTo1Patches[patchi] != -1)
        {
            from1ToAllPatches[fromAllTo1Patches[patchi]] = filteredPatchi;
        }
    }

    // Copy unique patches of mesh1.
    forAll(from1ToAllPatches, patchi)
    {
        label allPatchi = from1ToAllPatches[patchi];

        if (allPatchi >= patches0.size())
        {
            // Patch has not been merged with any mesh0 patch.

            label filteredPatchi;

            if
            (
                nFaces[allPatchi] == 0
             && isA<processorPolyPatch>(patches1[patchi])
            )
            {
                //Pout<< "Removing zero sized mesh1 patch "
                //    << patches1[patchi].name() << endl;
                filteredPatchi = -1;
            }
            else
            {
                filteredPatchi = allPatches.size();

                allPatches.append
                (
                    patches1[patchi].clone
                    (
                        allBoundaryMesh,
                        filteredPatchi,
                        nFaces[allPatchi],
                        startFacei
                    ).ptr()
                );
                startFacei += nFaces[allPatchi];
            }

            from1ToAllPatches[patchi] = filteredPatchi;
        }
    }

    allPatches.shrink();

    return allPatches;
}


Foam::labelList Foam::polyMeshAdder::getFaceOrder
(
    const cellList& cells,
    const label nInternalFaces,
    const labelList& owner,
    const labelList& neighbour
)
{
    labelList oldToNew(owner.size(), -1);

    // Leave boundary faces in order
    for (label facei = nInternalFaces; facei < owner.size(); ++facei)
    {
        oldToNew[facei] = facei;
    }

    // First unassigned face
    label newFacei = 0;

    forAll(cells, celli)
    {
        const labelList& cFaces = cells[celli];

        SortableList<label> nbr(cFaces.size());

        forAll(cFaces, i)
        {
            const label facei = cFaces[i];

            label nbrCelli = neighbour[facei];

            if (nbrCelli != -1)
            {
                // Internal face. Get cell on other side.
                if (nbrCelli == celli)
                {
                    nbrCelli = owner[facei];
                }

                if (celli < nbrCelli)
                {
                    // Celli is master
                    nbr[i] = nbrCelli;
                }
                else
                {
                    // nbrCell is master. Let it handle this face.
                    nbr[i] = -1;
                }
            }
            else
            {
                // External face. Do later.
                nbr[i] = -1;
            }
        }

        nbr.sort();

        forAll(nbr, i)
        {
            if (nbr[i] != -1)
            {
                oldToNew[cFaces[nbr.indices()[i]]] = newFacei++;
            }
        }
    }


    // Check done all faces.
    forAll(oldToNew, facei)
    {
        if (oldToNew[facei] == -1)
        {
            FatalErrorInFunction
                << "Did not determine new position"
                << " for face " << facei
                << abort(FatalError);
        }
    }

    return oldToNew;
}


// Extends face f with split points. cutEdgeToPoints gives for every
// edge the points introduced inbetween the endpoints.
void Foam::polyMeshAdder::insertVertices
(
    const edgeLookup& cutEdgeToPoints,
    const Map<label>& meshToMaster,
    const labelList& masterToCutPoints,
    const face& masterF,

    DynamicList<label>& workFace,
    face& allF
)
{
    workFace.clear();

    // Check any edge for being cut (check on the cut so takes account
    // for any point merging on the cut)

    forAll(masterF, fp)
    {
        label v0 = masterF[fp];
        label v1 = masterF.nextLabel(fp);

        // Copy existing face point
        workFace.append(allF[fp]);

        // See if any edge between v0,v1

        const auto v0Fnd = meshToMaster.cfind(v0);
        if (v0Fnd.found())
        {
            const auto v1Fnd = meshToMaster.cfind(v1);
            if (v1Fnd.found())
            {
                // Get edge in cutPoint numbering
                edge cutEdge
                (
                    masterToCutPoints[v0Fnd()],
                    masterToCutPoints[v1Fnd()]
                );

                const auto iter = cutEdgeToPoints.cfind(cutEdge);

                if (iter.found())
                {
                    const edge& e = iter.key();
                    const labelList& addedPoints = iter.val();

                    // cutPoints first in allPoints so no need for renumbering
                    if (e[0] == cutEdge[0])
                    {
                        forAll(addedPoints, i)
                        {
                            workFace.append(addedPoints[i]);
                        }
                    }
                    else
                    {
                        forAllReverse(addedPoints, i)
                        {
                            workFace.append(addedPoints[i]);
                        }
                    }
                }
            }
        }
    }

    if (workFace.size() != allF.size())
    {
        allF.transfer(workFace);
    }
}


// Adds primitives (cells, faces, points)
// Cells:
//  - all of mesh0
//  - all of mesh1
// Faces:
//  - all uncoupled of mesh0
//  - all coupled faces
//  - all uncoupled of mesh1
// Points:
//  - all coupled
//  - all uncoupled of mesh0
//  - all uncoupled of mesh1
void Foam::polyMeshAdder::mergePrimitives
(
    const polyMesh& mesh0,
    const polyMesh& mesh1,
    const faceCoupleInfo& coupleInfo,

    const label nAllPatches,                // number of patches in the new mesh
    const labelList& fromAllTo1Patches,
    const labelList& from1ToAllPatches,

    pointField& allPoints,
    labelList& from0ToAllPoints,
    labelList& from1ToAllPoints,

    faceList& allFaces,
    labelList& allOwner,
    labelList& allNeighbour,
    label& nInternalFaces,
    labelList& nFacesPerPatch,
    label& nCells,

    labelList& from0ToAllFaces,
    labelList& from1ToAllFaces,
    labelList& from1ToAllCells
)
{
    const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
    const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();

    const primitiveFacePatch& cutFaces = coupleInfo.cutFaces();
    const indirectPrimitivePatch& masterPatch = coupleInfo.masterPatch();
    const indirectPrimitivePatch& slavePatch = coupleInfo.slavePatch();


    // Points
    // ~~~~~~

    // Storage for new points
    allPoints.setSize(mesh0.nPoints() + mesh1.nPoints());
    label allPointi = 0;

    from0ToAllPoints.setSize(mesh0.nPoints());
    from0ToAllPoints = -1;
    from1ToAllPoints.setSize(mesh1.nPoints());
    from1ToAllPoints = -1;

    // Copy coupled points (on cut)
    {
        const pointField& cutPoints = coupleInfo.cutPoints();

        //const labelListList& cutToMasterPoints =
        //   coupleInfo.cutToMasterPoints();
        labelListList cutToMasterPoints
        (
            invertOneToMany
            (
                cutPoints.size(),
                coupleInfo.masterToCutPoints()
            )
        );

        //const labelListList& cutToSlavePoints =
        //    coupleInfo.cutToSlavePoints();
        labelListList cutToSlavePoints
        (
            invertOneToMany
            (
                cutPoints.size(),
                coupleInfo.slaveToCutPoints()
            )
        );

        forAll(cutPoints, i)
        {
            allPoints[allPointi] = cutPoints[i];

            // Mark all master and slave points referring to this point.

            const labelList& masterPoints = cutToMasterPoints[i];

            forAll(masterPoints, j)
            {
                label mesh0Pointi = masterPatch.meshPoints()[masterPoints[j]];
                from0ToAllPoints[mesh0Pointi] = allPointi;
            }

            const labelList& slavePoints = cutToSlavePoints[i];

            forAll(slavePoints, j)
            {
                label mesh1Pointi = slavePatch.meshPoints()[slavePoints[j]];
                from1ToAllPoints[mesh1Pointi] = allPointi;
            }
            allPointi++;
        }
    }

    // Add uncoupled mesh0 points
    forAll(mesh0.points(), pointi)
    {
        if (from0ToAllPoints[pointi] == -1)
        {
            allPoints[allPointi] = mesh0.points()[pointi];
            from0ToAllPoints[pointi] = allPointi;
            allPointi++;
        }
    }

    // Add uncoupled mesh1 points
    forAll(mesh1.points(), pointi)
    {
        if (from1ToAllPoints[pointi] == -1)
        {
            allPoints[allPointi] = mesh1.points()[pointi];
            from1ToAllPoints[pointi] = allPointi;
            allPointi++;
        }
    }

    allPoints.setSize(allPointi);


    // Faces
    // ~~~~~

    // Sizes per patch
    nFacesPerPatch.setSize(nAllPatches);
    nFacesPerPatch = 0;

    // Storage for faces and owner/neighbour
    allFaces.setSize(mesh0.nFaces() + mesh1.nFaces());
    allOwner.setSize(allFaces.size());
    allOwner = -1;
    allNeighbour.setSize(allFaces.size());
    allNeighbour = -1;
    label allFacei = 0;

    from0ToAllFaces.setSize(mesh0.nFaces());
    from0ToAllFaces = -1;
    from1ToAllFaces.setSize(mesh1.nFaces());
    from1ToAllFaces = -1;

    // Copy mesh0 internal faces (always uncoupled)
    for (label facei = 0; facei < mesh0.nInternalFaces(); facei++)
    {
        allFaces[allFacei] = renumber(from0ToAllPoints, mesh0.faces()[facei]);
        allOwner[allFacei] = mesh0.faceOwner()[facei];
        allNeighbour[allFacei] = mesh0.faceNeighbour()[facei];
        from0ToAllFaces[facei] = allFacei++;
    }

    // Copy coupled faces. Every coupled face has an equivalent master and
    // slave. Also uncount as boundary faces all the newly coupled faces.
    const labelList& cutToMasterFaces = coupleInfo.cutToMasterFaces();
    const labelList& cutToSlaveFaces = coupleInfo.cutToSlaveFaces();

    forAll(cutFaces, i)
    {
        label masterFacei = cutToMasterFaces[i];

        label mesh0Facei = masterPatch.addressing()[masterFacei];

        if (from0ToAllFaces[mesh0Facei] == -1)
        {
            // First occurrence of face
            from0ToAllFaces[mesh0Facei] = allFacei;

            // External face becomes internal so uncount
            label patch0 = patches0.whichPatch(mesh0Facei);
            nFacesPerPatch[patch0]--;
        }

        label slaveFacei = cutToSlaveFaces[i];

        label mesh1Facei = slavePatch.addressing()[slaveFacei];

        if (from1ToAllFaces[mesh1Facei] == -1)
        {
            from1ToAllFaces[mesh1Facei] = allFacei;

            label patch1 = patches1.whichPatch(mesh1Facei);
            nFacesPerPatch[from1ToAllPatches[patch1]]--;
        }

        // Copy cut face (since cutPoints are copied first no renumbering
        // necessary)
        allFaces[allFacei] = cutFaces[i];
        allOwner[allFacei] = mesh0.faceOwner()[mesh0Facei];
        allNeighbour[allFacei] = mesh1.faceOwner()[mesh1Facei] + mesh0.nCells();

        allFacei++;
    }

    // Copy mesh1 internal faces (always uncoupled)
    for (label facei = 0; facei < mesh1.nInternalFaces(); facei++)
    {
        allFaces[allFacei] = renumber(from1ToAllPoints, mesh1.faces()[facei]);
        allOwner[allFacei] = mesh1.faceOwner()[facei] + mesh0.nCells();
        allNeighbour[allFacei] = mesh1.faceNeighbour()[facei] + mesh0.nCells();
        from1ToAllFaces[facei] = allFacei++;
    }

    nInternalFaces = allFacei;

    // Copy (unmarked/uncoupled) external faces in new order.
    for (label allPatchi = 0; allPatchi < nAllPatches; allPatchi++)
    {
        if (allPatchi < patches0.size())
        {
            // Patch is present in mesh0
            const polyPatch& pp = patches0[allPatchi];

            nFacesPerPatch[allPatchi] += pp.size();

            label facei = pp.start();

            forAll(pp, i)
            {
                if (from0ToAllFaces[facei] == -1)
                {
                    // Is uncoupled face since has not yet been dealt with
                    allFaces[allFacei] = renumber
                    (
                        from0ToAllPoints,
                        mesh0.faces()[facei]
                    );
                    allOwner[allFacei] = mesh0.faceOwner()[facei];
                    allNeighbour[allFacei] = -1;

                    from0ToAllFaces[facei] = allFacei++;
                }
                facei++;
            }
        }
        if (fromAllTo1Patches[allPatchi] != -1)
        {
            // Patch is present in mesh1
            const polyPatch& pp = patches1[fromAllTo1Patches[allPatchi]];

            nFacesPerPatch[allPatchi] += pp.size();

            label facei = pp.start();

            forAll(pp, i)
            {
                if (from1ToAllFaces[facei] == -1)
                {
                    // Is uncoupled face
                    allFaces[allFacei] = renumber
                    (
                        from1ToAllPoints,
                        mesh1.faces()[facei]
                    );
                    allOwner[allFacei] =
                        mesh1.faceOwner()[facei]
                      + mesh0.nCells();
                    allNeighbour[allFacei] = -1;

                    from1ToAllFaces[facei] = allFacei++;
                }
                facei++;
            }
        }
    }
    allFaces.setSize(allFacei);
    allOwner.setSize(allFacei);
    allNeighbour.setSize(allFacei);


    // So now we have all ok for one-to-one mapping.
    // For split slace faces:
    // - mesh consistent with slave side
    // - mesh not consistent with owner side. It is not zipped up, the
    //   original faces need edges split.

    // Use brute force to prevent having to calculate addressing:
    // - get map from master edge to split edges.
    // - check all faces to find any edge that is split.
    {
        // From two cut-points to labels of cut-points inbetween.
        // (in order: from e[0] to e[1]
        const edgeLookup& cutEdgeToPoints = coupleInfo.cutEdgeToPoints();

        // Get map of master face (in mesh labels) that are in cut. These faces
        // do not need to be renumbered.
        labelHashSet masterCutFaces(cutToMasterFaces.size());
        forAll(cutToMasterFaces, i)
        {
            label meshFacei = masterPatch.addressing()[cutToMasterFaces[i]];

            masterCutFaces.insert(meshFacei);
        }

        DynamicList<label> workFace(100);

        forAll(from0ToAllFaces, face0)
        {
            if (!masterCutFaces.found(face0))
            {
                label allFacei = from0ToAllFaces[face0];

                insertVertices
                (
                    cutEdgeToPoints,
                    masterPatch.meshPointMap(),
                    coupleInfo.masterToCutPoints(),
                    mesh0.faces()[face0],

                    workFace,
                    allFaces[allFacei]
                );
            }
        }

        // Same for slave face

        labelHashSet slaveCutFaces(cutToSlaveFaces.size());
        forAll(cutToSlaveFaces, i)
        {
            label meshFacei = slavePatch.addressing()[cutToSlaveFaces[i]];

            slaveCutFaces.insert(meshFacei);
        }

        forAll(from1ToAllFaces, face1)
        {
            if (!slaveCutFaces.found(face1))
            {
                label allFacei = from1ToAllFaces[face1];

                insertVertices
                (
                    cutEdgeToPoints,
                    slavePatch.meshPointMap(),
                    coupleInfo.slaveToCutPoints(),
                    mesh1.faces()[face1],

                    workFace,
                    allFaces[allFacei]
                );
            }
        }
    }

    // Now we have a full facelist and owner/neighbour addressing.


    // Cells
    // ~~~~~

    from1ToAllCells.setSize(mesh1.nCells());
    from1ToAllCells = -1;

    forAll(mesh1.cells(), i)
    {
        from1ToAllCells[i] = i + mesh0.nCells();
    }

    // Make cells (= cell-face addressing)
    nCells = mesh0.nCells() + mesh1.nCells();
    cellList allCells(nCells);
    primitiveMesh::calcCells(allCells, allOwner, allNeighbour, nCells);

    // Reorder faces for upper-triangular order.
    labelList oldToNew
    (
        getFaceOrder
        (
            allCells,
            nInternalFaces,
            allOwner,
            allNeighbour
        )
    );

    inplaceReorder(oldToNew, allFaces);
    inplaceReorder(oldToNew, allOwner);
    inplaceReorder(oldToNew, allNeighbour);
    inplaceRenumber(oldToNew, from0ToAllFaces);
    inplaceRenumber(oldToNew, from1ToAllFaces);
}


void Foam::polyMeshAdder::mergePointZones
(
    const label nAllPoints,
    const pointZoneMesh& pz0,
    const pointZoneMesh& pz1,
    const labelList& from0ToAllPoints,
    const labelList& from1ToAllPoints,

    DynamicList<word>& zoneNames,
    labelList& from1ToAll,
    List<DynamicList<label>>& pzPoints
)
{
    zoneNames.setCapacity(pz0.size() + pz1.size());
    zoneNames.append(pz0.names());

    from1ToAll.setSize(pz1.size());

    forAll(pz1, zoneI)
    {
        from1ToAll[zoneI] = zoneIndex(pz1[zoneI].name(), zoneNames);
    }
    zoneNames.shrink();



    // Zone(s) per point. Two levels: if only one zone
    // stored in pointToZone. Any extra stored in additionalPointToZones.
    // This is so we only allocate labelLists per point if absolutely
    // necessary.
    labelList pointToZone(nAllPoints, -1);
    labelListList addPointToZones(nAllPoints);

    // mesh0 zones kept
    forAll(pz0, zoneI)
    {
        const pointZone& pz = pz0[zoneI];

        forAll(pz, i)
        {
            label point0 = pz[i];
            label allPointi = from0ToAllPoints[point0];

            if (pointToZone[allPointi] == -1)
            {
                pointToZone[allPointi] = zoneI;
            }
            else if (pointToZone[allPointi] != zoneI)
            {
                labelList& pZones = addPointToZones[allPointi];
                pZones.appendUniq(zoneI);
            }
        }
    }

    // mesh1 zones renumbered
    forAll(pz1, zoneI)
    {
        const pointZone& pz = pz1[zoneI];
        const label allZoneI = from1ToAll[zoneI];

        forAll(pz, i)
        {
            label point1 = pz[i];
            label allPointi = from1ToAllPoints[point1];

            if (pointToZone[allPointi] == -1)
            {
                pointToZone[allPointi] = allZoneI;
            }
            else if (pointToZone[allPointi] != allZoneI)
            {
                labelList& pZones = addPointToZones[allPointi];
                pZones.appendUniq(allZoneI);
            }
        }
    }

    // Extract back into zones

    // 1. Count
    labelList nPoints(zoneNames.size(), Zero);
    forAll(pointToZone, allPointi)
    {
        label zoneI = pointToZone[allPointi];
        if (zoneI != -1)
        {
            nPoints[zoneI]++;
        }
    }
    forAll(addPointToZones, allPointi)
    {
        const labelList& pZones = addPointToZones[allPointi];
        forAll(pZones, i)
        {
            nPoints[pZones[i]]++;
        }
    }

    // 2. Fill
    pzPoints.setSize(zoneNames.size());
    forAll(pzPoints, zoneI)
    {
        pzPoints[zoneI].setCapacity(nPoints[zoneI]);
    }
    forAll(pointToZone, allPointi)
    {
        label zoneI = pointToZone[allPointi];
        if (zoneI != -1)
        {
            pzPoints[zoneI].append(allPointi);
        }
    }
    forAll(addPointToZones, allPointi)
    {
        const labelList& pZones = addPointToZones[allPointi];
        forAll(pZones, i)
        {
            pzPoints[pZones[i]].append(allPointi);
        }
    }
    forAll(pzPoints, i)
    {
        pzPoints[i].shrink();
        stableSort(pzPoints[i]);
    }
}


void Foam::polyMeshAdder::mergeFaceZones
(
    const labelList& allOwner,