/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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-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/>.

Application
    splitMeshRegions

Group
    grpMeshManipulationUtilities

Description
    Splits mesh into multiple regions.

    Each region is defined as a domain whose cells can all be reached by
    cell-face-cell walking without crossing
    - boundary faces
    - additional faces from faceset (-blockedFaces faceSet).
    - any face between differing cellZones (-cellZones)

    Output is:
    - volScalarField with regions as different scalars (-detectOnly)
            or
    - mesh with multiple regions and mapped patches. These patches
      either cover the whole interface between two region (default) or
      only part according to faceZones (-useFaceZones)
            or
    - mesh with cells put into cellZones (-makeCellZones)

    Note:

    - multiple cellZones can be combined into a single cellZone (cluster)
    for further analysis using the 'combineZones' option:
        -combineZones '((zoneA "zoneB.*")(none otherZone))
    This can be combined with e.g. 'cellZones' or 'cellZonesOnly'. The
    corresponding region name will be the names of the cellZones combined e.g.
        zoneA_zoneB0_zoneB1

    - cellZonesOnly does not do a walk and uses the cellZones only. Use
    this if you don't mind having disconnected domains in a single region.
    This option requires all cells to be in one (and one only) cellZone.

    - cellZonesFileOnly behaves like -cellZonesOnly but reads the cellZones
    from the specified file. This allows one to explicitly specify the region
    distribution and still have multiple cellZones per region.

    - prefixRegion prefixes all normal patches with region name (interface
    (patches already have region name prefix)

    - Should work in parallel.
    cellZones can differ on either side of processor boundaries in which case
    the faces get moved from processor patch to mapped patch. Not very well
    tested.

    - If a cell zone gets split into more than one region it can detect
    the largest matching region (-sloppyCellZones). This will accept any
    region that covers more than 50% of the zone. It has to be a subset
    so cannot have any cells in any other zone.

    - If explicitly a single region has been selected (-largestOnly or
      -insidePoint) its region name will be either
        - name of a cellZone it matches to or
        - "largestOnly" respectively "insidePoint" or
        - polyMesh::defaultRegion if additionally -overwrite
          (so it will overwrite the input mesh!)

    - writes maps like decomposePar back to original mesh:
        - pointRegionAddressing : for every point in this region the point in
        the original mesh
        - cellRegionAddressing  :   ,,      cell                ,,  cell    ,,
        - faceRegionAddressing  :   ,,      face                ,,  face in
        the original mesh + 'turning index'. For a face in the same orientation
        this is the original facelabel+1, for a turned face this is -facelabel-1
        - boundaryRegionAddressing : for every patch in this region the
        patch in the original mesh (or -1 if added patch)

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

#include "SortableList.H"
#include "argList.H"
#include "regionSplit.H"
#include "fvMeshSubset.H"
#include "IOobjectList.H"
#include "volFields.H"
#include "faceSet.H"
#include "cellSet.H"
#include "polyTopoChange.H"
#include "removeCells.H"
#include "edgeHashes.H"
#include "syncTools.H"
#include "ReadFields.H"
#include "mappedWallPolyPatch.H"
#include "fvMeshTools.H"
#include "zeroGradientFvPatchFields.H"
#include "processorMeshes.H"

using namespace Foam;

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

// Prepend prefix to selected patches.
void renamePatches
(
    fvMesh& mesh,
    const word& prefix,
    const labelList& patchesToRename
)
{
    polyBoundaryMesh& polyPatches =
        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
    forAll(patchesToRename, i)
    {
        label patchi = patchesToRename[i];
        polyPatch& pp = polyPatches[patchi];

        if (isA<coupledPolyPatch>(pp))
        {
            WarningInFunction
                << "Encountered coupled patch " << pp.name()
                << ". Will only rename the patch itself,"
                << " not any referred patches."
                << " This might have to be done by hand."
                << endl;
        }

        pp.name() = prefix + '_' + pp.name();
    }
    // Recalculate any demand driven data (e.g. group to name lookup)
    polyPatches.updateMesh();
}


template<class GeoField>
void subsetVolFields
(
    const fvMesh& mesh,
    const fvMesh& subMesh,
    const labelList& cellMap,
    const labelList& faceMap,
    const labelHashSet& addedPatches
)
{
    const labelList patchMap(identity(mesh.boundaryMesh().size()));

    HashTable<const GeoField*> fields
    (
        mesh.objectRegistry::lookupClass<GeoField>()
    );
    forAllConstIters(fields, iter)
    {
        const GeoField& fld = *iter.val();

        Info<< "Mapping field " << fld.name() << endl;

        tmp<GeoField> tSubFld
        (
            fvMeshSubset::interpolate
            (
                fld,
                subMesh,
                patchMap,
                cellMap,
                faceMap
            )
        );

        // Hack: set value to 0 for introduced patches (since don't
        //       get initialised.
        forAll(tSubFld().boundaryField(), patchi)
        {
            if (addedPatches.found(patchi))
            {
                tSubFld.ref().boundaryFieldRef()[patchi] ==
                    typename GeoField::value_type(Zero);
            }
        }

        // Store on subMesh
        GeoField* subFld = tSubFld.ptr();
        subFld->rename(fld.name());
        subFld->writeOpt(IOobject::AUTO_WRITE);
        subFld->store();
    }
}


template<class GeoField>
void subsetSurfaceFields
(
    const fvMesh& mesh,
    const fvMesh& subMesh,
    const labelList& cellMap,
    const labelList& faceMap,
    const labelHashSet& addedPatches
)
{
    const labelList patchMap(identity(mesh.boundaryMesh().size()));

    HashTable<const GeoField*> fields
    (
        mesh.objectRegistry::lookupClass<GeoField>()
    );
    forAllConstIters(fields, iter)
    {
        const GeoField& fld = *iter.val();

        Info<< "Mapping field " << fld.name() << endl;

        tmp<GeoField> tSubFld
        (
            fvMeshSubset::interpolate
            (
                fld,
                subMesh,
                patchMap,
                cellMap,
                faceMap
            )
        );

        // Hack: set value to 0 for introduced patches (since don't
        //       get initialised.
        forAll(tSubFld().boundaryField(), patchi)
        {
            if (addedPatches.found(patchi))
            {
                tSubFld.ref().boundaryFieldRef()[patchi] ==
                    typename GeoField::value_type(Zero);
            }
        }

        // Store on subMesh
        GeoField* subFld = tSubFld.ptr();
        subFld->rename(fld.name());
        subFld->writeOpt(IOobject::AUTO_WRITE);
        subFld->store();
    }
}

// Select all cells not in the region
labelList getNonRegionCells(const labelList& cellRegion, const label regionI)
{
    DynamicList<label> nonRegionCells(cellRegion.size());
    forAll(cellRegion, celli)
    {
        if (cellRegion[celli] != regionI)
        {
            nonRegionCells.append(celli);
        }
    }
    return nonRegionCells.shrink();
}


void addToInterface
(
    const polyMesh& mesh,
    const label zoneID,
    const label ownRegion,
    const label neiRegion,
    EdgeMap<Map<label>>& regionsToSize
)
{
    edge interface
    (
        min(ownRegion, neiRegion),
        max(ownRegion, neiRegion)
    );

    auto iter = regionsToSize.find(interface);

    if (iter.found())
    {
        // Check if zone present
        auto zoneIter = iter().find(zoneID);
        if (zoneIter.found())
        {
            // Found zone. Increment count.
            ++(*zoneIter);
        }
        else
        {
            // New or no zone. Insert with count 1.
            iter().insert(zoneID, 1);
        }
    }
    else
    {
        // Create new interface of size 1.
        Map<label> zoneToSize;
        zoneToSize.insert(zoneID, 1);
        regionsToSize.insert(interface, zoneToSize);
    }
}


// Get region-region interface name and sizes.
// Returns interfaces as straight list for looping in identical order.
void getInterfaceSizes
(
    const polyMesh& mesh,
    const bool useFaceZones,
    const labelList& cellRegion,
    const wordList& regionNames,

    edgeList& interfaces,
    List<Pair<word>>& interfaceNames,
    labelList& interfaceSizes,
    labelList& faceToInterface
)
{
    // From region-region to faceZone (or -1) to number of faces.

    EdgeMap<Map<label>> regionsToSize;


    // Internal faces
    // ~~~~~~~~~~~~~~

    forAll(mesh.faceNeighbour(), facei)
    {
        label ownRegion = cellRegion[mesh.faceOwner()[facei]];
        label neiRegion = cellRegion[mesh.faceNeighbour()[facei]];

        if (ownRegion != neiRegion)
        {
            addToInterface
            (
                mesh,
                (useFaceZones ? mesh.faceZones().whichZone(facei) : -1),
                ownRegion,
                neiRegion,
                regionsToSize
            );
        }
    }

    // Boundary faces
    // ~~~~~~~~~~~~~~

    // Neighbour cellRegion.
    labelList coupledRegion(mesh.nBoundaryFaces());

    forAll(coupledRegion, i)
    {
        label celli = mesh.faceOwner()[i+mesh.nInternalFaces()];
        coupledRegion[i] = cellRegion[celli];
    }
    syncTools::swapBoundaryFaceList(mesh, coupledRegion);

    forAll(coupledRegion, i)
    {
        label facei = i+mesh.nInternalFaces();
        label ownRegion = cellRegion[mesh.faceOwner()[facei]];
        label neiRegion = coupledRegion[i];

        if (ownRegion != neiRegion)
        {
            addToInterface
            (
                mesh,
                (useFaceZones ? mesh.faceZones().whichZone(facei) : -1),
                ownRegion,
                neiRegion,
                regionsToSize
            );
        }
    }


    if (Pstream::parRun())
    {
        if (Pstream::master())
        {
            // Receive and add to my sizes
            for (const int slave : Pstream::subProcs())
            {
                IPstream fromSlave(Pstream::commsTypes::blocking, slave);

                EdgeMap<Map<label>> slaveSizes(fromSlave);

                forAllConstIters(slaveSizes, slaveIter)
                {
                    const Map<label>& slaveInfo = *slaveIter;

                    auto masterIter = regionsToSize.find(slaveIter.key());

                    if (masterIter.found())
                    {
                        // Same inter-region
                        Map<label>& masterInfo = *masterIter;

                        forAllConstIters(slaveInfo, iter)
                        {
                            const label zoneID = iter.key();
                            const label slaveSize = iter.val();

                            auto zoneIter = masterInfo.find(zoneID);
                            if (zoneIter.found())
                            {
                                *zoneIter += slaveSize;
                            }
                            else
                            {
                                masterInfo.insert(zoneID, slaveSize);
                            }
                        }
                    }
                    else
                    {
                        regionsToSize.insert(slaveIter.key(), slaveInfo);
                    }
                }
            }
        }
        else
        {
            // Send to master
            {
                OPstream toMaster
                (
                    Pstream::commsTypes::blocking,
                    Pstream::masterNo()
                );
                toMaster << regionsToSize;
            }
        }
    }

    // Rework

    Pstream::scatter(regionsToSize);



    // Now we have the global sizes of all inter-regions.
    // Invert this on master and distribute.
    label nInterfaces = 0;
    forAllConstIters(regionsToSize, iter)
    {
        const Map<label>& info = iter.val();
        nInterfaces += info.size();
    }

    interfaces.setSize(nInterfaces);
    interfaceNames.setSize(nInterfaces);
    interfaceSizes.setSize(nInterfaces);
    EdgeMap<Map<label>> regionsToInterface(nInterfaces);

    nInterfaces = 0;
    forAllConstIters(regionsToSize, iter)
    {
        const edge& e = iter.key();
        const Map<label>& info = iter.val();

        const word& name0 = regionNames[e[0]];
        const word& name1 = regionNames[e[1]];

        forAllConstIters(info, infoIter)
        {
            interfaces[nInterfaces] = iter.key();
            label zoneID = infoIter.key();
            if (zoneID == -1)
            {
                interfaceNames[nInterfaces] = Pair<word>
                (
                    name0 + "_to_" + name1,
                    name1 + "_to_" + name0
                );
            }
            else
            {
                const word& zoneName = mesh.faceZones()[zoneID].name();
                interfaceNames[nInterfaces] = Pair<word>
                (
                    zoneName + "_" + name0 + "_to_" + name1,
                    zoneName + "_" + name1 + "_to_" + name0
                );
            }
            interfaceSizes[nInterfaces] = infoIter();

            if (regionsToInterface.found(e))
            {
                regionsToInterface[e].insert(zoneID, nInterfaces);
            }
            else
            {
                Map<label> zoneAndInterface;
                zoneAndInterface.insert(zoneID, nInterfaces);
                regionsToInterface.insert(e, zoneAndInterface);
            }
            nInterfaces++;
        }
    }


    // Now all processor have consistent interface information

    Pstream::scatter(interfaces);
    Pstream::scatter(interfaceNames);
    Pstream::scatter(interfaceSizes);
    Pstream::scatter(regionsToInterface);

    // Mark all inter-region faces.
    faceToInterface.setSize(mesh.nFaces(), -1);

    forAll(mesh.faceNeighbour(), facei)
    {
        label ownRegion = cellRegion[mesh.faceOwner()[facei]];
        label neiRegion = cellRegion[mesh.faceNeighbour()[facei]];

        if (ownRegion != neiRegion)
        {
            label zoneID = -1;
            if (useFaceZones)
            {
                zoneID = mesh.faceZones().whichZone(facei);
            }

            edge interface
            (
                min(ownRegion, neiRegion),
                max(ownRegion, neiRegion)
            );

            faceToInterface[facei] = regionsToInterface[interface][zoneID];
        }
    }
    forAll(coupledRegion, i)
    {
        label facei = i+mesh.nInternalFaces();
        label ownRegion = cellRegion[mesh.faceOwner()[facei]];
        label neiRegion = coupledRegion[i];

        if (ownRegion != neiRegion)
        {
            label zoneID = -1;
            if (useFaceZones)
            {
                zoneID = mesh.faceZones().whichZone(facei);
            }

            edge interface
            (
                min(ownRegion, neiRegion),
                max(ownRegion, neiRegion)
            );

            faceToInterface[facei] = regionsToInterface[interface][zoneID];
        }
    }
}


// Create mesh for region.
autoPtr<mapPolyMesh> createRegionMesh
(
    const fvMesh& mesh,
    // Region info
    const labelList& cellRegion,
    const label regionI,
    const word& regionName,
    // Interface info
    const labelList& interfacePatches,
    const labelList& faceToInterface,

    autoPtr<fvMesh>& newMesh
)
{
    // Create dummy system/fv*
    fvMeshTools::createDummyFvMeshFiles(mesh, regionName, true);

    // Neighbour cellRegion.
    labelList coupledRegion(mesh.nBoundaryFaces());

    forAll(coupledRegion, i)
    {
        label celli = mesh.faceOwner()[i+mesh.nInternalFaces()];
        coupledRegion[i] = cellRegion[celli];
    }
    syncTools::swapBoundaryFaceList(mesh, coupledRegion);


    // Topology change container. Start off from existing mesh.
    polyTopoChange meshMod(mesh);

    // Cell remover engine
    removeCells cellRemover(mesh);

    // Select all but region cells
    labelList cellsToRemove(getNonRegionCells(cellRegion, regionI));

    // Find out which faces will get exposed. Note that this
    // gets faces in mesh face order. So both regions will get same
    // face in same order (important!)
    labelList exposedFaces = cellRemover.getExposedFaces(cellsToRemove);

    labelList exposedPatchIDs(exposedFaces.size());
    forAll(exposedFaces, i)
    {
        label facei = exposedFaces[i];
        label interfacei = faceToInterface[facei];

        label ownRegion = cellRegion[mesh.faceOwner()[facei]];
        label neiRegion = -1;

        if (mesh.isInternalFace(facei))
        {
            neiRegion = cellRegion[mesh.faceNeighbour()[facei]];
        }
        else
        {
            neiRegion = coupledRegion[facei-mesh.nInternalFaces()];
        }


        // Check which side is being kept - determines which of the two
        // patches will be used.

        label otherRegion = -1;

        if (ownRegion == regionI && neiRegion != regionI)
        {
            otherRegion = neiRegion;
        }
        else if (ownRegion != regionI && neiRegion == regionI)
        {
            otherRegion = ownRegion;
        }
        else
        {
            FatalErrorInFunction
                << "Exposed face:" << facei
                << " fc:" << mesh.faceCentres()[facei]
                << " has owner region " << ownRegion
                << " and neighbour region " << neiRegion
                << " when handling region:" << regionI
                << exit(FatalError);
        }

        // Find the patch.
        if (regionI < otherRegion)
        {
            exposedPatchIDs[i] = interfacePatches[interfacei];
        }
        else
        {
            exposedPatchIDs[i] = interfacePatches[interfacei]+1;
        }
    }

    // Remove faces
    cellRemover.setRefinement
    (
        cellsToRemove,
        exposedFaces,
        exposedPatchIDs,
        meshMod
    );

    autoPtr<mapPolyMesh> map = meshMod.makeMesh
    (
        newMesh,
        IOobject
        (
            regionName,
            mesh.time().timeName(),
            mesh.time(),
            IOobject::READ_IF_PRESENT,  // read fv* if present
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    return map;
}


void createAndWriteRegion
(
    const fvMesh& mesh,
    const labelList& cellRegion,
    const wordList& regionNames,
    const bool prefixRegion,
    const labelList& faceToInterface,
    const labelList& interfacePatches,
    const label regionI,
    const word& newMeshInstance
)
{
    Info<< "Creating mesh for region " << regionI
        << ' ' << regionNames[regionI] << endl;

    autoPtr<fvMesh> newMesh;
    autoPtr<mapPolyMesh> map = createRegionMesh
    (
        mesh,
        cellRegion,
        regionI,
        regionNames[regionI],
        interfacePatches,
        faceToInterface,
        newMesh
    );


    // Make map of all added patches
    labelHashSet addedPatches(2*interfacePatches.size());
    forAll(interfacePatches, interfacei)
    {
        addedPatches.insert(interfacePatches[interfacei]);
        addedPatches.insert(interfacePatches[interfacei]+1);
    }


    Info<< "Mapping fields" << endl;

    // Map existing fields
    newMesh().updateMesh(map());

    // Add subsetted fields
    subsetVolFields<volScalarField>
    (
        mesh,
        newMesh(),
        map().cellMap(),
        map().faceMap(),
        addedPatches
    );
    subsetVolFields<volVectorField>
    (
        mesh,
        newMesh(),
        map().cellMap(),
        map().faceMap(),
        addedPatches
    );
    subsetVolFields<volSphericalTensorField>
    (
        mesh,
        newMesh(),
        map().cellMap(),
        map().faceMap(),
        addedPatches
    );
    subsetVolFields<volSymmTensorField>
    (
        mesh,
        newMesh(),
        map().cellMap(),
        map().faceMap(),
        addedPatches
    );
    subsetVolFields<volTensorField>
    (
        mesh,
        newMesh(),
        map().cellMap(),
        map().faceMap(),
        addedPatches
    );

    subsetSurfaceFields<surfaceScalarField>
    (
        mesh,
        newMesh(),
        map().cellMap(),
        map().faceMap(),
        addedPatches
    );
    subsetSurfaceFields<surfaceVectorField>
    (
        mesh,
        newMesh(),
        map().cellMap(),
        map().faceMap(),
        addedPatches
    );
    subsetSurfaceFields<surfaceSphericalTensorField>
    (
        mesh,
        newMesh(),
        map().cellMap(),
        map().faceMap(),
        addedPatches
    );
    subsetSurfaceFields<surfaceSymmTensorField>
    (
        mesh,
        newMesh(),
        map().cellMap(),
        map().faceMap(),
        addedPatches
    );
    subsetSurfaceFields<surfaceTensorField>
    (
        mesh,
        newMesh(),
        map().cellMap(),
        map().faceMap(),
        addedPatches
    );


    const polyBoundaryMesh& newPatches = newMesh().boundaryMesh();
    newPatches.checkParallelSync(true);

    // Delete empty patches
    // ~~~~~~~~~~~~~~~~~~~~

    // Create reordering list to move patches-to-be-deleted to end
    labelList oldToNew(newPatches.size(), -1);
    DynamicList<label> sharedPatches(newPatches.size());
    label newI = 0;

    Info<< "Deleting empty patches" << endl;

    // Assumes all non-proc boundaries are on all processors!
    forAll(newPatches, patchi)
    {
        const polyPatch& pp = newPatches[patchi];

        if (!isA<processorPolyPatch>(pp))
        {
            if (returnReduce(pp.size(), sumOp<label>()) > 0)
            {
                oldToNew[patchi] = newI;
                if (!addedPatches.found(patchi))
                {
                    sharedPatches.append(newI);
                }
                newI++;
            }
        }
    }

    // Same for processor patches (but need no reduction)
    forAll(newPatches, patchi)
    {
        const polyPatch& pp = newPatches[patchi];

        if (isA<processorPolyPatch>(pp) && pp.size())
        {
            oldToNew[patchi] = newI++;
        }
    }

    const label nNewPatches = newI;

    // Move all deleteable patches to the end
    forAll(oldToNew, patchi)
    {
        if (oldToNew[patchi] == -1)
        {
            oldToNew[patchi] = newI++;
        }
    }

    //reorderPatches(newMesh(), oldToNew, nNewPatches);
    fvMeshTools::reorderPatches(newMesh(), oldToNew, nNewPatches, true);

    // Rename shared patches with region name
    if (prefixRegion)
    {
        Info<< "Prefixing patches with region name" << endl;

        renamePatches(newMesh(), regionNames[regionI], sharedPatches);
    }


    Info<< "Writing new mesh" << endl;

    newMesh().setInstance(newMeshInstance);
    newMesh().write();
    topoSet::removeFiles(newMesh());
    processorMeshes::removeFiles(newMesh());

    // Write addressing files like decomposePar
    Info<< "Writing addressing to base mesh" << endl;

    labelIOList pointProcAddressing
    (
        IOobject
        (
            "pointRegionAddressing",
            newMesh().facesInstance(),
            newMesh().meshSubDir,
            newMesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        map().pointMap()
    );
    Info<< "Writing map " << pointProcAddressing.name()
        << " from region" << regionI
        << " points back to base mesh." << endl;
    pointProcAddressing.write();

    labelIOList faceProcAddressing
    (
        IOobject
        (
            "faceRegionAddressing",
            newMesh().facesInstance(),
            newMesh().meshSubDir,
            newMesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        newMesh().nFaces()
    );
    forAll(faceProcAddressing, facei)
    {
        // face + turning index. (see decomposePar)
        // Is the face pointing in the same direction?
        label oldFacei = map().faceMap()[facei];

        if
        (
            map().cellMap()[newMesh().faceOwner()[facei]]
         == mesh.faceOwner()[oldFacei]
        )
        {
            faceProcAddressing[facei] = oldFacei+1;
        }
        else
        {
            faceProcAddressing[facei] = -(oldFacei+1);
        }
    }
    Info<< "Writing map " << faceProcAddressing.name()
        << " from region" << regionI
        << " faces back to base mesh." << endl;
    faceProcAddressing.write();

    labelIOList cellProcAddressing
    (
        IOobject
        (
            "cellRegionAddressing",
            newMesh().facesInstance(),
            newMesh().meshSubDir,
            newMesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        map().cellMap()
    );
    Info<< "Writing map " <<cellProcAddressing.name()
        << " from region" << regionI
        << " cells back to base mesh." << endl;
    cellProcAddressing.write();

    labelIOList boundaryProcAddressing
    (
        IOobject
        (
            "boundaryRegionAddressing",
            newMesh().facesInstance(),
            newMesh().meshSubDir,
            newMesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        labelList(nNewPatches, -1)
    );
    forAll(oldToNew, i)
    {
        if (!addedPatches.found(i))
        {
            label newI = oldToNew[i];
            if (newI >= 0 && newI < nNewPatches)
            {
                boundaryProcAddressing[oldToNew[i]] = i;
            }
        }
    }
    Info<< "Writing map " << boundaryProcAddressing.name()
        << " from region" << regionI
        << " boundary back to base mesh." << endl;
    boundaryProcAddressing.write();
}


// Create for every region-region interface with non-zero size two patches.
// First one is for minimumregion to maximumregion.
// Note that patches get created in same order on all processors (if parallel)
// since looping over synchronised 'interfaces'.
labelList addRegionPatches
(
    fvMesh& mesh,
    const wordList& regionNames,
    const edgeList& interfaces,
    const List<Pair<word>>& interfaceNames
)
{
    Info<< nl << "Adding patches" << nl << endl;

    labelList interfacePatches(interfaces.size());

    forAll(interfaces, interI)
    {
        const edge& e = interfaces[interI];
        const Pair<word>& names = interfaceNames[interI];

        //Info<< "For interface " << interI
        //    << " between regions " << e
        //    << " trying to add patches " << names << endl;


        mappedWallPolyPatch patch1
        (
            names[0],
            0,                  // overridden
            0,                  // overridden
            0,                  // overridden
            regionNames[e[1]],  // sampleRegion
            mappedPatchBase::NEARESTPATCHFACE,
            names[1],           // samplePatch
            point::zero,        // offset
            mesh.boundaryMesh()
        );

        interfacePatches[interI] = fvMeshTools::addPatch
        (
            mesh,
            patch1,
            dictionary(),   //optional per field value
            calculatedFvPatchField<scalar>::typeName,
            true            //validBoundary
        );

        mappedWallPolyPatch patch2
        (
            names[1],
            0,
            0,
            0,
            regionNames[e[0]],  // sampleRegion
            mappedPatchBase::NEARESTPATCHFACE,
            names[0],
            point::zero,        // offset
            mesh.boundaryMesh()
        );
        fvMeshTools::addPatch
        (
            mesh,
            patch2,
            dictionary(),   //optional per field value
            calculatedFvPatchField<scalar>::typeName,
            true            //validBoundary
        );

        Info<< "For interface between region " << regionNames[e[0]]
            << " and " << regionNames[e[1]] << " added patches" << endl
            << "    " << interfacePatches[interI]
            << "\t" << mesh.boundaryMesh()[interfacePatches[interI]].name()
            << endl
            << "    " << interfacePatches[interI]+1
            << "\t" << mesh.boundaryMesh()[interfacePatches[interI]+1].name()
            << endl;
    }
    return interfacePatches;
}


// Find region that covers most of cell zone
label findCorrespondingRegion
(
    const labelList& existingZoneID,    // per cell the (unique) zoneID
    const labelList& cellRegion,
    const label nCellRegions,
    const label zoneI,
    const label minOverlapSize
)
{
    // Per region the number of cells in zoneI
    labelList cellsInZone(nCellRegions, Zero);

    forAll(cellRegion, celli)
    {
        if (existingZoneID[celli] == zoneI)
        {
            cellsInZone[cellRegion[celli]]++;
        }
    }

    Pstream::listCombineGather(cellsInZone, plusEqOp<label>());
    Pstream::listCombineScatter(cellsInZone);

    // Pick region with largest overlap of zoneI
    label regionI = findMax(cellsInZone);


    if (cellsInZone[regionI] < minOverlapSize)
    {
        // Region covers too little of zone. Not good enough.
        regionI = -1;
    }
    else
    {
        // Check that region contains no cells that aren't in cellZone.
        forAll(cellRegion, celli)
        {
            if (cellRegion[celli] == regionI && existingZoneID[celli] != zoneI)
            {
                // celli in regionI but not in zoneI
                regionI = -1;
                break;
            }
        }
        // If one in error, all should be in error. Note that branch gets taken
        // on all procs.
        reduce(regionI, minOp<label>());
    }

    return regionI;
}


void getClusterID
(
    const polyMesh& mesh,
    const cellZoneMesh& cellZones,
    const wordList& clusterNames,
    const labelListList& clusterToZones,
    labelList& clusterID,
    labelList& neiClusterID
)
{
    // Existing zoneID
    clusterID.setSize(mesh.nCells());
    clusterID = -1;

    forAll(clusterToZones, clusterI)
    {
        for (const label zoneI : clusterToZones[clusterI])
        {
            const cellZone& cz = cellZones[zoneI];

            forAll(cz, i)
            {
                label celli = cz[i];
                if (clusterID[celli] == -1)
                {
                    clusterID[celli] = clusterI;
                }
                else
                {
                    FatalErrorInFunction
                        << "Cell " << celli << " with cell centre "
                        << mesh.cellCentres()[celli]
                        << " is multiple zones. This is not allowed." << endl
                        << "It is in zone " << clusterNames[clusterID[celli]]
                        << " and in zone " << clusterNames[clusterI]
                        << exit(FatalError);
                }
            }
        }
    }

    // Neighbour zoneID.
    syncTools::swapBoundaryCellList(mesh, clusterID, neiClusterID);
}


word makeRegionName
(
    const cellZoneMesh& czs,
    const label regioni,
    const labelList& zoneIDs
)
{
    // Synthesise region name. Equals the zone name if cluster consist of only
    // one zone

    if (zoneIDs.empty())
    {
        return word("domain") + Foam::name(regioni);
    }
    else
    {
        // Synthesize name from appended cellZone names
        word regionName(czs[zoneIDs[0]].name());
        for (label i = 1; i < zoneIDs.size(); i++)
        {
            regionName += "_" + czs[zoneIDs[i]].name();
        }
        return regionName;
    }
}


void makeClusters
(
    const List<wordRes>& zoneClusters,
    const cellZoneMesh& cellZones,
    wordList& clusterNames,
    labelListList& clusterToZones,
    labelList& zoneToCluster
)
{
    // Check if there are clustering for zones. If none every zone goes into
    // its own cluster.

    clusterNames.clear();
    clusterToZones.clear();
    zoneToCluster.setSize(cellZones.size());
    zoneToCluster = -1;

    if (zoneClusters.size())
    {
        forAll(zoneClusters, clusteri)
        {
            const labelList zoneIDs(cellZones.indices(zoneClusters[clusteri]));
            UIndirectList<label>(zoneToCluster, zoneIDs) = clusteri;
            clusterNames.append(makeRegionName(cellZones, clusteri, zoneIDs));
            clusterToZones.append(std::move(zoneIDs));
        }

        // Unclustered zone
        forAll(zoneToCluster, zonei)
        {
            if (zoneToCluster[zonei] == -1)
            {
                clusterNames.append(cellZones[zonei].name());
                clusterToZones.append(labelList(1, zonei));
                zoneToCluster[zonei] = clusterToZones.size();
            }
        }
    }
    else
    {
        for (const auto& cellZone : cellZones)
        {
            const label nClusters = clusterToZones.size();
            clusterNames.append(cellZone.name());
            clusterToZones.append(labelList(1, cellZone.index()));
            zoneToCluster[cellZone.index()] = nClusters;
        }
    }
}


void matchRegions
(
    const bool sloppyCellZones,
    const polyMesh& mesh,

    const wordList& clusterNames,
    const labelListList& clusterToZones,
    const labelList& clusterID,

    const label nCellRegions,
    const labelList& cellRegion,

    labelListList& regionToZones,
    wordList& regionNames,
    labelList& zoneToRegion
)
{
    const cellZoneMesh& cellZones = mesh.cellZones();

    regionToZones.setSize(nCellRegions);
    regionToZones = labelList();
    regionNames.setSize(nCellRegions);
    regionNames = word::null;
    zoneToRegion.setSize(cellZones.size(), -1);


    // Sizes per cluster
    labelList clusterSizes(clusterToZones.size(), Zero);
    forAll(clusterToZones, clusterI)
    {
        for (const label zoneI : clusterToZones[clusterI])
        {
            clusterSizes[clusterI] += cellZones[zoneI].size();
        }
        reduce(clusterSizes[clusterI], sumOp<label>());
    }

    if (sloppyCellZones)
    {
        Info<< "Trying to match regions to existing cell zones;"
            << " region can be subset of cell zone." << nl << endl;

        forAll(clusterToZones, clusterI)
        {
            label regionI = findCorrespondingRegion
            (
                clusterID,
                cellRegion,
                nCellRegions,
                clusterI,
                label(0.5*clusterSizes[clusterI]) // minimum overlap
            );

            if (regionI != -1)
            {
                Info<< "Sloppily matched region " << regionI
                    //<< " size " << regionSizes[regionI]
                    << " to cluster " << clusterI
                    << " size " << clusterSizes[clusterI]
                    << endl;
                UIndirectList<label>
                (
                    zoneToRegion,
                    clusterToZones[clusterI]
                ) = regionI;
                regionToZones[regionI] = clusterToZones[clusterI];
                regionNames[regionI] = clusterNames[clusterI];
            }
        }
    }
    else
    {
        Info<< "Trying to match regions to existing cell zones." << nl << endl;

        forAll(clusterToZones, clusterI)
        {
            label regionI = findCorrespondingRegion
            (
                clusterID,
                cellRegion,
                nCellRegions,
                clusterI,
                clusterSizes[clusterI]  // require exact match
            );

            if (regionI != -1)
            {
                UIndirectList<label>
                (
                    zoneToRegion,
                    clusterToZones[clusterI]
                ) = regionI;
                regionToZones[regionI] = clusterToZones[clusterI];
                regionNames[regionI] = clusterNames[clusterI];
            }
        }
    }
    // Allocate region names for unmatched regions.
    forAll(regionNames, regionI)
    {
        regionNames[regionI] = makeRegionName
        (
            cellZones,
            regionI,
            regionToZones[regionI]
        );
    }
}


void writeCellToRegion(const fvMesh& mesh, const labelList& cellRegion)
{
    // Write to manual decomposition option
    {
        labelIOList cellToRegion
        (
            IOobject
            (
                "cellToRegion",
                mesh.facesInstance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            cellRegion
        );
        cellToRegion.write();

        Info<< "Writing region per cell file (for manual decomposition) to "
            << cellToRegion.objectPath() << nl << endl;
    }
    // Write for postprocessing
    {
        volScalarField cellToRegion
        (
            IOobject
            (
                "cellToRegion",
                mesh.time().timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            mesh,
            dimensionedScalar(dimless, Zero),
            zeroGradientFvPatchScalarField::typeName
        );
        forAll(cellRegion, celli)
        {
            cellToRegion[celli] = cellRegion[celli];
        }
        cellToRegion.write();

        Info<< "Writing region per cell as volScalarField to "
            << cellToRegion.objectPath() << nl << endl;
    }
}




int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Split mesh into multiple regions (detected by walking across faces)"
    );
    #include "addRegionOption.H"
    #include "addOverwriteOption.H"
    argList::addBoolOption
    (
        "cellZones",
        "Additionally split cellZones off into separate regions"
    );
    argList::addBoolOption
    (
        "cellZonesOnly",
        "Use cellZones only to split mesh into regions; do not use walking"
    );
    argList::addOption
    (
        "cellZonesFileOnly",
        "file",
        "Like -cellZonesOnly, but use specified file"
    );
    argList::addOption
    (
        "combineZones",
        "lists of zones",
        "Combine zones in follow-on analysis"
    );
    argList::addOption
    (
        "blockedFaces",
        "faceSet",
        "Specify additional region boundaries that walking does not cross"
    );
    argList::addBoolOption
    (
        "makeCellZones",
        "Place cells into cellZones instead of splitting mesh"
    );
    argList::addBoolOption
    (
        "largestOnly",
        "Only write largest region"
    );
    argList::addOption
    (
        "insidePoint",
        "point",
        "Only write region containing point"
    );
    argList::addBoolOption
    (
        "detectOnly",
        "Do not write mesh"
    );
    argList::addBoolOption
    (
        "sloppyCellZones",
        "Try to match heuristically regions to existing cell zones"
    );
    argList::addBoolOption
    (
        "useFaceZones",
        "Use faceZones to patch inter-region faces instead of single patch"
    );
    argList::addBoolOption
    (
        "prefixRegion",
        "Prefix region name to all patches, not just coupling patches"
    );

    argList::noFunctionObjects();  // Never use function objects

    #include "setRootCase.H"
    #include "createTime.H"
    #include "createNamedMesh.H"

    const word oldInstance = mesh.pointsInstance();

    word blockedFacesName;
    if (args.readIfPresent("blockedFaces", blockedFacesName))
    {
        Info<< "Reading blocked internal faces from faceSet "
            << blockedFacesName << nl << endl;
    }

    const bool makeCellZones    = args.found("makeCellZones");
    const bool largestOnly      = args.found("largestOnly");
    const bool insidePoint      = args.found("insidePoint");
    const bool useCellZones     = args.found("cellZones");
    const bool useCellZonesOnly = args.found("cellZonesOnly");
    const bool useCellZonesFile = args.found("cellZonesFileOnly");
    const bool combineZones     = args.found("combineZones");
    const bool overwrite        = args.found("overwrite");
    const bool detectOnly       = args.found("detectOnly");
    const bool sloppyCellZones  = args.found("sloppyCellZones");
    const bool useFaceZones     = args.found("useFaceZones");
    const bool prefixRegion     = args.found("prefixRegion");


    if
    (
        (useCellZonesOnly || useCellZonesFile)
     && (useCellZones || blockedFacesName.size())
    )
    {
        FatalErrorInFunction
            << "You cannot specify both -cellZonesOnly or -cellZonesFileOnly"
            << " (which specify complete split)"
            << " in combination with -blockedFaces or -cellZones"
            << " (which imply a split based on topology)"
            << exit(FatalError);
    }


    if (useFaceZones)
    {
        Info<< "Using current faceZones to divide inter-region interfaces"
            << " into multiple patches."
            << nl << endl;
    }
    else
    {
        Info<< "Creating single patch per inter-region interface."
            << nl << endl;
    }



    if (insidePoint && largestOnly)
    {
        FatalErrorInFunction
            << "You cannot specify both -largestOnly"
            << " (keep region with most cells)"
            << " and -insidePoint (keep region containing point)"
            << exit(FatalError);
    }


    // Make sure cellZone names consistent across processors
    mesh.cellZones().checkParallelSync(true);

    List<wordRes> zoneClusters;
    if (combineZones)
    {
        zoneClusters = args.get<List<wordRes>>("combineZones");
    }


    // Determine per cell the region it belongs to
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    // cellRegion is the labelList with the region per cell.
    labelList cellRegion;
    // Region to zone(s)
    labelListList regionToZones;
    // Name of region
    wordList regionNames;
    // Zone to region
    labelList zoneToRegion;

    label nCellRegions = 0;
    if (useCellZonesOnly)
    {
        Info<< "Using current cellZones to split mesh into regions."
            << " This requires all"
            << " cells to be in one and only one cellZone." << nl << endl;

        // Collect sets of zones into clusters. If no cluster is just an identity
        // list (cluster 0 is cellZone 0 etc.)
        wordList clusterNames;
        labelListList clusterToZones;
        labelList zoneToCluster;
        makeClusters
        (
            zoneClusters,
            mesh.cellZones(),
            clusterNames,
            clusterToZones,
            zoneToCluster
        );

        // Existing clusterID
        labelList clusterID(mesh.nCells(), -1);
        // Neighbour clusterID.
        labelList neiClusterID(mesh.nBoundaryFaces());
        getClusterID
        (
            mesh,
            mesh.cellZones(),
            clusterNames,
            clusterToZones,
            clusterID,
            neiClusterID
        );

        label unzonedCelli = clusterID.find(-1);
        if (unzonedCelli != -1)
        {
            FatalErrorInFunction
                << "For the cellZonesOnly option all cells "
                << "have to be in a cellZone." << endl
                << "Cell " << unzonedCelli
                << " at" << mesh.cellCentres()[unzonedCelli]
                << " is not in a cellZone. There might be more unzoned cells."
                << exit(FatalError);
        }
        cellRegion = clusterID;
        nCellRegions = gMax(cellRegion)+1;
        zoneToRegion = zoneToCluster;
        regionToZones = clusterToZones;
        regionNames = clusterNames;
    }
    else if (useCellZonesFile)
    {
        const word zoneFile(args["cellZonesFileOnly"]);
        Info<< "Reading split from cellZones file " << zoneFile << endl
            << "This requires all"
            << " cells to be in one and only one cellZone." << nl << endl;

        cellZoneMesh newCellZones
        (
            IOobject
            (
                zoneFile,
                mesh.facesInstance(),
                polyMesh::meshSubDir,
                mesh,
                IOobject::MUST_READ,
                IOobject::NO_WRITE,
                false
            ),
            mesh
        );

        wordList clusterNames;
        labelListList clusterToZones;
        labelList zoneToCluster;
        makeClusters
        (
            zoneClusters,
            newCellZones,
            clusterNames,
            clusterToZones,
            zoneToCluster
        );


        // Existing clusterID
        labelList clusterID(mesh.nCells(), -1);
        // Neighbour clusterID.
        labelList neiClusterID(mesh.nBoundaryFaces());
        getClusterID
        (
            mesh,
            newCellZones,
            clusterNames,
            clusterToZones,
            clusterID,
            neiClusterID
        );


        label unzonedCelli = clusterID.find(-1);
        if (unzonedCelli != -1)
        {
            FatalErrorInFunction
                << "For the cellZonesFileOnly option all cells "
                << "have to be in a cellZone." << endl
                << "Cell " << unzonedCelli
                << " at" << mesh.cellCentres()[unzonedCelli]
                << " is not in a cellZone. There might be more unzoned cells."
                << exit(FatalError);
        }
        cellRegion = clusterID;
        nCellRegions = gMax(cellRegion)+1;
        zoneToRegion = zoneToCluster;
        regionToZones = clusterToZones;
        regionNames = clusterNames;
    }
    else
    {
        // Determine connected regions
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~

        // Mark additional faces that are blocked
        boolList blockedFace;

        // Read from faceSet
        if (blockedFacesName.size())
        {
            faceSet blockedFaceSet(mesh, blockedFacesName);
            Info<< "Read "
                << returnReduce(blockedFaceSet.size(), sumOp<label>())
                << " blocked faces from set " << blockedFacesName << nl << endl;

            blockedFace.setSize(mesh.nFaces(), false);

            for (const label facei : blockedFaceSet)
            {
                blockedFace[facei] = true;
            }
        }

        // Collect sets of zones into clusters. If no cluster is just an
        // identity list (cluster 0 is cellZone 0 etc.)
        wordList clusterNames;
        labelListList clusterToZones;
        labelList zoneToCluster;
        makeClusters
        (
            zoneClusters,
            mesh.cellZones(),
            clusterNames,
            clusterToZones,
            zoneToCluster
        );

        // Existing clusterID
        labelList clusterID(mesh.nCells(), -1);
        // Neighbour clusterID.
        labelList neiClusterID(mesh.nBoundaryFaces());
        getClusterID
        (
            mesh,
            mesh.cellZones(),
            clusterNames,
            clusterToZones,
            clusterID,
            neiClusterID
        );


        // Imply from differing cellZones
        if (useCellZones)
        {
            blockedFace.setSize(mesh.nFaces(), false);

            for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
            {
                label ownCluster = clusterID[mesh.faceOwner()[facei]];
                label neiCluster = clusterID[mesh.faceNeighbour()[facei]];
                
                if (ownCluster != neiCluster)
                {
                    blockedFace[facei] = true;
                }
            }

            // Different cellZones on either side of processor patch.
            forAll(neiClusterID, i)
            {
                label facei = i+mesh.nInternalFaces();
                label ownCluster = clusterID[mesh.faceOwner()[facei]];
                label neiCluster = neiClusterID[i];

                if (ownCluster != neiCluster)
                {
                    blockedFace[facei] = true;
                }
            }
        }

        // Do a topological walk to determine regions
        regionSplit regions(mesh, blockedFace);
        nCellRegions = regions.nRegions();
        cellRegion.transfer(regions);

        // Match regions to zones
        matchRegions
        (
            sloppyCellZones,
            mesh,

            // cluster-to-zone and cluster-to-cell addressing
            clusterNames,
            clusterToZones,
            clusterID,

            // cell-to-region addressing
            nCellRegions,
            cellRegion,

            // matched zones
            regionToZones,
            regionNames,
            zoneToRegion
        );

        // Override any default region names if single region selected
        if (largestOnly || insidePoint)
        {
            forAll(regionToZones, regionI)
            {
                if (regionToZones[regionI].empty())
                {
                    if (overwrite)
                    {
                        regionNames[regionI] = polyMesh::defaultRegion;
                    }
                    else if (insidePoint)
                    {
                        regionNames[regionI] = "insidePoint";
                    }
                    else if (largestOnly)
                    {
                        regionNames[regionI] = "largestOnly";
                    }
                }
            }
        }
    }

    Info<< endl << "Number of regions:" << nCellRegions << nl << endl;


    // Write decomposition to file
    writeCellToRegion(mesh, cellRegion);



    // Sizes per region
    // ~~~~~~~~~~~~~~~~

    labelList regionSizes(nCellRegions, Zero);

    forAll(cellRegion, celli)
    {
        regionSizes[cellRegion[celli]]++;
    }
    forAll(regionSizes, regionI)
    {
        reduce(regionSizes[regionI], sumOp<label>());
    }

    Info<< "Region\tCells" << nl
        << "------\t-----" << endl;

    forAll(regionSizes, regionI)
    {
        Info<< regionI << "\t\t" << regionSizes[regionI] << nl;
    }
    Info<< endl;



    // Print region to zone
    Info<< "Region\tZone\tName" << nl
        << "------\t----\t----" << endl;
    forAll(regionToZones, regionI)
    {
        Info<< regionI << "\t\t" << flatOutput(regionToZones[regionI])
            << '\t'
            << regionNames[regionI] << nl;
    }
    Info<< endl;



    // Since we're going to mess with patches and zones make sure all
    // is synchronised
    mesh.boundaryMesh().checkParallelSync(true);
    mesh.faceZones().checkParallelSync(true);


    // Interfaces
    // ----------
    // per interface:
    // - the two regions on either side
    // - the name
    // - the (global) size
    edgeList interfaces;
    List<Pair<word>> interfaceNames;
    labelList interfaceSizes;
    // per face the interface
    labelList faceToInterface;

    getInterfaceSizes
    (
        mesh,
        useFaceZones,
        cellRegion,
        regionNames,

        interfaces,
        interfaceNames,
        interfaceSizes,
        faceToInterface
    );

    Info<< "Sizes of interfaces between regions:" << nl << nl
        << "Interface\tRegion\tRegion\tFaces" << nl
        << "---------\t------\t------\t-----" << endl;

    forAll(interfaces, interI)
    {
        const edge& e = interfaces[interI];

        Info<< interI
            << "\t\t\t" << e[0] << "\t\t" << e[1]
            << "\t\t" << interfaceSizes[interI] << nl;
    }
    Info<< endl;


    if (detectOnly)
    {
        Info<< "End\n" << endl;

        return 0;
    }


    // Read objects in time directory
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    IOobjectList objects(mesh, runTime.timeName());

    // Read vol fields.

    PtrList<volScalarField> vsFlds;
    ReadFields(mesh, objects, vsFlds);

    PtrList<volVectorField> vvFlds;
    ReadFields(mesh, objects, vvFlds);

    PtrList<volSphericalTensorField> vstFlds;
    ReadFields(mesh, objects, vstFlds);

    PtrList<volSymmTensorField> vsymtFlds;
    ReadFields(mesh, objects, vsymtFlds);

    PtrList<volTensorField> vtFlds;
    ReadFields(mesh, objects, vtFlds);

    // Read surface fields.

    PtrList<surfaceScalarField> ssFlds;
    ReadFields(mesh, objects, ssFlds);

    PtrList<surfaceVectorField> svFlds;
    ReadFields(mesh, objects, svFlds);

    PtrList<surfaceSphericalTensorField> sstFlds;
    ReadFields(mesh, objects, sstFlds);

    PtrList<surfaceSymmTensorField> ssymtFlds;
    ReadFields(mesh, objects, ssymtFlds);

    PtrList<surfaceTensorField> stFlds;
    ReadFields(mesh, objects, stFlds);

    Info<< endl;


    // Remove any demand-driven fields ('S', 'V' etc)
    mesh.clearOut();


    if (nCellRegions == 1)
    {
        Info<< "Only one region. Doing nothing." << endl;
    }
    else if (makeCellZones)
    {
        Info<< "Putting cells into cellZones instead of splitting mesh."
            << endl;

        // Check if region overlaps with existing zone. If so keep.

        for (label regionI = 0; regionI < nCellRegions; regionI++)
        {
            const labelList& zones = regionToZones[regionI];

            if (zones.size() == 1 && zones[0] != -1)
            {
                // Existing zone
                const label zoneI = zones[0];
                Info<< "    Region " << regionI << " : corresponds to existing"
                    << " cellZone "
                    << zoneI << ' ' << mesh.cellZones()[zoneI].name() << endl;
            }
            else
            {
                // Create new cellZone.
                const labelList regionCells(findIndices(cellRegion, regionI));

                const word zoneName = "region" + Foam::name(regionI);

                label zoneI = mesh.cellZones().findZoneID(zoneName);

                if (zoneI == -1)
                {
                    zoneI = mesh.cellZones().size();
                    mesh.cellZones().setSize(zoneI+1);
                    mesh.cellZones().set
                    (
                        zoneI,
                        new cellZone
                        (
                            zoneName,           //name
                            std::move(regionCells), //addressing
                            zoneI,              //index
                            mesh.cellZones()    //cellZoneMesh
                        )
                    );
                }
                else
                {
                    mesh.cellZones()[zoneI].clearAddressing();
                    mesh.cellZones()[zoneI] = regionCells;
                }
                Info<< "    Region " << regionI << " : created new cellZone "
                    << zoneI << ' ' << mesh.cellZones()[zoneI].name() << endl;
            }
        }
        mesh.cellZones().writeOpt(IOobject::AUTO_WRITE);

        if (!overwrite)
        {
            ++runTime;
            mesh.setInstance(runTime.timeName());
        }
        else
        {
            mesh.setInstance(oldInstance);
        }

        Info<< "Writing cellZones as new mesh to time " << runTime.timeName()
            << nl << endl;

        mesh.write();


        // Write cellSets for convenience
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        Info<< "Writing cellSets corresponding to cellZones." << nl << endl;

        for (const auto& cz : mesh.cellZones())
        {
            cellSet(mesh, cz.name(), cz).write();
        }
    }
    else
    {
        // Add patches for interfaces
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~

        // Add all possible patches. Empty ones get filtered later on.
        Info<< nl << "Adding patches" << nl << endl;

        labelList interfacePatches
        (
            addRegionPatches
            (
                mesh,
                regionNames,
                interfaces,
                interfaceNames
            )
        );


        if (!overwrite)
        {
            ++runTime;
        }


        // Create regions
        // ~~~~~~~~~~~~~~

        if (insidePoint)
        {
            const point insidePoint = args.get<point>("insidePoint");

            label regionI = -1;

            (void)mesh.tetBasePtIs();

            label celli = mesh.findCell(insidePoint);

            Info<< nl << "Found point " << insidePoint << " in cell " << celli
                << endl;

            if (celli != -1)
            {
                regionI = cellRegion[celli];
            }

            reduce(regionI, maxOp<label>());

            Info<< nl
                << "Subsetting region " << regionI
                << " containing point " << insidePoint << endl;

            if (regionI == -1)
            {
                FatalErrorInFunction
                    << "Point " << insidePoint
                    << " is not inside the mesh." << nl
                    << "Bounding box of the mesh:" << mesh.bounds()
                    << exit(FatalError);
            }

            createAndWriteRegion
            (
                mesh,
                cellRegion,
                regionNames,
                prefixRegion,
                faceToInterface,
                interfacePatches,
                regionI,
                (overwrite ? oldInstance : runTime.timeName())
            );
        }
        else if (largestOnly)
        {
            label regionI = findMax(regionSizes);

            Info<< nl
                << "Subsetting region " << regionI
                << " of size " << regionSizes[regionI]
                << " as named region " << regionNames[regionI] << endl;

            createAndWriteRegion
            (
                mesh,
                cellRegion,
                regionNames,
                prefixRegion,
                faceToInterface,
                interfacePatches,
                regionI,
                (overwrite ? oldInstance : runTime.timeName())
            );
        }
        else
        {
            // Split all
            for (label regionI = 0; regionI < nCellRegions; regionI++)
            {
                Info<< nl
                    << "Region " << regionI << nl
                    << "-------- " << endl;

                createAndWriteRegion
                (
                    mesh,
                    cellRegion,
                    regionNames,
                    prefixRegion,
                    faceToInterface,
                    interfacePatches,
                    regionI,
                    (overwrite ? oldInstance : runTime.timeName())
                );
            }
        }
    }

    Info<< "End\n" << endl;

    return 0;
}


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