Skip to content
Snippets Groups Projects
splitMeshRegions.C 52.2 KiB
Newer Older
  • Learn to ignore specific revisions
  •             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();
    }
    
    
    
    mattijs's avatar
    mattijs committed
    // 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'.
    
    EdgeMap<label> addRegionPatches
    (
        fvMesh& mesh,
    
        const labelList& cellRegion,
        const label nCellRegions,
    
    mattijs's avatar
    mattijs committed
        const edgeList& interfaces,
    
        const EdgeMap<label>& interfaceSizes,
        const wordList& regionNames
    )
    {
        // Check that all patches are present in same order.
        mesh.boundaryMesh().checkParallelSync(true);
    
        Info<< nl << "Adding patches" << nl << endl;
    
    
        EdgeMap<label> interfaceToPatch(nCellRegions);
    
    mattijs's avatar
    mattijs committed
        forAll(interfaces, interI)
    
    mattijs's avatar
    mattijs committed
            const edge& e = interfaces[interI];
    
    mattijs's avatar
    mattijs committed
            if (interfaceSizes[e] > 0)
            {
    
                const word inter1 = regionNames[e[0]] + "_to_" + regionNames[e[1]];
                const word inter2 = regionNames[e[1]] + "_to_" + regionNames[e[0]];
    
                directMappedWallPolyPatch patch1
    
                    inter1,
                    0,                  // overridden
                    0,                  // overridden
                    0,                  // overridden
                    regionNames[e[1]],  // sampleRegion
                    directMappedPatchBase::NEARESTPATCHFACE,
                    inter2,             // samplePatch
                    point::zero,        // offset
                    mesh.boundaryMesh()
    
    
                label patchI = addPatch(mesh, patch1);
    
                directMappedWallPolyPatch patch2
    
                    inter2,
                    0,
                    0,
                    0,
                    regionNames[e[0]],  // sampleRegion
                    directMappedPatchBase::NEARESTPATCHFACE,
                    inter1,
                    point::zero,        // offset
                    mesh.boundaryMesh()
    
    
                Info<< "For interface between region " << e[0]
                    << " and " << e[1] << " added patch " << patchI
                    << " " << mesh.boundaryMesh()[patchI].name()
                    << endl;
    
    
    mattijs's avatar
    mattijs committed
                interfaceToPatch.insert(e, patchI);
    
    //// Checks if regionI in cellRegion is subset of existing cellZone. Returns -1
    //// if no zone found, zone otherwise
    //label findCorrespondingSubZone
    //(
    //    const cellZoneMesh& cellZones,
    //    const labelList& existingZoneID,
    //    const labelList& cellRegion,
    //    const label regionI
    //)
    //{
    //    // Zone corresponding to region. No corresponding zone.
    //    label zoneI = labelMax;
    //
    //    labelList regionCells = findIndices(cellRegion, regionI);
    //
    //    if (regionCells.empty())
    //    {
    //        // My local portion is empty. Maps to any empty cellZone. Mark with
    //        // special value which can get overwritten by other processors.
    //        zoneI = -1;
    //    }
    //    else
    //    {
    //        // Get zone for first element.
    //        zoneI = existingZoneID[regionCells[0]];
    //
    //        if (zoneI == -1)
    //        {
    //            zoneI = labelMax;
    //        }
    //        else
    //        {
    //            // 1. All regionCells in zoneI?
    //            forAll(regionCells, i)
    //            {
    //                if (existingZoneID[regionCells[i]] != zoneI)
    //                {
    //                    zoneI = labelMax;
    //                    break;
    //                }
    //            }
    //        }
    //    }
    //
    //    // Determine same zone over all processors.
    //    reduce(zoneI, maxOp<label>());
    //
    //    if (zoneI == labelMax)
    //    {
    //        // Cells in region that are not in zoneI
    //        zoneI = -1;
    //    }
    //
    //    return zoneI;
    //}
    
    
    // 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, 0);
    
        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;
    
            // 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;
    
            // If one in error, all should be in error. Note that branch gets taken
            // on all procs.
    
            reduce(regionI, minOp<label>());
    
    
    
    //// Checks if cellZone has corresponding cellRegion.
    //label findCorrespondingRegion
    //(
    //    const cellZoneMesh& cellZones,
    //    const labelList& existingZoneID,    // per cell the (unique) zoneID
    
    //    const labelList& cellRegion,
    //    const label nCellRegions,
    
    //    const label zoneI
    //)
    //{
    //    // Region corresponding to zone. Start off with special value: no
    //    // corresponding region.
    //    label regionI = labelMax;
    //
    //    const cellZone& cz = cellZones[zoneI];
    //
    //    if (cz.empty())
    //    {
    //        // My local portion is empty. Maps to any empty cellZone. Mark with
    //        // special value which can get overwritten by other processors.
    //        regionI = -1;
    //    }
    //    else
    //    {
    //        regionI = cellRegion[cz[0]];
    //
    //        forAll(cz, i)
    //        {
    //            if (cellRegion[cz[i]] != regionI)
    //            {
    //                regionI = labelMax;
    //                break;
    //            }
    //        }
    //    }
    //
    //    // Determine same zone over all processors.
    //    reduce(regionI, maxOp<label>());
    //
    //
    //    // 2. All of region present?
    //
    //    if (regionI == labelMax)
    //    {
    //        regionI = -1;
    //    }
    //    else if (regionI != -1)
    //    {
    //        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;
    //}
    
    
    
    // Main program:
    
    int main(int argc, char *argv[])
    {
    
    #   include "addOverwriteOption.H"
    
        argList::addBoolOption("cellZones");
        argList::addBoolOption("cellZonesOnly");
        argList::addOption("blockedFaces", "faceSet");
        argList::addBoolOption("makeCellZones");
        argList::addBoolOption("largestOnly");
        argList::addOption("insidePoint", "point");
        argList::addBoolOption("detectOnly");
        argList::addBoolOption("sloppyCellZones");
    
    
    #   include "setRootCase.H"
    #   include "createTime.H"
    
    mattijs's avatar
    mattijs committed
        runTime.functionObjects().off();
    
    #   include "createMesh.H"
    
    mattijs's avatar
    mattijs committed
        const word oldInstance = mesh.pointsInstance();
    
        if (args.optionReadIfPresent("blockedFaces", blockedFacesName))
    
        {
            Info<< "Reading blocked internal faces from faceSet "
                << blockedFacesName << nl << endl;
        }
    
    
        const bool makeCellZones    = args.optionFound("makeCellZones");
        const bool largestOnly      = args.optionFound("largestOnly");
        const bool insidePoint      = args.optionFound("insidePoint");
        const bool useCellZones     = args.optionFound("cellZones");
        const bool useCellZonesOnly = args.optionFound("cellZonesOnly");
        const bool overwrite        = args.optionFound("overwrite");
        const bool detectOnly       = args.optionFound("detectOnly");
        const bool sloppyCellZones  = args.optionFound("sloppyCellZones");
    
    
        if (insidePoint && largestOnly)
        {
            FatalErrorIn(args.executable())
                << "You cannot specify both -largestOnly"
                << " (keep region with most cells)"
                << " and -insidePoint (keep region containing point)"
                << exit(FatalError);
        }
    
    
        const cellZoneMesh& cellZones = mesh.cellZones();
    
    
    
        // Collect zone per cell
        // ~~~~~~~~~~~~~~~~~~~~~
        // - non-unique zoning
        // - coupled zones
    
    
        // Existing zoneID
        labelList zoneID(mesh.nCells(), -1);
    
        forAll(cellZones, zoneI)
        {
    
            const cellZone& cz = cellZones[zoneI];
    
            forAll(cz, i)
            {
                label cellI = cz[i];
                if (zoneID[cellI] == -1)
                {
                    zoneID[cellI] = zoneI;
                }
                else
                {
                    FatalErrorIn(args.executable())
                        << "Cell " << cellI << " with cell centre "
                        << mesh.cellCentres()[cellI]
                        << " is multiple zones. This is not allowed." << endl
                        << "It is in zone " << cellZones[zoneID[cellI]].name()
                        << " and in zone " << cellZones[zoneI].name()
                        << exit(FatalError);
                }
            }
    
        }
    
        // Neighbour zoneID.
        labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
    
        forAll(neiZoneID, i)
        {
            neiZoneID[i] = zoneID[mesh.faceOwner()[i+mesh.nInternalFaces()]];
        }
        syncTools::swapBoundaryFaceList(mesh, neiZoneID, false);
    
    
        // 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);
    
            forAllConstIter(faceSet, blockedFaceSet, iter)
            {
                blockedFace[iter.key()] = true;
            }
        }
    
        // Imply from differing cellZones
        if (useCellZones)
        {
            blockedFace.setSize(mesh.nFaces(), false);
    
            for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
            {
                label own = mesh.faceOwner()[faceI];
                label nei = mesh.faceNeighbour()[faceI];
    
                if (zoneID[own] != zoneID[nei])
                {
                    blockedFace[faceI] = true;
                }
            }
    
    
            // Different cellZones on either side of processor patch.
    
            forAll(neiZoneID, i)
            {
                label faceI = i+mesh.nInternalFaces();
    
                if (zoneID[mesh.faceOwner()[faceI]] != neiZoneID[i])
                {
    
                    blockedFace[faceI] = true;
    
        // Determine per cell the region it belongs to
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
        // cellRegion is the labelList with the region per cell.
        labelList cellRegion;
        label nCellRegions = 0;
        if (useCellZonesOnly)
        {
            label unzonedCellI = findIndex(zoneID, -1);
            if (unzonedCellI != -1)
            {
                FatalErrorIn(args.executable())
                    << "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 = zoneID;
            nCellRegions = gMax(cellRegion)+1;
        }
        else
        {
            // Do a topological walk to determine regions
            regionSplit regions(mesh, blockedFace);
            nCellRegions = regions.nRegions();
            cellRegion.transfer(regions);
        }
    
        Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
    
    
    
        // 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,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE,
                    false
                ),
                mesh,
    
                dimensionedScalar("zero", dimless, 0),
                zeroGradientFvPatchScalarField::typeName
    
            );
            forAll(cellRegion, cellI)
            {
                cellToRegion[cellI] = cellRegion[cellI];
            }
            cellToRegion.write();
    
            Info<< "Writing region per cell as volScalarField to "
                << cellToRegion.objectPath() << nl << endl;
        }
    
    
    
        // Sizes per region
        // ~~~~~~~~~~~~~~~~
    
    
        labelList regionSizes(nCellRegions, 0);
    
    
        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' << regionSizes[regionI] << nl;
        }
        Info<< endl;
    
    
    
        // Sizes per cellzone
        // ~~~~~~~~~~~~~~~~~~
    
        labelList zoneSizes(cellZones.size(), 0);
        if (useCellZones || makeCellZones || sloppyCellZones)
        {
            List<wordList> zoneNames(Pstream::nProcs());
            zoneNames[Pstream::myProcNo()] = cellZones.names();
            Pstream::gatherList(zoneNames);
            Pstream::scatterList(zoneNames);
    
            forAll(zoneNames, procI)
            {
                if (zoneNames[procI] != zoneNames[0])
                {
                    FatalErrorIn(args.executable())
                        << "cellZones not synchronised across processors." << endl
                        << "Master has cellZones " << zoneNames[0] << endl
                        << "Processor " << procI
                        << " has cellZones " << zoneNames[procI]
                        << exit(FatalError);
                }
            }
    
            forAll(cellZones, zoneI)
            {
                zoneSizes[zoneI] = returnReduce
                (
                    cellZones[zoneI].size(),
                    sumOp<label>()
                );
            }
        }
    
    
    
        // Whether region corresponds to a cellzone
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
        // Region per zone
    
        labelList regionToZone(nCellRegions, -1);
    
        wordList regionNames(nCellRegions);
    
        // Zone to region
        labelList zoneToRegion(cellZones.size(), -1);
    
        if (sloppyCellZones)
    
            Info<< "Trying to match regions to existing cell zones;"
                << " region can be subset of cell zone." << nl << endl;
    
            forAll(cellZones, zoneI)
    
                label regionI = findCorrespondingRegion
                (
                    zoneID,
                    cellRegion,
    
                    zoneI,
                    label(0.5*zoneSizes[zoneI]) // minimum overlap
                );
    
                if (regionI != -1)
                {
                    Info<< "Sloppily matched region " << regionI
                        << " size " << regionSizes[regionI]
                        << " to zone " << zoneI << " size " << zoneSizes[zoneI]
                        << endl;
                    zoneToRegion[zoneI] = regionI;
                    regionToZone[regionI] = zoneI;
                    regionNames[regionI] = cellZones[zoneI].name();
                }
    
        }
        else
        {
            Info<< "Trying to match regions to existing cell zones." << nl << endl;
    
            forAll(cellZones, zoneI)
            {
                label regionI = findCorrespondingRegion
                (
                    zoneID,
                    cellRegion,
    
                    zoneI,
                    1               // minimum overlap
                );
    
                if (regionI != -1)
                {
                    zoneToRegion[zoneI] = regionI;
                    regionToZone[regionI] = zoneI;
                    regionNames[regionI] = cellZones[zoneI].name();
                }
            }
        }
        // Allocate region names for unmatched regions.
        forAll(regionToZone, regionI)
        {
            if (regionToZone[regionI] == -1)
    
            {
                regionNames[regionI] = "domain" + Foam::name(regionI);
            }
    
        // Print region to zone
        Info<< "Region\tZone\tName" << nl
            << "------\t----\t----" << endl;
        forAll(regionToZone, regionI)
        {
    
            Info<< regionI << '\t' << regionToZone[regionI] << '\t'
                << regionNames[regionI] << nl;
        }
        Info<< endl;
    
    
        //// Print zone to region
        //Info<< "Zone\tName\tRegion" << nl
        //    << "----\t----\t------" << endl;
        //forAll(zoneToRegion, zoneI)
        //{
        //    Info<< zoneI << '\t' << cellZones[zoneI].name() << '\t'
        //        << zoneToRegion[zoneI] << nl;
        //}
        //Info<< endl;
    
    
    
    
        // Since we're going to mess with patches make sure all non-processor ones
        // are on all processors.
        mesh.boundaryMesh().checkParallelSync(true);
    
    
        // Sizes of interface between regions. From pair of regions to number of
        // faces.
    
    mattijs's avatar
    mattijs committed
        edgeList interfaces;
        EdgeMap<label> interfaceSizes;
        getInterfaceSizes
    
    mattijs's avatar
    mattijs committed
            mesh,
            cellRegion,
            true,      // sum in parallel?
    
            interfaces,
            interfaceSizes
    
        Info<< "Sizes inbetween regions:" << nl << nl
            << "Region\tRegion\tFaces" << nl
    
            << "------\t------\t-----" << endl;
    
    
    mattijs's avatar
    mattijs committed
        forAll(interfaces, interI)
    
    mattijs's avatar
    mattijs committed
            const edge& e = interfaces[interI];
    
    mattijs's avatar
    mattijs committed
            Info<< e[0] << '\t' << e[1] << '\t' << interfaceSizes[e] << nl;
    
    
    
        // 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++)
    
            {
                label zoneI = regionToZone[regionI];
    
                if (zoneI != -1)
                {
                    Info<< "    Region " << regionI << " : corresponds to existing"
                        << " cellZone "
                        << zoneI << ' ' << cellZones[zoneI].name() << endl;
                }
                else
                {
                    // Create new cellZone.
                    labelList regionCells = findIndices(cellRegion, regionI);
    
                    word zoneName = "region" + Foam::name(regionI);
    
                    zoneI = cellZones.findZoneID(zoneName);
    
                    if (zoneI == -1)
                    {
                        zoneI = cellZones.size();
                        mesh.cellZones().setSize(zoneI+1);
                        mesh.cellZones().set
                        (
                            zoneI,
                            new cellZone
                            (
                                zoneName,           //name
                                regionCells,        //addressing
                                zoneI,              //index
                                cellZones           //cellZoneMesh
                            )
                        );
                    }
                    else
                    {
                        mesh.cellZones()[zoneI].clearAddressing();
                        mesh.cellZones()[zoneI] = regionCells;
                    }
                    Info<< "    Region " << regionI << " : created new cellZone "
                        << zoneI << ' ' << cellZones[zoneI].name() << endl;
                }
            }
            mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
    
    Andrew Heather's avatar
    Andrew Heather committed
    
            if (!overwrite)
            {
                runTime++;
    
    mattijs's avatar
    mattijs committed
                mesh.setInstance(runTime.timeName());
            }
            else
            {
                mesh.setInstance(oldInstance);
    
    Andrew Heather's avatar
    Andrew Heather committed
            }
    
    
            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;
    
            forAll(cellZones, zoneI)
            {
                const cellZone& cz = cellZones[zoneI];
    
                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;
    
            EdgeMap<label> interfaceToPatch
            (
                addRegionPatches
                (
                    mesh,
                    cellRegion,
    
    mattijs's avatar
    mattijs committed
                    interfaces,
    
    Andrew Heather's avatar
    Andrew Heather committed
            if (!overwrite)
            {
                runTime++;
            }
    
    
    
            // Create regions
            // ~~~~~~~~~~~~~~
    
            if (insidePoint)
            {
    
                const point insidePoint = args.optionRead<point>("insidePoint");
    
    
                label regionI = -1;
    
                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)
                {
                    FatalErrorIn(args.executable())
                        << "Point " << insidePoint
                        << " is not inside the mesh." << nl
                        << "Bounding box of the mesh:" << mesh.bounds()
                        << exit(FatalError);
                }
    
                createAndWriteRegion
                (
                    mesh,
                    cellRegion,
                    regionNames,
                    interfaceToPatch,
    
    mattijs's avatar
    mattijs committed
                    regionI,
                    (overwrite ? oldInstance : runTime.timeName())
    
                );
            }
            else if (largestOnly)
            {
                label regionI = findMax(regionSizes);
    
                Info<< nl
                    << "Subsetting region " << regionI
                    << " of size " << regionSizes[regionI] << endl;
    
                createAndWriteRegion
                (
                    mesh,
                    cellRegion,
                    regionNames,
                    interfaceToPatch,
    
    mattijs's avatar
    mattijs committed
                    regionI,
                    (overwrite ? oldInstance : runTime.timeName())
    
                for (label regionI = 0; regionI < nCellRegions; regionI++)
    
                {
                    Info<< nl
                        << "Region " << regionI << nl
                        << "-------- " << endl;
    
                    createAndWriteRegion
                    (
                        mesh,
                        cellRegion,
                        regionNames,
                        interfaceToPatch,
    
    mattijs's avatar
    mattijs committed
                        regionI,
                        (overwrite ? oldInstance : runTime.timeName())
    
                    );
                }
            }
        }
    
        Info<< "End\n" << endl;
    
        return 0;
    }
    
    
    // ************************************************************************* //