openfoam/src/dynamicMesh/fvMeshSubset/fvMeshSubset.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
    Copyright (C) 2015-2018 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 "fvMeshSubset.H"
#include "boolList.H"
#include "BitOps.H"
#include "pointIndList.H"
#include "Pstream.H"
#include "emptyPolyPatch.H"
#include "demandDrivenData.H"
#include "cyclicPolyPatch.H"
#include "removeCells.H"
#include "polyTopoChange.H"
#include "mapPolyMesh.H"

// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //

namespace
{

// Mark faces/points (with 0) in labelList
inline void markUsed
(
    const Foam::labelUList& locations,
    Foam::labelList& map
)
{
    for (auto idx : locations)
    {
        map[idx] = 0;
    }
}

} // End anonymous namespace


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

Foam::word Foam::fvMeshSubset::exposedPatchName("oldInternalFaces");


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

bool Foam::fvMeshSubset::checkCellSubset() const
{
    if (!fvMeshSubsetPtr_)
    {
        FatalErrorInFunction
            << "setCellSubset()" << nl
            << "before attempting to access subset data"
            << abort(FatalError);

        return false;
    }

    return true;
}


void Foam::fvMeshSubset::calcFaceFlipMap() const
{
    const labelList& subToBaseFace = faceMap();
    const labelList& subToBaseCell = cellMap();

    faceFlipMapPtr_.reset(new labelList(subToBaseFace.size()));
    auto& faceFlipMap = *faceFlipMapPtr_;

    // Only exposed internal faces might be flipped (since we don't do
    // any cell renumbering, just compacting)
    const label subInt = subMesh().nInternalFaces();

    const labelList& subOwn = subMesh().faceOwner();
    const labelList& own = baseMesh_.faceOwner();

    for (label subFaceI = 0; subFaceI < subInt; ++subFaceI)
    {
        faceFlipMap[subFaceI] = subToBaseFace[subFaceI]+1;
    }
    for (label subFaceI = subInt; subFaceI < subOwn.size(); ++subFaceI)
    {
        const label faceI = subToBaseFace[subFaceI];
        if (subToBaseCell[subOwn[subFaceI]] == own[faceI])
        {
            faceFlipMap[subFaceI] = faceI+1;
        }
        else
        {
            faceFlipMap[subFaceI] = -faceI-1;
        }
    }
}


void Foam::fvMeshSubset::doCoupledPatches
(
    const bool syncPar,
    labelList& nCellsUsingFace
) const
{
    // Synchronize facesToSubset on both sides of coupled patches.
    // Marks faces that become 'uncoupled' with 3.

    const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();

    label nUncoupled = 0;

    if (syncPar && Pstream::parRun())
    {
        PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);

        // Send face usage across processor patches
        for (const polyPatch& pp : oldPatches)
        {
            if (isA<processorPolyPatch>(pp))
            {
                const processorPolyPatch& procPatch =
                    refCast<const processorPolyPatch>(pp);

                UOPstream toNeighbour(procPatch.neighbProcNo(), pBufs);

                if (!nCellsUsingFace.empty())
                {
                    toNeighbour <<
                        SubList<label>(nCellsUsingFace, pp.size(), pp.start());
                }
                else
                {
                    toNeighbour << labelList();
                }
            }
        }

        pBufs.finishedSends();

        // Receive face usage count and check for faces that become uncoupled.
        for (const polyPatch& pp : oldPatches)
        {
            if (isA<processorPolyPatch>(pp))
            {
                const processorPolyPatch& procPatch =
                    refCast<const processorPolyPatch>(pp);

                UIPstream fromNeighbour(procPatch.neighbProcNo(), pBufs);

                const labelList nbrList(fromNeighbour);

                // Combine with this side.

                if (!nCellsUsingFace.empty())
                {
                    const labelList& nbrCellsUsingFace(nbrList);

                    // Combine with this side.

                    forAll(pp, i)
                    {
                        if
                        (
                            nCellsUsingFace[pp.start()+i] == 1
                         && nbrCellsUsingFace[i] == 0
                        )
                        {
                            // Face's neighbour is no longer there. Mark face
                            // off as coupled
                            nCellsUsingFace[pp.start()+i] = 3;
                            ++nUncoupled;
                        }
                    }
                }
            }
        }
    }

    // Do same for cyclics.
    for (const polyPatch& pp : oldPatches)
    {
        if (isA<cyclicPolyPatch>(pp))
        {
            const cyclicPolyPatch& cycPatch =
                refCast<const cyclicPolyPatch>(pp);

            if (!nCellsUsingFace.empty())
            {
                forAll(cycPatch, i)
                {
                    label thisFacei = cycPatch.start() + i;
                    label otherFacei = cycPatch.transformGlobalFace(thisFacei);

                    if
                    (
                        nCellsUsingFace[thisFacei] == 1
                     && nCellsUsingFace[otherFacei] == 0
                    )
                    {
                        nCellsUsingFace[thisFacei] = 3;
                        ++nUncoupled;
                    }
                }
            }
        }
    }

    if (syncPar)
    {
        reduce(nUncoupled, sumOp<label>());
    }

    if (nUncoupled > 0)
    {
        Info<< "Uncoupled " << nUncoupled << " faces on coupled patches. "
            << "(processorPolyPatch, cyclicPolyPatch)" << endl;
    }
}


void Foam::fvMeshSubset::removeCellsImpl
(
    const bitSet& cellsToRemove,
    const labelList& exposedFaces,
    const labelList& patchIDs,
    const bool syncCouples
)
{
    // Mesh changing engine.
    polyTopoChange meshMod(baseMesh());

    removeCells cellRemover(baseMesh(), syncCouples);

    cellRemover.setRefinement
    (
        cellsToRemove,
        exposedFaces,
        patchIDs,
        meshMod
    );

    // Create mesh, return map from old to new mesh.
    autoPtr<mapPolyMesh> map = meshMod.makeMesh
    (
        fvMeshSubsetPtr_,
        IOobject
        (
            baseMesh().name(),
            baseMesh().time().timeName(),
            baseMesh().time(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        baseMesh(),
        syncCouples
    );

    pointMap_ = map().pointMap();
    faceMap_ = map().faceMap();
    cellMap_ = map().cellMap();
    patchMap_ = identity(baseMesh().boundaryMesh().size());
}


Foam::labelList Foam::fvMeshSubset::subsetSubset
(
    const label nElems,
    const labelUList& selectedElements,
    const labelUList& subsetMap
)
{
    // Mark selected elements.
    const bitSet selected(nElems, selectedElements);

    // Count subset of selected elements
    label n = 0;
    forAll(subsetMap, i)
    {
        if (selected[subsetMap[i]])
        {
            ++n;
        }
    }

    // Collect selected elements
    labelList subsettedElements(n);
    n = 0;

    forAll(subsetMap, i)
    {
        if (selected[subsetMap[i]])
        {
            subsettedElements[n] = i;
            ++n;
        }
    }

    return subsettedElements;
}


void Foam::fvMeshSubset::subsetZones()
{
    // Keep all zones, even if zero size.

    // PointZones

    const pointZoneMesh& pointZones = baseMesh().pointZones();

    List<pointZone*> pZonePtrs(pointZones.size());

    forAll(pointZones, zonei)
    {
        const pointZone& pz = pointZones[zonei];

        pZonePtrs[zonei] = new pointZone
        (
            pz.name(),
            subsetSubset(baseMesh().nPoints(), pz, pointMap()),
            zonei,
            fvMeshSubsetPtr_().pointZones()
        );
    }


    // FaceZones
    // Do we need to remove zones where the side we're interested in
    // no longer exists? Guess not.

    const faceZoneMesh& faceZones = baseMesh().faceZones();

    List<faceZone*> fZonePtrs(faceZones.size());

    forAll(faceZones, zonei)
    {
        const faceZone& fz = faceZones[zonei];

        // Expand faceZone to full mesh
        // +1 : part of faceZone, flipped
        // -1 :    ,,           , unflipped
        //  0 : not part of faceZone
        labelList zone(baseMesh().nFaces(), Zero);
        forAll(fz, j)
        {
            if (fz.flipMap()[j])
            {
                zone[fz[j]] = 1;
            }
            else
            {
                zone[fz[j]] = -1;
            }
        }

        // Select faces
        label nSub = 0;
        forAll(faceMap(), j)
        {
            if (zone[faceMap()[j]] != 0)
            {
                ++nSub;
            }
        }
        labelList subAddressing(nSub);
        boolList subFlipStatus(nSub);
        nSub = 0;
        forAll(faceMap(), subFacei)
        {
            const label meshFacei = faceMap()[subFacei];
            if (zone[meshFacei] != 0)
            {
                subAddressing[nSub] = subFacei;
                const label subOwner = subMesh().faceOwner()[subFacei];
                const label baseOwner = baseMesh().faceOwner()[meshFacei];
                // If subowner is the same cell as the base keep the flip status
                const bool sameOwner = (cellMap()[subOwner] == baseOwner);
                const bool flip = (zone[meshFacei] == 1);
                subFlipStatus[nSub] = (sameOwner == flip);

                ++nSub;
            }
        }

        fZonePtrs[zonei] = new faceZone
        (
            fz.name(),
            subAddressing,
            subFlipStatus,
            zonei,
            fvMeshSubsetPtr_().faceZones()
        );
    }

    // Cell Zones

    const cellZoneMesh& cellZones = baseMesh().cellZones();

    List<cellZone*> cZonePtrs(cellZones.size());

    forAll(cellZones, zonei)
    {
        const cellZone& cz = cellZones[zonei];

        cZonePtrs[zonei] = new cellZone
        (
            cz.name(),
            subsetSubset(baseMesh().nCells(), cz, cellMap()),
            zonei,
            fvMeshSubsetPtr_().cellZones()
        );
    }


    // Add the zones
    fvMeshSubsetPtr_().addZones(pZonePtrs, fZonePtrs, cZonePtrs);
}


Foam::bitSet Foam::fvMeshSubset::getCellsToRemove
(
    const bitSet& selectedCells
) const
{
    // Work on a copy
    bitSet cellsToRemove(selectedCells);

    // Ensure we have the full range
    cellsToRemove.resize(baseMesh().nCells(), false);

    // Invert the selection
    cellsToRemove.flip();

    return cellsToRemove;
}


Foam::bitSet Foam::fvMeshSubset::getCellsToRemove
(
    const label regioni,
    const labelUList& regions
) const
{
    return BitSetOps::create
    (
        baseMesh().nCells(),
        regioni,
        regions,
        false  // on=false: invert return cells to remove
    );
}


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

Foam::fvMeshSubset::fvMeshSubset(const fvMesh& baseMesh)
:
    baseMesh_(baseMesh),
    fvMeshSubsetPtr_(nullptr),
    faceFlipMapPtr_(nullptr),
    pointMap_(),
    faceMap_(),
    cellMap_(),
    patchMap_()
{}


Foam::fvMeshSubset::fvMeshSubset
(
    const fvMesh& baseMesh,
    const bitSet& selectedCells,
    const label patchID,
    const bool syncCouples
)
:
    fvMeshSubset(baseMesh)
{
    setCellSubset(selectedCells, patchID, syncCouples);
}


Foam::fvMeshSubset::fvMeshSubset
(
    const fvMesh& baseMesh,
    const labelHashSet& selectedCells,
    const label patchID,
    const bool syncCouples
)
:
    fvMeshSubset(baseMesh)
{
    setCellSubset(selectedCells, patchID, syncCouples);
}


Foam::fvMeshSubset::fvMeshSubset
(
    const fvMesh& baseMesh,
    const labelUList& selectedCells,
    const label patchID,
    const bool syncCouples
)
:
    fvMeshSubset(baseMesh)
{
    setCellSubset(selectedCells, patchID, syncCouples);
}


Foam::fvMeshSubset::fvMeshSubset
(
    const fvMesh& baseMesh,
    const label regioni,
    const labelUList& regions,
    const label patchID,
    const bool syncCouples
)
:
    fvMeshSubset(baseMesh)
{
    setCellSubset(regioni, regions, patchID, syncCouples);
}


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

void Foam::fvMeshSubset::clear()
{
    fvMeshSubsetPtr_.reset(nullptr);
    faceFlipMapPtr_.reset(nullptr);

    pointMap_.clear();
    faceMap_.clear();
    cellMap_.clear();
    patchMap_.clear();
}


void Foam::fvMeshSubset::setCellSubset
(
    const bitSet& selectedCells,
    const label patchID,
    const bool syncPar
)
{
    const cellList& oldCells = baseMesh().cells();
    const faceList& oldFaces = baseMesh().faces();
    const pointField& oldPoints = baseMesh().points();
    const labelList& oldOwner = baseMesh().faceOwner();
    const labelList& oldNeighbour = baseMesh().faceNeighbour();
    const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
    const label oldNInternalFaces = baseMesh().nInternalFaces();

    // Initial check on patches before doing anything time consuming.

    label wantedPatchID = patchID;

    if (wantedPatchID == -1)
    {
        // No explicit patch specified. Put in oldInternalFaces patch.
        // Check if patch with this name already exists.
        wantedPatchID = oldPatches.findPatchID(exposedPatchName);
    }
    else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.size())
    {
        FatalErrorInFunction
            << "Non-existing patch index " << wantedPatchID << endl
            << "Should be between 0 and " << oldPatches.size()-1
            << abort(FatalError);
    }

    // Clear all old maps and pointers
    clear();

    // The selected cells - sorted in ascending order
    cellMap_ = selectedCells.sortedToc();

    // The selectedCells should normally be in the [0,nCells) range,
    // but it is better not to trust that.
    {
        label len = cellMap_.size();
        for
        (
            label i = len-1;
            i >= 0 && (cellMap_[i] >= oldCells.size());
            --i
        )
        {
            len = i;
        }
        cellMap_.resize(len);
    }


    // Mark all used faces. Count number of cells using them
    // 0: face not used anymore
    // 1: face used by one cell, face becomes/stays boundary face
    // 2: face still used and remains internal face
    // 3: face coupled and used by one cell only (so should become normal,
    //    non-coupled patch face)
    //
    // Note that this is not really necessary - but means we can size things
    // correctly. Also makes handling coupled faces much easier.

    labelList nCellsUsingFace(oldFaces.size(), Zero);

    label nFacesInSet = 0;
    forAll(oldFaces, oldFacei)
    {
        bool faceUsed = false;

        if (selectedCells.test(oldOwner[oldFacei]))
        {
            ++nCellsUsingFace[oldFacei];
            faceUsed = true;
        }

        if
        (
            baseMesh().isInternalFace(oldFacei)
         && selectedCells.test(oldNeighbour[oldFacei])
        )
        {
            ++nCellsUsingFace[oldFacei];
            faceUsed = true;
        }

        if (faceUsed)
        {
            ++nFacesInSet;
        }
    }
    faceMap_.setSize(nFacesInSet);

    // Handle coupled faces. Modifies patch faces to be uncoupled to 3.
    doCoupledPatches(syncPar, nCellsUsingFace);


    // See which patch to use for exposed internal faces.
    label oldInternalPatchID = 0;

    // Insert faces before which patch
    label nextPatchID = oldPatches.size();

    // old to new patches
    labelList globalPatchMap(oldPatches.size());

    // New patch size
    label nbSize = oldPatches.size();

    if (wantedPatchID == -1)
    {
        // Create 'oldInternalFaces' patch at the end (or before
        // processorPatches)
        // and put all exposed internal faces in there.

        forAll(oldPatches, patchi)
        {
            if (isA<processorPolyPatch>(oldPatches[patchi]))
            {
                nextPatchID = patchi;
                break;
            }
            ++oldInternalPatchID;
        }

        ++nbSize;

        // adapt old to new patches for inserted patch
        for (label oldPatchi = 0; oldPatchi < nextPatchID; oldPatchi++)
        {
            globalPatchMap[oldPatchi] = oldPatchi;
        }
        for
        (
            label oldPatchi = nextPatchID;
            oldPatchi < oldPatches.size();
            oldPatchi++
        )
        {
            globalPatchMap[oldPatchi] = oldPatchi+1;
        }
    }
    else
    {
        oldInternalPatchID = wantedPatchID;
        nextPatchID = wantedPatchID+1;

        // old to new patches
        globalPatchMap = identity(oldPatches.size());
    }

    labelList boundaryPatchSizes(nbSize, Zero);


    // Make a global-to-local point map
    labelList globalPointMap(oldPoints.size(), -1);
    labelList globalFaceMap(oldFaces.size(), -1);

    label facei = 0;

    // 1. Pick up all preserved internal faces.
    for (label oldFacei = 0; oldFacei < oldNInternalFaces; ++oldFacei)
    {
        if (nCellsUsingFace[oldFacei] == 2)
        {
            globalFaceMap[oldFacei] = facei;
            faceMap_[facei++] = oldFacei;

            // Mark all points from the face
            markUsed(oldFaces[oldFacei], globalPointMap);
        }
    }

    // These are all the internal faces in the mesh.
    const label nInternalFaces = facei;

    // 2. Boundary faces up to where we want to insert old internal faces
    for
    (
        label oldPatchi = 0;
        oldPatchi < oldPatches.size()
     && oldPatchi < nextPatchID;
        oldPatchi++
    )
    {
        const polyPatch& oldPatch = oldPatches[oldPatchi];

        label oldFacei = oldPatch.start();

        forAll(oldPatch, i)
        {
            if (nCellsUsingFace[oldFacei] == 1)
            {
                // Boundary face is kept.

                // Mark face and increment number of points in set
                globalFaceMap[oldFacei] = facei;
                faceMap_[facei++] = oldFacei;

                // Mark all points from the face
                markUsed(oldFaces[oldFacei], globalPointMap);

                // Increment number of patch faces
                ++boundaryPatchSizes[globalPatchMap[oldPatchi]];
            }
            ++oldFacei;
        }
    }

    // 3a. old internal faces that have become exposed.
    for (label oldFacei = 0; oldFacei < oldNInternalFaces; ++oldFacei)
    {
        if (nCellsUsingFace[oldFacei] == 1)
        {
            globalFaceMap[oldFacei] = facei;
            faceMap_[facei++] = oldFacei;

            // Mark all points from the face
            markUsed(oldFaces[oldFacei], globalPointMap);

            // Increment number of patch faces
            ++boundaryPatchSizes[oldInternalPatchID];
        }
    }

    // 3b. coupled patch faces that have become uncoupled.
    for
    (
        label oldFacei = oldNInternalFaces;
        oldFacei < oldFaces.size();
        oldFacei++
    )
    {
        if (nCellsUsingFace[oldFacei] == 3)
        {
            globalFaceMap[oldFacei] = facei;
            faceMap_[facei++] = oldFacei;

            // Mark all points from the face
            markUsed(oldFaces[oldFacei], globalPointMap);

            // Increment number of patch faces
            ++boundaryPatchSizes[oldInternalPatchID];
        }
    }

    // 4. Remaining boundary faces
    for
    (
        label oldPatchi = nextPatchID;
        oldPatchi < oldPatches.size();
        oldPatchi++
    )
    {
        const polyPatch& oldPatch = oldPatches[oldPatchi];

        label oldFacei = oldPatch.start();

        forAll(oldPatch, i)
        {
            if (nCellsUsingFace[oldFacei] == 1)
            {
                // Boundary face is kept.

                // Mark face and increment number of points in set
                globalFaceMap[oldFacei] = facei;
                faceMap_[facei++] = oldFacei;

                // Mark all points from the face
                markUsed(oldFaces[oldFacei], globalPointMap);

                // Increment number of patch faces
                ++boundaryPatchSizes[globalPatchMap[oldPatchi]];
            }
            ++oldFacei;
        }
    }

    if (facei != nFacesInSet)
    {
        FatalErrorInFunction
            << "Problem" << abort(FatalError);
    }


    // Grab the points map
    label nPointsInSet = 0;

    forAll(globalPointMap, pointi)
    {
        if (globalPointMap[pointi] != -1)
        {
            ++nPointsInSet;
        }
    }
    pointMap_.setSize(nPointsInSet);

    nPointsInSet = 0;

    forAll(globalPointMap, pointi)
    {
        if (globalPointMap[pointi] != -1)
        {
            pointMap_[nPointsInSet] = pointi;
            globalPointMap[pointi] = nPointsInSet;
            ++nPointsInSet;
        }
    }


    //Pout<< "Number of points,faces,cells in new mesh : "
    //    << pointMap_.size() << ' '
    //    << faceMap_.size() << ' '
    //    << cellMap_.size() << nl;


    // Make a new mesh

    //
    // Create new points
    //
    pointField newPoints(pointUIndList(oldPoints, pointMap_));


    //
    // Create new faces
    //

    faceList newFaces(faceMap_.size());
    {
        auto iter = newFaces.begin();
        const auto& renumbering = globalPointMap;

        // Internal faces
        for (label facei = 0; facei < nInternalFaces; ++facei)
        {
            face& newItem = *iter;
            ++iter;

            const face& oldItem = oldFaces[faceMap_[facei]];

            // Copy relabelled
            newItem.resize(oldItem.size());

            forAll(oldItem, i)
            {
                newItem[i] = renumbering[oldItem[i]];
            }
        }

        // Boundary faces - may need to be flipped
        for (label facei = nInternalFaces; facei < faceMap_.size(); ++facei)
        {
            const label oldFacei = faceMap_[facei];

            face& newItem = *iter;
            ++iter;

            const face oldItem =
            (
                (
                    baseMesh().isInternalFace(oldFacei)
                 && selectedCells.test(oldNeighbour[oldFacei])
                 && !selectedCells.test(oldOwner[oldFacei])
                )
                // Face flipped to point outwards
              ? oldFaces[oldFacei].reverseFace()
              : oldFaces[oldFacei]
            );

            // Copy relabelled
            newItem.resize(oldItem.size());

            forAll(oldItem, i)
            {
                newItem[i] = renumbering[oldItem[i]];
            }
        }
    }


    //
    // Create new cells
    //

    cellList newCells(cellMap_.size());
    {
        auto iter = newCells.begin();
        const auto& renumbering = globalFaceMap;

        for (const label oldCelli : cellMap_)
        {
            cell& newItem = *iter;
            ++iter;

            const labelList& oldItem = oldCells[oldCelli];

            // Copy relabelled
            newItem.resize(oldItem.size());

            forAll(oldItem, i)
            {
                newItem[i] = renumbering[oldItem[i]];
            }
        }
    }


    // Make a new mesh
    //
    // Note that mesh gets registered with same name as original mesh.
    // This is not proper but cannot be avoided since otherwise
    // surfaceInterpolation cannot find its fvSchemes.
    // It will try to read, for example.  "system/region0SubSet/fvSchemes"
    //
    fvMeshSubsetPtr_ = autoPtr<fvMesh>::New
    (
        IOobject
        (
            baseMesh().name(),
            baseMesh().time().timeName(),
            baseMesh().time(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        std::move(newPoints),
        std::move(newFaces),
        std::move(newCells),
        syncPar           // parallel synchronisation
    );


    // Add old patches
    List<polyPatch*> newBoundary(nbSize);
    patchMap_.setSize(nbSize);
    label nNewPatches = 0;
    label patchStart = nInternalFaces;


    // For parallel: only remove patch if none of the processors has it.
    // This only gets done for patches before the one being inserted
    // (so patches < nextPatchID)

    // Get sum of patch sizes. Zero if patch can be deleted.
    labelList globalPatchSizes(boundaryPatchSizes);
    globalPatchSizes.setSize(nextPatchID);

    if (syncPar && Pstream::parRun())
    {
        // Get patch names (up to nextPatchID)
        List<wordList> patchNames(Pstream::nProcs());
        patchNames[Pstream::myProcNo()] = oldPatches.names();
        patchNames[Pstream::myProcNo()].setSize(nextPatchID);
        Pstream::gatherList(patchNames);
        Pstream::scatterList(patchNames);

        // Get patch sizes (up to nextPatchID).
        // Note that up to nextPatchID the globalPatchMap is an identity so
        // no need to index through that.
        Pstream::listCombineGather(globalPatchSizes, plusEqOp<label>());
        Pstream::listCombineScatter(globalPatchSizes);

        // Now all processors have all the patchnames.
        // Decide: if all processors have the same patch names and size is zero
        // everywhere remove the patch.
        bool samePatches = true;

        for (label proci = 1; proci < patchNames.size(); ++proci)
        {
            if (patchNames[proci] != patchNames[0])
            {
                samePatches = false;
                break;
            }
        }

        if (!samePatches)
        {
            // Patchnames not sync on all processors so disable removal of
            // zero sized patches.
            globalPatchSizes = labelMax;
        }
    }


    // Old patches

    for
    (
        label oldPatchi = 0;
        oldPatchi < oldPatches.size()
     && oldPatchi < nextPatchID;
        oldPatchi++
    )
    {
        const label newSize = boundaryPatchSizes[globalPatchMap[oldPatchi]];

        // Clone (even if 0 size)
        newBoundary[nNewPatches] = oldPatches[oldPatchi].clone
        (
            fvMeshSubsetPtr_().boundaryMesh(),
            nNewPatches,
            newSize,
            patchStart
        ).ptr();

        patchStart += newSize;
        patchMap_[nNewPatches] = oldPatchi;    // compact patchMap
        ++nNewPatches;
    }

    // Inserted patch

    if (wantedPatchID == -1)
    {
        label oldInternalSize = boundaryPatchSizes[oldInternalPatchID];

        if (syncPar)
        {
            reduce(oldInternalSize, sumOp<label>());
        }

        // Newly created patch so is at end. Check if any faces in it.
        if (oldInternalSize > 0)
        {
            newBoundary[nNewPatches] = new emptyPolyPatch
            (
                exposedPatchName,
                boundaryPatchSizes[oldInternalPatchID],
                patchStart,
                nNewPatches,
                fvMeshSubsetPtr_().boundaryMesh(),
                emptyPolyPatch::typeName
            );

            //Pout<< "    " << exposedPatchName << " : "
            //    << boundaryPatchSizes[oldInternalPatchID] << endl;

            // The index for the first patch is -1 as it originates from
            // the internal faces
            patchStart += boundaryPatchSizes[oldInternalPatchID];
            patchMap_[nNewPatches] = -1;
            ++nNewPatches;
        }
    }

    // Old patches

    for
    (
        label oldPatchi = nextPatchID;
        oldPatchi < oldPatches.size();
        oldPatchi++
    )
    {
        const label newSize = boundaryPatchSizes[globalPatchMap[oldPatchi]];

        // Patch still exists. Add it
        newBoundary[nNewPatches] = oldPatches[oldPatchi].clone
        (
            fvMeshSubsetPtr_().boundaryMesh(),
            nNewPatches,
            newSize,
            patchStart
        ).ptr();

        //Pout<< "    " << oldPatches[oldPatchi].name() << " : "
        //    << newSize << endl;

        patchStart += newSize;
        patchMap_[nNewPatches] = oldPatchi;    // compact patchMap
        ++nNewPatches;
    }


    // Reset the patch lists
    newBoundary.setSize(nNewPatches);
    patchMap_.setSize(nNewPatches);


    // Add the fvPatches
    fvMeshSubsetPtr_().addFvPatches(newBoundary, syncPar);

    // Subset and add any zones
    subsetZones();
}


void Foam::fvMeshSubset::setCellSubset
(
    const labelUList& selectedCells,
    const label patchID,
    const bool syncPar
)
{
    setCellSubset
    (
        BitSetOps::create(baseMesh().nCells(), selectedCells),
        patchID,
        syncPar
    );
}


void Foam::fvMeshSubset::setCellSubset
(
    const labelHashSet& selectedCells,
    const label patchID,
    const bool syncPar
)
{
    setCellSubset
    (
        BitSetOps::create(baseMesh().nCells(), selectedCells),
        patchID,
        syncPar
    );
}


void Foam::fvMeshSubset::setCellSubset
(
    const label regioni,
    const labelUList& regions,
    const label patchID,
    const bool syncPar
)
{
    setCellSubset
    (
        BitSetOps::create(baseMesh().nCells(), regioni, regions),
        patchID,
        syncPar
    );
}


Foam::labelList Foam::fvMeshSubset::getExposedFaces
(
    const bitSet& selectedCells,
    const bool syncCouples
) const
{
    return
        removeCells(baseMesh(), syncCouples)
       .getExposedFaces(getCellsToRemove(selectedCells));
}


Foam::labelList Foam::fvMeshSubset::getExposedFaces
(
    const label regioni,
    const labelUList& regions,
    const bool syncCouples
) const
{
    return
        removeCells(baseMesh(), syncCouples)
       .getExposedFaces(getCellsToRemove(regioni, regions));
}


void Foam::fvMeshSubset::setCellSubset
(
    const bitSet& selectedCells,
    const labelList& exposedFaces,
    const labelList& patchIDs,
    const bool syncCouples
)
{
    removeCellsImpl
    (
        getCellsToRemove(selectedCells),
        exposedFaces,
        patchIDs,
        syncCouples
    );
}


void Foam::fvMeshSubset::setCellSubset
(
    const label selectRegion,
    const labelList& regions,
    const labelList& exposedFaces,
    const labelList& patchIDs,
    const bool syncCouples
)
{
    removeCellsImpl
    (
        getCellsToRemove(selectRegion, regions),
        exposedFaces,
        patchIDs,
        syncCouples
    );
}


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