Skip to content
Snippets Groups Projects
splitMeshRegions.C 60.3 KiB
Newer Older
  • Learn to ignore specific revisions
  • /*---------------------------------------------------------------------------*\
      =========                 |
      \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       \\    /   O peration     |
    
        \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
    
         \\/     M anipulation  |
    -------------------------------------------------------------------------------
    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/>.
    
    mattijs's avatar
    mattijs committed
    Application
        splitMeshRegions
    
        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 inbetween 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)
    
    
        - 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.
    
        - useCellZonesOnly 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.
    
    
        - 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
    
        - 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!)
    
    
    mattijs's avatar
    mattijs committed
        - 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
    
    \*---------------------------------------------------------------------------*/
    
    #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 "EdgeMap.H"
    #include "syncTools.H"
    #include "ReadFields.H"
    
    #include "mappedWallPolyPatch.H"
    
    #include "zeroGradientFvPatchFields.H"
    
    
    using namespace Foam;
    
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
    template<class GeoField>
    void addPatchFields(fvMesh& mesh, const word& patchFieldType)
    {
        HashTable<const GeoField*> flds
        (
            mesh.objectRegistry::lookupClass<GeoField>()
        );
    
    
        forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
    
        {
            const GeoField& fld = *iter();
    
            typename GeoField::GeometricBoundaryField& bfld =
                const_cast<typename GeoField::GeometricBoundaryField&>
                (
                    fld.boundaryField()
                );
    
            label sz = bfld.size();
            bfld.setSize(sz+1);
            bfld.set
            (
                sz,
                GeoField::PatchFieldType::New
                (
                    patchFieldType,
                    mesh.boundary()[sz],
                    fld.dimensionedInternalField()
                )
            );
        }
    }
    
    
    // Remove last patch field
    template<class GeoField>
    void trimPatchFields(fvMesh& mesh, const label nPatches)
    {
        HashTable<const GeoField*> flds
        (
            mesh.objectRegistry::lookupClass<GeoField>()
        );
    
    
        forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
    
        {
            const GeoField& fld = *iter();
    
            const_cast<typename GeoField::GeometricBoundaryField&>
            (
                fld.boundaryField()
            ).setSize(nPatches);
        }
    }
    
    
    // Reorder patch field
    template<class GeoField>
    void reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
    {
        HashTable<const GeoField*> flds
        (
            mesh.objectRegistry::lookupClass<GeoField>()
        );
    
    
        forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
    
        {
            const GeoField& fld = *iter();
    
            typename GeoField::GeometricBoundaryField& bfld =
                const_cast<typename GeoField::GeometricBoundaryField&>
                (
                    fld.boundaryField()
                );
    
            bfld.reorder(oldToNew);
        }
    }
    
    
    // Adds patch if not yet there. Returns patchID.
    
    label addPatch(fvMesh& mesh, const polyPatch& patch)
    
    {
        polyBoundaryMesh& polyPatches =
            const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
    
    
        label patchI = polyPatches.findPatchID(patch.name());
    
            if (polyPatches[patchI].type() == patch.type())
    
            {
                // Already there
                return patchI;
            }
    
                FatalErrorIn("addPatch(fvMesh&, const polyPatch*)")
                    << "Already have patch " << patch.name()
                    << " but of type " << patch.type()
    
        }
    
    
        label insertPatchI = polyPatches.size();
        label startFaceI = mesh.nFaces();
    
        forAll(polyPatches, patchI)
        {
            const polyPatch& pp = polyPatches[patchI];
    
            if (isA<processorPolyPatch>(pp))
            {
                insertPatchI = patchI;
                startFaceI = pp.start();
                break;
            }
        }
    
    
        // Below is all quite a hack. Feel free to change once there is a better
        // mechanism to insert and reorder patches.
    
        // Clear local fields and e.g. polyMesh parallelInfo.
        mesh.clearOut();
    
        label sz = polyPatches.size();
    
        fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
    
        // Add polyPatch at the end
        polyPatches.setSize(sz+1);
        polyPatches.set
        (
            sz,
    
                polyPatches,
                insertPatchI,   //index
                0,              //size
                startFaceI      //start
    
            )
        );
        fvPatches.setSize(sz+1);
        fvPatches.set
        (
            sz,
            fvPatch::New
            (
                polyPatches[sz],  // point to newly added polyPatch
                mesh.boundary()
            )
        );
    
        addPatchFields<volScalarField>
        (
            mesh,
            calculatedFvPatchField<scalar>::typeName
        );
        addPatchFields<volVectorField>
        (
            mesh,
            calculatedFvPatchField<vector>::typeName
        );
        addPatchFields<volSphericalTensorField>
        (
            mesh,
            calculatedFvPatchField<sphericalTensor>::typeName
        );
        addPatchFields<volSymmTensorField>
        (
            mesh,
            calculatedFvPatchField<symmTensor>::typeName
        );
        addPatchFields<volTensorField>
        (
            mesh,
            calculatedFvPatchField<tensor>::typeName
        );
    
        // Surface fields
    
        addPatchFields<surfaceScalarField>
        (
            mesh,
            calculatedFvPatchField<scalar>::typeName
        );
        addPatchFields<surfaceVectorField>
        (
            mesh,
            calculatedFvPatchField<vector>::typeName
        );
        addPatchFields<surfaceSphericalTensorField>
        (
            mesh,
            calculatedFvPatchField<sphericalTensor>::typeName
        );
        addPatchFields<surfaceSymmTensorField>
        (
            mesh,
            calculatedFvPatchField<symmTensor>::typeName
        );
        addPatchFields<surfaceTensorField>
        (
            mesh,
            calculatedFvPatchField<tensor>::typeName
        );
    
        // Create reordering list
        // patches before insert position stay as is
        labelList oldToNew(sz+1);
        for (label i = 0; i < insertPatchI; i++)
        {
            oldToNew[i] = i;
        }
        // patches after insert position move one up
        for (label i = insertPatchI; i < sz; i++)
        {
            oldToNew[i] = i+1;
        }
        // appended patch gets moved to insert position
        oldToNew[sz] = insertPatchI;
    
        // Shuffle into place
        polyPatches.reorder(oldToNew);
        fvPatches.reorder(oldToNew);
    
        reorderPatchFields<volScalarField>(mesh, oldToNew);
        reorderPatchFields<volVectorField>(mesh, oldToNew);
        reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
        reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
        reorderPatchFields<volTensorField>(mesh, oldToNew);
        reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
        reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
        reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
        reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
        reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
    
        return insertPatchI;
    }
    
    
    // Reorder and delete patches.
    void reorderPatches
    (
        fvMesh& mesh,
        const labelList& oldToNew,
        const label nNewPatches
    )
    {
        polyBoundaryMesh& polyPatches =
            const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
        fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
    
        // Shuffle into place
        polyPatches.reorder(oldToNew);
        fvPatches.reorder(oldToNew);
    
        reorderPatchFields<volScalarField>(mesh, oldToNew);
        reorderPatchFields<volVectorField>(mesh, oldToNew);
        reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
        reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
        reorderPatchFields<volTensorField>(mesh, oldToNew);
        reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
        reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
        reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
        reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
        reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
    
        // Remove last.
        polyPatches.setSize(nNewPatches);
        fvPatches.setSize(nNewPatches);
        trimPatchFields<volScalarField>(mesh, nNewPatches);
        trimPatchFields<volVectorField>(mesh, nNewPatches);
        trimPatchFields<volSphericalTensorField>(mesh, nNewPatches);
        trimPatchFields<volSymmTensorField>(mesh, nNewPatches);
        trimPatchFields<volTensorField>(mesh, nNewPatches);
        trimPatchFields<surfaceScalarField>(mesh, nNewPatches);
        trimPatchFields<surfaceVectorField>(mesh, nNewPatches);
        trimPatchFields<surfaceSphericalTensorField>(mesh, nNewPatches);
        trimPatchFields<surfaceSymmTensorField>(mesh, nNewPatches);
        trimPatchFields<surfaceTensorField>(mesh, nNewPatches);
    }
    
    
    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>()
        );
        forAllConstIter(typename HashTable<const GeoField*>, fields, iter)
        {
            const GeoField& fld = *iter();
    
            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().boundaryField()[patchI] ==
                        pTraits<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& faceMap,
        const labelHashSet& addedPatches
    
    )
    {
        const labelList patchMap(identity(mesh.boundaryMesh().size()));
    
        HashTable<const GeoField*> fields
        (
            mesh.objectRegistry::lookupClass<GeoField>()
        );
        forAllConstIter(typename HashTable<const GeoField*>, fields, iter)
        {
            const GeoField& fld = *iter();
    
            Info<< "Mapping field " << fld.name() << endl;
    
            tmp<GeoField> tSubFld
            (
                fvMeshSubset::interpolate
                (
                    fld,
                    subMesh,
                    patchMap,
                    faceMap
                )
            );
    
            // Hack: set value to 0 for introduced patches (since don't
            //       get initialised.
            forAll(tSubFld().boundaryField(), patchI)
            {
    
                if (addedPatches.found(patchI))
    
                {
                    tSubFld().boundaryField()[patchI] ==
                        pTraits<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)
        );
    
        EdgeMap<Map<label> >::iterator iter = regionsToSize.find
        (
            interface
        );
    
        if (iter != regionsToSize.end())
        {
            // Check if zone present
            Map<label>::iterator zoneFnd = iter().find(zoneID);
            if (zoneFnd != iter().end())
            {
                // Found zone. Increment count.
                zoneFnd()++;
            }
            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.
    
    mattijs's avatar
    mattijs committed
    // Returns interfaces as straight list for looping in identical order.
    void getInterfaceSizes
    
    (
        const polyMesh& mesh,
    
        const labelList& cellRegion,
    
    mattijs's avatar
    mattijs committed
    
        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)
            {
    
                    mesh,
                    (useFaceZones ? mesh.faceZones().whichZone(faceI) : -1),
                    ownRegion,
                    neiRegion,
                    regionsToSize
    
        // Boundary faces
        // ~~~~~~~~~~~~~~
    
        // Neighbour cellRegion.
        labelList coupledRegion(mesh.nFaces()-mesh.nInternalFaces());
    
        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)
            {
    
                    mesh,
                    (useFaceZones ? mesh.faceZones().whichZone(faceI) : -1),
                    ownRegion,
                    neiRegion,
                    regionsToSize
    
        {
            if (Pstream::master())
            {
                // Receive and add to my sizes
                for
                (
                    int slave=Pstream::firstSlave();
                    slave<=Pstream::lastSlave();
                    slave++
                )
                {
                    IPstream fromSlave(Pstream::blocking, slave);
    
    
                    EdgeMap<Map<label> > slaveSizes(fromSlave);
    
                    forAllConstIter(EdgeMap<Map<label> >, slaveSizes, slaveIter)
    
                        EdgeMap<Map<label> >::iterator masterIter =
                            regionsToSize.find(slaveIter.key());
    
                        if (masterIter != regionsToSize.end())
    
                            // Same inter-region
                            const Map<label>& slaveInfo = slaveIter();
                            Map<label>& masterInfo = masterIter();
    
                            forAllConstIter(Map<label>, slaveInfo, iter)
                            {
                                label zoneID = iter.key();
                                label slaveSize = iter();
    
                                Map<label>::iterator zoneFnd = masterInfo.find
                                (
                                    zoneID
                                );
                                if (zoneFnd != masterInfo.end())
                                {
                                    zoneFnd() += slaveSize;
                                }
                                else
                                {
                                    masterInfo.insert(zoneID, slaveSize);
                                }
                            }
    
                            regionsToSize.insert(slaveIter.key(), slaveIter());
    
    mattijs's avatar
    mattijs committed
                    OPstream toMaster(Pstream::blocking, Pstream::masterNo());
    
    andy's avatar
    andy committed
        // Rework
    
    
        Pstream::scatter(regionsToSize);
    
    
    
        // Now we have the global sizes of all inter-regions.
        // Invert this on master and distribute.
        label nInterfaces = 0;
        forAllConstIter(EdgeMap<Map<label> >, regionsToSize, iter)
        {
            const Map<label>& info = iter();
            nInterfaces += info.size();
        }
    
        interfaces.setSize(nInterfaces);
        interfaceNames.setSize(nInterfaces);
        interfaceSizes.setSize(nInterfaces);
        EdgeMap<Map<label> > regionsToInterface(nInterfaces);
    
        nInterfaces = 0;
        forAllConstIter(EdgeMap<Map<label> >, regionsToSize, iter)
        {
            const edge& e = iter.key();
            const word& name0 = regionNames[e[0]];
            const word& name1 = regionNames[e[1]];
    
            const Map<label>& info = iter();
            forAllConstIter(Map<label>, 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();
    
                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)
    
    mattijs's avatar
    mattijs committed
        {
    
            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];
            }
    
    mattijs's avatar
    mattijs committed
        }
    
    }
    
    
    // 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*
        {
            IOobject io
            (
                "fvSchemes",
                mesh.time().system(),
                regionName,
    
    mattijs's avatar
    mattijs committed
                mesh,
    
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            );
    
            Info<< "Testing:" << io.objectPath() << endl;
    
            if (!io.headerOk())
    
            // if (!exists(io.objectPath()))
    
            {
                Info<< "Writing dummy " << regionName/io.name() << endl;
                dictionary dummyDict;
                dictionary divDict;
                dummyDict.add("divSchemes", divDict);
                dictionary gradDict;
                dummyDict.add("gradSchemes", gradDict);
                dictionary laplDict;
                dummyDict.add("laplacianSchemes", laplDict);
    
                IOdictionary(io, dummyDict).regIOobject::write();
            }
        }
        {
            IOobject io
            (
                "fvSolution",
                mesh.time().system(),
                regionName,
    
    mattijs's avatar
    mattijs committed
                mesh,
    
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            );
    
            if (!io.headerOk())
            //if (!exists(io.objectPath()))
            {
                Info<< "Writing dummy " << regionName/io.name() << endl;
                dictionary dummyDict;
                IOdictionary(io, dummyDict).regIOobject::write();
            }
        }
    
    
    
        // Neighbour cellRegion.
        labelList coupledRegion(mesh.nFaces()-mesh.nInternalFaces());
    
        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
            {
                FatalErrorIn("createRegionMesh(..)")
                    << "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;
    
    Andrew Heather's avatar
    Andrew Heather committed
            }
    
        }
    
        // Remove faces
        cellRemover.setRefinement
        (
            cellsToRemove,
            exposedFaces,
            exposedPatchIDs,
            meshMod
        );
    
        autoPtr<mapPolyMesh> map = meshMod.makeMesh
        (
            newMesh,
            IOobject
            (
                regionName,
                mesh.time().timeName(),
                mesh.time(),
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh
        );
    
        return map;
    }
    
    
    void createAndWriteRegion
    (
        const fvMesh& mesh,
    
        const labelList& cellRegion,
    
        const wordList& regionNames,
    
        const labelList& faceToInterface,
        const labelList& interfacePatches,
    
    mattijs's avatar
    mattijs committed
        const label regionI,
        const word& newMeshInstance
    
    )
    {
        Info<< "Creating mesh for region " << regionI
            << ' ' << regionNames[regionI] << endl;
    
        autoPtr<fvMesh> newMesh;
        autoPtr<mapPolyMesh> map = createRegionMesh
        (
            mesh,