ENH: fvMeshSubset improvements (issue #951)

- what was previously termed 'setLargeCellSubset()' is now simply
  'setCellSubset()' and supports memory efficient interfaces.

  The new parameter ordering avoids ambiguities caused by default
  parameters.

  Old parameter order:

      setLargeCellSubset
      (
          const labelList& region,
          const label currentRegion,
          const label patchID = -1,
          const bool syncCouples = true
      );

  New parameter order:

      setCellSubset
      (
          const label regioni,
          const labelUList& regions,
          const label patchID = -1,
          const bool syncCouples = true
      );

   And without ambiguity:

      setCellSubset
      (
          const labelUList& selectedCells,
          const label patchID = -1,
          const bool syncCouples = true
      );

- support bitSet directly for specifying the selectedCells for
  memory efficiency and ease of use.

- Additional constructors to perform setCellSubset() immediately,
  which simplifies coding.

  For example,

      meshParts.set
      (
          zonei,
          new fvMeshSubset(mesh, selectedCells)
      );

  Or even

      return autoPtr<fvMeshSubset>::New(mesh, selectedCells);

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H

@@ -25,10 +25,15 @@ labelList receptorNeigCell(mesh.nInternalFaces(), -1);
     PtrList<fvMeshSubset> meshParts(nZones);
     labelList nCellsPerZone(nZones, 0);
 
-    forAll(nCellsPerZone, zoneI)
+    // A mesh subset for each zone
+    forAll(meshParts, zonei)
     {
-        meshParts.set(zoneI, new fvMeshSubset(mesh));
-        meshParts[zoneI].setLargeCellSubset(zoneID, zoneI);
+        meshParts.set
+        (
+            zonei,
+            // Select cells where the zoneID == zonei
+            new fvMeshSubset(mesh, zonei, zoneID)
+        );
     }
 
     for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)

applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C

@@ -51,6 +51,7 @@ Description
 #include "fvMeshSubset.H"
 #include "argList.H"
 #include "cellSet.H"
+#include "BitOps.H"
 #include "IOobjectList.H"
 #include "volFields.H"
 #include "mapPolyMesh.H"
@@ -774,19 +775,19 @@ int main(int argc, char *argv[])
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     //
 
-    // Create mesh subsetting engine
-    fvMeshSubset subsetter(mesh);
-
-    {
-
-        cellSet blockedCells(mesh, blockedSetName);
-
-        // invert
-        blockedCells.invert(mesh.nCells());
-
-        // Create subsetted mesh.
-        subsetter.setLargeCellSubset(blockedCells, defaultPatchi, true);
-    }
+    // Mesh subsetting engine
+    fvMeshSubset subsetter
+    (
+        mesh,
+        BitSetOps::create
+        (
+            mesh.nCells(),
+            cellSet(mesh, blockedSetName), // Blocked cells as labelHashSet
+            false  // on=false: invert logic => retain the unblocked cells
+        ),
+        defaultPatchi,
+        true
+    );
 
 
     // Subset wantedPatch. Note that might also include boundary faces

applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C

@@ -570,17 +570,16 @@ labelList regionRenumber
 
     label celli = 0;
 
-    forAll(regionToCells, regionI)
+    forAll(regionToCells, regioni)
     {
-        Info<< "    region " << regionI << " starts at " << celli << endl;
+        Info<< "    region " << regioni << " starts at " << celli << endl;
 
         // Make sure no parallel comms
-        bool oldParRun = UPstream::parRun();
+        const bool oldParRun = UPstream::parRun();
         UPstream::parRun() = false;
 
         // Per region do a reordering.
-        fvMeshSubset subsetter(mesh);
-        subsetter.setLargeCellSubset(cellToRegion, regionI);
+        fvMeshSubset subsetter(mesh, regioni, cellToRegion);
 
         const fvMesh& subMesh = subsetter.subMesh();
 

applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,6 +45,7 @@ Description
 #include "volFields.H"
 #include "topoDistanceData.H"
 #include "FaceCellWave.H"
+#include "BitOps.H"
 #include "cellSet.H"
 #include "faceSet.H"
 #include "pointSet.H"
@@ -311,7 +312,7 @@ void subsetTopoSets
         label nSet = 0;
         forAll(map, i)
         {
-            if (isSet[map[i]])
+            if (isSet.test(map[i]))
             {
                 ++nSet;
             }
@@ -325,7 +326,7 @@ void subsetTopoSets
         TopoSet& subSet = subSets[i];
         forAll(map, i)
         {
-            if (isSet[map[i]])
+            if (isSet.test(map[i]))
             {
                 subSet.insert(i);
             }
@@ -371,7 +372,7 @@ int main(int argc, char *argv[])
 
     #include "createNamedMesh.H"
 
-    const word setName = args[1];
+    const word selectionName = args[1];
 
     word meshInstance = mesh.pointsInstance();
     word fieldsInstance = runTime.timeName();
@@ -389,73 +390,91 @@ int main(int argc, char *argv[])
     }
 
 
-    Info<< "Reading cell set from " << setName << nl << endl;
+    Info<< "Reading cell set from " << selectionName << nl << endl;
 
-    // Create mesh subsetting engine
-    fvMeshSubset subsetter(mesh);
 
-    labelList exposedPatchIDs;
+    // Default exposed patch id
+    labelList exposedPatchIDs(one(), -1);
 
-    if (args.found("patch"))
+    if (args.found("patches"))
     {
-        const word patchName = args["patch"];
+        const wordRes patchNames(args.readList<wordRe>("patches"));
 
-        exposedPatchIDs = { mesh.boundaryMesh().findPatchID(patchName) };
+        exposedPatchIDs = mesh.boundaryMesh().patchSet(patchNames).sortedToc();
 
-        if (exposedPatchIDs[0] == -1)
+        Info<< "Adding exposed internal faces to nearest of patches "
+            << patchNames << nl << endl;
+
+        if (exposedPatchIDs.empty())
         {
             FatalErrorInFunction
-                << nl << "Valid patches are " << mesh.boundaryMesh().names()
+                << nl << "No patches matched. Patches: "
+                << mesh.boundaryMesh().names()
                 << exit(FatalError);
         }
 
-        Info<< "Adding exposed internal faces to patch " << patchName
-            << nl << endl;
     }
-    else if (args.found("patches"))
+    else if (args.found("patch"))
     {
-        const wordRes patchNames(args.readList<wordRe>("patches"));
+        const word patchName = args["patch"];
 
-        exposedPatchIDs = mesh.boundaryMesh().patchSet(patchNames).sortedToc();
+        exposedPatchIDs.first() = mesh.boundaryMesh().findPatchID(patchName);
 
-        Info<< "Adding exposed internal faces to nearest of patches "
-            << patchNames << nl << endl;
+        Info<< "Adding exposed internal faces to patch " << patchName
+            << nl << endl;
+
+        if (exposedPatchIDs.first() == -1)
+        {
+            FatalErrorInFunction
+                << nl << "No such patch. Patches: "
+                << mesh.boundaryMesh().names()
+                << exit(FatalError);
+        }
     }
     else
     {
         Info<< "Adding exposed internal faces to a patch called"
             << " \"oldInternalFaces\" (created if necessary)" << endl
             << endl;
-        exposedPatchIDs = { -1 };
     }
 
 
-    cellSet currentSet(mesh, setName);
+    // Mesh subsetting engine
+    fvMeshSubset subsetter(mesh);
+
+    cellSet currentSet(mesh, selectionName);
 
-    if (exposedPatchIDs.size() == 1)
-    {
-        subsetter.setLargeCellSubset(currentSet, exposedPatchIDs[0], true);
-    }
-    else
     {
-        // Find per face the nearest patch
-        labelList nearestExposedPatch(nearestPatch(mesh, exposedPatchIDs));
+        bitSet selectedCells = BitSetOps::create(mesh.nCells(), currentSet);
 
-        labelList region(mesh.nCells(), 0);
-        forAllConstIter(cellSet, currentSet, iter)
+        if (exposedPatchIDs.size() == 1)
         {
-            region[iter.key()] = 1;
+            // Single patch for exposed faces
+            subsetter.setCellSubset
+            (
+                selectedCells,
+                exposedPatchIDs.first(),
+                true
+            );
         }
+        else
+        {
+            // The nearest patch per face
+            labelList nearestExposedPatch(nearestPatch(mesh, exposedPatchIDs));
 
-        labelList exposedFaces(subsetter.getExposedFaces(region, 1, true));
-        subsetter.setLargeCellSubset
-        (
-            region,
-            1,
-            exposedFaces,
-            labelUIndList(nearestExposedPatch, exposedFaces)(),
-            true
-        );
+            labelList exposedFaces
+            (
+                subsetter.getExposedFaces(selectedCells, true)
+            );
+
+            subsetter.setCellSubset
+            (
+                selectedCells,
+                exposedFaces,
+                labelUIndList(nearestExposedPatch, exposedFaces)(),
+                true
+            );
+        }
     }
 
 

applications/utilities/parallelProcessing/redistributePar/redistributePar.C

@@ -872,30 +872,23 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
             // Find last non-processor patch.
             const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-            label nonProcI = -1;
+            const label nonProcI = (patches.nNonProcessor() - 1);
 
-            forAll(patches, patchI)
-            {
-                if (isA<processorPolyPatch>(patches[patchI]))
-                {
-                    break;
-                }
-                nonProcI++;
-            }
-
-            if (nonProcI == -1)
+            if (nonProcI < 0)
             {
                 FatalErrorInFunction
                     << "Cannot find non-processor patch on processor "
-                    << Pstream::myProcNo() << endl
+                    << Pstream::myProcNo() << nl
                     << " Current patches:" << patches.names()
                     << abort(FatalError);
             }
 
-            // Subset 0 cells, no parallel comms. This is used to create
-            // zero-sized fields.
-            subsetterPtr.reset(new fvMeshSubset(mesh));
-            subsetterPtr().setLargeCellSubset(labelHashSet(0), nonProcI, false);
+            // Subset 0 cells, no parallel comms.
+            // This is used to create zero-sized fields.
+            subsetterPtr.reset
+            (
+                new fvMeshSubset(mesh, bitSet(), nonProcI, false)
+            );
         }
 
 

applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamMesh.C

@@ -423,8 +423,7 @@ void Foam::vtkPVFoam::convertMeshCellZones()
 
         if (!vtkgeom)
         {
-            fvMeshSubset subsetter(mesh);
-            subsetter.setLargeCellSubset(zMesh[zoneId]);
+            fvMeshSubset subsetter(mesh, zMesh[zoneId]);
 
             vtkgeom = vtuData.subset(subsetter, this->decomposePoly_);
         }
@@ -490,8 +489,7 @@ void Foam::vtkPVFoam::convertMeshCellSets()
 
         if (!vtkgeom)
         {
-            fvMeshSubset subsetter(mesh);
-            subsetter.setLargeCellSubset(cellSet(mesh, partName));
+            fvMeshSubset subsetter(mesh, cellSet(mesh, partName));
 
             vtkgeom = vtuData.subset(subsetter, this->decomposePoly_);
         }

src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C

@@ -2022,14 +2022,12 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
             //OPstream str(Pstream::commsTypes::blocking, recvProc);
             UOPstream str(recvProc, pBufs);
 
-            // Mesh subsetting engine
-            fvMeshSubset subsetter(mesh_);
-
-            // Subset the cells of the current domain.
-            subsetter.setLargeCellSubset
+            // Mesh subsetting engine - subset the cells of the current domain.
+            fvMeshSubset subsetter
             (
-                distribution,
+                mesh_,
                 recvProc,
+                distribution,
                 oldInternalPatchi,  // oldInternalFaces patch
                 false               // no parallel sync
             );

src/dynamicMesh/fvMeshSubset/fvMeshSubset.H

@@ -25,18 +25,15 @@ Class
     Foam::fvMeshSubset
 
 Description
-    Post-processing mesh subset tool.  Given the original mesh and the
-    list of selected cells, it creates the mesh consisting only of the
-    desired cells, with the mapping list for points, faces, and cells.
+    Given the original mesh and the list of selected cells, it creates the
+    mesh consisting only of the desired cells, with the mapping list for
+    points, faces, and cells.
 
     Puts all exposed internal faces into either
     - a user supplied patch
     - a newly created patch "oldInternalFaces"
 
-    - setLargeCellSubset is for largish subsets (>10% of mesh).
-      Uses labelLists instead.
-
-    - setLargeCellSubset does coupled patch subsetting as well. If it detects
+    - setCellSubset does coupled patch subsetting as well. If it detects
       a face on a coupled patch 'losing' its neighbour it will move the
       face into the oldInternalFaces patch.
 
@@ -58,6 +55,7 @@ SourceFiles
 #include "fvMesh.H"
 #include "pointMesh.H"
 #include "GeometricField.H"
+#include "bitSet.H"
 #include "HashSet.H"
 #include "surfaceMesh.H"
 
@@ -77,9 +75,12 @@ class fvMeshSubset
         //- Mesh to subset from
         const fvMesh& baseMesh_;
 
-        //- Subset mesh pointer
+        //- Demand-driven subset mesh pointer
         autoPtr<fvMesh> fvMeshSubsetPtr_;
 
+        //- Optional face mapping array with flip encoded (-1/+1)
+        mutable autoPtr<labelList> faceFlipMapPtr_;
+
         //- Point mapping array
         labelList pointMap_;
 
@@ -92,21 +93,14 @@ class fvMeshSubset
         //- Patch mapping array
         labelList patchMap_;
 
-        //- Optional face mapping array with flip encoded
-        mutable autoPtr<labelList> faceFlipMapPtr_;
-
 
     // Private Member Functions
 
         //- Check if subset has been performed
         bool checkCellSubset() const;
 
-        //- Mark points (with 0) in labelList
-        static void markPoints
-        (
-            const labelUList& curPoints,
-            labelList& pointMap
-        );
+        //- Calculate face flip map
+        void calcFaceFlipMap() const;
 
         //- Adapt nCellsUsingFace for coupled faces becoming 'uncoupled'.
         void doCoupledPatches
@@ -115,24 +109,41 @@ class fvMeshSubset
             labelList& nCellsUsingFace
         ) const;
 
+        //- Forwarding to Foam::removeCells
+        void removeCellsImpl
+        (
+            const bitSet& cellsToRemove,
+            const labelList& exposedFaces,
+            const labelList& patchIDs,
+            const bool syncCouples
+        );
+
+
         //- Subset of subset
-        static labelList subset
+        static labelList subsetSubset
         (
             const label nElems,
-            const labelList& selectedElements,
-            const labelList& subsetMap
+            const labelUList& selectedElements,
+            const labelUList& subsetMap
         );
 
         //- Create zones for submesh
         void subsetZones();
 
         //- Helper: extract cells-to-remove from cells-to-keep
-        labelList getCellsToRemove
+        bitSet getCellsToRemove
         (
-            const labelList& region,
-            const label currentRegion
+            const bitSet& selectedCells
         ) const;
 
+        //- Helper: extract cells-to-remove from cells-to-keep
+        bitSet getCellsToRemove
+        (
+            const label regioni,
+            const labelUList& regions
+        ) const;
+
+
         //- No copy construct
         fvMeshSubset(const fvMeshSubset&) = delete;
 
@@ -145,173 +156,316 @@ public:
     // Constructors
 
         //- Construct given a mesh to subset
-        explicit fvMeshSubset(const fvMesh&);
+        explicit fvMeshSubset(const fvMesh& baseMesh);
+
+        //- Construct for a cell-subset of the given mesh
+        //  See setCellSubset() for more details.
+        fvMeshSubset
+        (
+            const fvMesh& baseMesh,
+            const bitSet& selectedCells,
+            const label patchID = -1,
+            const bool syncPar = true
+        );
+
+        //- Construct for a cell-subset of the given mesh
+        //  See setCellSubset() for more details.
+        fvMeshSubset
+        (
+            const fvMesh& baseMesh,
+            const labelUList& selectedCells,
+            const label patchID = -1,
+            const bool syncPar = true
+        );
+
+        //- Construct for a cell-subset of the given mesh
+        //  See setCellSubset() for more details.
+        fvMeshSubset
+        (
+            const fvMesh& baseMesh,
+            const labelHashSet& selectedCells,
+            const label patchID = -1,
+            const bool syncPar = true
+        );
+
+        //- Construct for a cell-subset of the given mesh
+        //  See setCellSubset() for more details.
+        fvMeshSubset
+        (
+            const fvMesh& baseMesh,
+            const label regioni,
+            const labelUList& regions,
+            const label patchID = -1,
+            const bool syncPar = true
+        );
 
 
     // Member Functions
 
-        // Edit
+    // Access
 
-            //- Set the subset from all cells with region == currentRegion.
-            //  Create "oldInternalFaces" patch for exposed
-            //  internal faces (patchID==-1) or use supplied patch.
-            //  Handles coupled patches by if necessary making coupled patch
-            //  face part of patchID (so uncoupled)
-            void setLargeCellSubset
-            (
-                const labelList& region,
-                const label currentRegion,
-                const label patchID = -1,
-                const bool syncCouples = true
-            );
+        //- Original mesh
+        inline const fvMesh& baseMesh() const;
 
-            //- setLargeCellSubset but only marking certain cells
-            void setLargeCellSubset
-            (
-                const labelUList& globalCellMap,
-                const label patchID = -1,
-                const bool syncPar = true
-            );
+        //- Have subMesh?
+        inline bool hasSubMesh() const;
 
-            //- setLargeCellSubset but with labelHashSet.
-            void setLargeCellSubset
-            (
-                const labelHashSet& globalCellMap,
-                const label patchID = -1,
-                const bool syncPar = true
-            );
+        //- Return reference to subset mesh
+        inline const fvMesh& subMesh() const;
 
+        //- Return reference to subset mesh
+        inline fvMesh& subMesh();
 
-            //- Two step subsetting
+        //- Return point map
+        inline const labelList& pointMap() const;
 
-                //- Get labels of exposed faces.
-                //  These are
-                //  - internal faces that become boundary faces
-                //  - coupled faces that become uncoupled (since one of the
-                //    sides gets deleted)
-                labelList getExposedFaces
-                (
-                    const labelList& region,
-                    const label currentRegion,
-                    const bool syncCouples = true
-                ) const;
+        //- Return face map
+        inline const labelList& faceMap() const;
 
-                //- For every exposed face (from above getExposedFaces)
-                //  used supplied (existing!) patch
-                void setLargeCellSubset
-                (
-                    const labelList& region,
-                    const label currentRegion,
-                    const labelList& exposedFaces,
-                    const labelList& patchIDs,
-                    const bool syncCouples = true
-                );
+        //- Return face map with sign to encode flipped faces
+        inline const labelList& faceFlipMap() const;
 
+        //- Return cell map
+        inline const labelList& cellMap() const;
 
-        // Access
+        //- Return patch map
+        inline const labelList& patchMap() const;
 
-            //- Original mesh
-            const fvMesh& baseMesh() const
-            {
-                return baseMesh_;
-            }
 
-            //- Have subMesh?
-            bool hasSubMesh() const;
+    // Edit
 
-            //- Return reference to subset mesh
-            const fvMesh& subMesh() const;
+        //- Reset maps and subsetting
+        void clear();
 
-            fvMesh& subMesh();
 
-            //- Return point map
-            const labelList& pointMap() const;
+        //- Define the cell sub