ENH: multiple zone selection for fvMeshSubsetProxy (#973)

- handle tmp fields in interpolate methods

- special method interpolateInternal() for creating a volume field
  with zero-gradient treatment for patches from an internal field.

  This method was previously also called interpolate(), but that
  masked the ability to subset the internal field only.

  Ensight output needs the volume field:
      uses interpolateInternal().

  VTK output has separate handling of internal and patch fields:
      uses interpolate().

ENH: added fvMeshSubset mesh() method for baseMesh or subMesh.

- simplies coding when the fvMeshSubset may or may not be in active use.

ENH: update foamToEnsight to use newer methods in wrapped form

- static interpolate functions with renaming for manual use with
  fvMeshSubset (when fvMeshSubsetProxy may be too limiting in
  functionality)

applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C

@@ -88,15 +88,49 @@ Note
 #include "ensightGeoFile.H"
 #include "ensightMesh.H"
 #include "ensightOutput.H"
+#include "fvMeshSubsetProxy.H"
 
 // local files
-#include "fvMeshSubsetProxy.H"
 #include "ensightOutputCloud.H"
 
 #include "memInfo.H"
 
 using namespace Foam;
 
+
+//- Get the field and subset it
+template<class GeoField>
+tmp<GeoField> getField(IOobject& io, const fvMeshSubsetProxy& proxy)
+{
+    auto tfield = tmp<GeoField>::New(io, proxy.baseMesh());
+    return proxy.interpolate(tfield);
+}
+
+
+//- Get the field and subset it, or return nullptr
+template<class GeoField>
+tmp<GeoField> getField(const IOobject* io, const fvMeshSubsetProxy& proxy)
+{
+    if (io)
+    {
+        auto tfield = tmp<GeoField>::New(*io, proxy.baseMesh());
+        return proxy.interpolate(tfield);
+    }
+
+    return nullptr;
+}
+
+
+//- Get internal field and make it a zero-gradient volume field with subsetting
+template<class GeoField>
+tmp<GeoField>
+getZeroGradInternalField(IOobject& io, const fvMeshSubsetProxy& proxy)
+{
+    auto tfield = tmp<typename GeoField::Internal>::New(io, proxy.baseMesh());
+    return proxy.interpolateInternal(tfield);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -256,12 +290,12 @@ int main(int argc, char *argv[])
     if (args.readIfPresent("cellZone", cellZoneName))
     {
         Info<< "Converting cellZone " << cellZoneName
-            << " only. Places new outside faces into \"oldInternalFaces\"."
+            << " only, with new outside faces as \"oldInternalFaces\"."
             << nl;
     }
 
-    // fvMeshSubset is ignored if cellZoneName is empty
+    // Ignored (unproxied) if cellZoneName is empty
-    fvMeshSubsetProxy meshRef(mesh, fvMeshSubsetProxy::ZONE, cellZoneName);
+    fvMeshSubsetProxy meshProxy(mesh, fvMeshSubsetProxy::ZONE, cellZoneName);
 
     //
     // Open new ensight case file, initialize header etc.
@@ -275,7 +309,7 @@ int main(int argc, char *argv[])
 
 
     // Construct the Ensight mesh
-    ensightMesh ensMesh(meshRef.mesh(), writeOpts);
+    ensightMesh ensMesh(meshProxy.mesh(), writeOpts);
 
     if (Pstream::master())
     {
@@ -344,7 +378,7 @@ int main(int argc, char *argv[])
         polyMesh::readUpdateState meshState = mesh.readUpdate();
         if (meshState != polyMesh::UNCHANGED)
         {
-            meshRef.correct();
+            meshProxy.correct();
             ensMesh.expire();
             ensMesh.correct();
         }
@@ -396,10 +430,9 @@ int main(int argc, char *argv[])
                         fieldName
                     );
 
-                    volScalarField vf(fieldObject, mesh);
                     wrote = ensightOutput::writeField<scalar>
                     (
-                        meshRef.interpolate(vf),
+                        getField<volScalarField>(fieldObject, meshProxy),
                         ensMesh,
                         os,
                         nodeValues
@@ -412,10 +445,9 @@ int main(int argc, char *argv[])
                         fieldName
                     );
 
-                    volVectorField vf(fieldObject, mesh);
                     wrote = ensightOutput::writeField<vector>
                     (
-                        meshRef.interpolate(vf),
+                        getField<volVectorField>(fieldObject, meshProxy),
                         ensMesh,
                         os,
                         nodeValues
@@ -428,10 +460,13 @@ int main(int argc, char *argv[])
                         fieldObject.name()
                     );
 
-                    volSphericalTensorField vf(fieldObject, mesh);
                     wrote = ensightOutput::writeField<sphericalTensor>
                     (
-                        meshRef.interpolate(vf),
+                        getField<volSphericalTensorField>
+                        (
+                            fieldObject,
+                            meshProxy
+                        ),
                         ensMesh,
                         os,
                         nodeValues
@@ -444,10 +479,13 @@ int main(int argc, char *argv[])
                         fieldName
                     );
 
-                    volSymmTensorField vf(fieldObject, mesh);
                     wrote = ensightOutput::writeField<symmTensor>
                     (
-                        meshRef.interpolate(vf),
+                        getField<volSymmTensorField>
+                        (
+                            fieldObject,
+                            meshProxy
+                        ),
                         ensMesh,
                         os,
                         nodeValues
@@ -460,10 +498,9 @@ int main(int argc, char *argv[])
                         fieldName
                     );
 
-                    volTensorField vf(fieldObject, mesh);
                     wrote = ensightOutput::writeField<tensor>
                     (
-                        meshRef.interpolate(vf),
+                        getField<volTensorField>(fieldObject, meshProxy),
                         ensMesh,
                         os,
                         nodeValues
@@ -481,14 +518,13 @@ int main(int argc, char *argv[])
                         fieldName
                     );
 
-                    volScalarField::Internal df
-                    (
-                        fieldObject,
-                        mesh
-                    );
                     wrote = ensightOutput::writeField<scalar>
                     (
-                        meshRef.interpolate<scalar>(df),
+                        getZeroGradInternalField<volScalarField>
+                        (
+                            fieldObject,
+                            meshProxy
+                        ),
                         ensMesh,
                         os,
                         nodeValues
@@ -505,14 +541,13 @@ int main(int argc, char *argv[])
                         fieldName
                     );
 
-                    volVectorField::Internal df
-                    (
-                        fieldObject,
-                        mesh
-                    );
                     wrote = ensightOutput::writeField<vector>
                     (
-                        meshRef.interpolate<vector>(df),
+                        getZeroGradInternalField<volVectorField>
+                        (
+                            fieldObject,
+                            meshProxy
+                        ),
                         ensMesh,
                         os,
                         nodeValues
@@ -529,14 +564,13 @@ int main(int argc, char *argv[])
                         fieldName
                     );
 
-                    volSphericalTensorField::Internal df
-                    (
-                        fieldObject,
-                        mesh
-                    );
                     wrote = ensightOutput::writeField<sphericalTensor>
                     (
-                        meshRef.interpolate<sphericalTensor>(df),
+                        getZeroGradInternalField<volSphericalTensorField>
+                        (
+                            fieldObject,
+                            meshProxy
+                        ),
                         ensMesh,
                         os,
                         nodeValues
@@ -553,14 +587,13 @@ int main(int argc, char *argv[])
                         fieldName
                     );
 
-                    volSymmTensorField::Internal df
-                    (
-                        fieldObject,
-                        mesh
-                    );
                     wrote = ensightOutput::writeField<symmTensor>
                     (
-                        meshRef.interpolate<symmTensor>(df),
+                        getZeroGradInternalField<volSymmTensorField>
+                        (
+                            fieldObject,
+                            meshProxy
+                        ),
                         ensMesh,
                         os,
                         nodeValues
@@ -577,14 +610,13 @@ int main(int argc, char *argv[])
                         fieldName
                     );
 
-                    volTensorField::Internal df
-                    (
-                        fieldObject,
-                        mesh
-                    );
                     wrote = ensightOutput::writeField<tensor>
                     (
-                        meshRef.interpolate<tensor>(df),
+                        getZeroGradInternalField<volTensorField>
+                        (
+                            fieldObject,
+                            meshProxy
+                        ),
                         ensMesh,
                         os,
                         nodeValues

applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C

@@ -554,7 +554,7 @@ int main(int argc, char *argv[])
     instantList timeDirs = timeSelector::select0(runTime, args);
 
     // Mesh wrapper: does subsetting and decomposition
-    fvMeshSubsetProxy meshRef(mesh, cellSubsetType, cellSubsetName);
+    fvMeshSubsetProxy meshProxy(mesh, cellSubsetType, cellSubsetName);
 
     // Collect decomposition information etc.
     vtk::vtuCells vtuMeshCells(fmtType, decomposePoly);
@@ -616,9 +616,9 @@ int main(int argc, char *argv[])
 
         // Check for new polyMesh/ and update mesh, fvMeshSubset and cell
         // decomposition.
-        polyMesh::readUpdateState meshState = meshRef.readUpdate();
+        polyMesh::readUpdateState meshState = meshProxy.readUpdate();
 
-        const fvMesh& mesh = meshRef.mesh();
+        const fvMesh& mesh = meshProxy.mesh();
         if
         (
             meshState == polyMesh::TOPO_CHANGE
@@ -649,7 +649,7 @@ int main(int argc, char *argv[])
 
             vtk::writeFaceSet
             (
-                meshRef.mesh(),
+                meshProxy.mesh(),
                 set,
                 outputName,
                 fmtType
@@ -677,7 +677,7 @@ int main(int argc, char *argv[])
 
             vtk::writePointSet
             (
-                meshRef.mesh(),
+                meshProxy.mesh(),
                 set,
                 outputName,
                 fmtType
@@ -718,8 +718,8 @@ int main(int argc, char *argv[])
         {
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 vScalarFld
@@ -728,8 +728,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 vVectorFld
@@ -738,8 +738,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 vSphTensorf
@@ -748,8 +748,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 vSymTensorFld
@@ -758,8 +758,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 vTensorFld
@@ -797,8 +797,8 @@ int main(int argc, char *argv[])
         {
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 dScalarFld
@@ -807,8 +807,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 dVectorFld
@@ -817,8 +817,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 dSphTensorFld
@@ -827,8 +827,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 dSymTensorFld
@@ -837,8 +837,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 dTensorFld
@@ -865,7 +865,7 @@ int main(int argc, char *argv[])
                 const bool throwing = FatalError.throwExceptions();
                 try
                 {
-                    aMeshPtr.reset(new faMesh(meshRef.baseMesh()));
+                    aMeshPtr.reset(new faMesh(meshProxy.baseMesh()));
                 }
                 catch (Foam::error& err)
                 {
@@ -992,11 +992,11 @@ int main(int argc, char *argv[])
             ).size()
         )
         {
-            const pointMesh& ptMesh = pointMesh::New(meshRef.baseMesh());
+            const pointMesh& ptMesh = pointMesh::New(meshProxy.baseMesh());
 
             readFields
             (
-                meshRef,
+                meshProxy,
                 ptMesh,
                 objects,
                 selectedFields,
@@ -1006,7 +1006,7 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
+                meshProxy,
                 ptMesh,
                 objects,
                 selectedFields,
@@ -1016,7 +1016,7 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
+                meshProxy,
                 ptMesh,
                 objects,
                 selectedFields,
@@ -1026,7 +1026,7 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
+                meshProxy,
                 ptMesh,
                 objects,
                 selectedFields,
@@ -1036,7 +1036,7 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
+                meshProxy,
                 ptMesh,
                 objects,
                 selectedFields,
@@ -1057,7 +1057,7 @@ int main(int argc, char *argv[])
             if (vtuMeshCells.empty())
             {
                 // subMesh or baseMesh
-                vtuMeshCells.reset(meshRef.mesh());
+                vtuMeshCells.reset(meshProxy.mesh());
             }
 
             // Create file and write header
@@ -1073,7 +1073,7 @@ int main(int argc, char *argv[])
             // Write mesh
             vtk::internalWriter writer
             (
-                meshRef.mesh(),
+                meshProxy.mesh(),
                 vtuMeshCells,
                 outputName,
                 fmtType
@@ -1147,8 +1147,8 @@ int main(int argc, char *argv[])
             PtrList<const surfaceScalarField> sScalarFld;
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 sScalarFld
@@ -1158,8 +1158,8 @@ int main(int argc, char *argv[])
             PtrList<const surfaceVectorField> sVectorFld;
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 sVectorFld
@@ -1199,7 +1199,7 @@ int main(int argc, char *argv[])
 
                 vtk::writeSurfFields
                 (
-                    meshRef.mesh(),
+                    meshProxy.mesh(),
                     outputName,
                     fmtType,
                     sVectorFld
@@ -1223,7 +1223,7 @@ int main(int argc, char *argv[])
             fileName outputName
             (
                 fvPath/"allPatches"
-              / (meshRef.useSubMesh() ? cellSubsetName : "allPatches")
+              / (meshProxy.useSubMesh() ? cellSubsetName : "allPatches")
               + "_"
               + timeDesc
             );
@@ -1232,7 +1232,7 @@ int main(int argc, char *argv[])
 
             vtk::patchWriter writer
             (
-                meshRef.mesh(),
+                meshProxy.mesh(),
                 outputName,
                 fmtType,
                 nearCellValue,
@@ -1293,7 +1293,7 @@ int main(int argc, char *argv[])
                 fileName outputName
                 (
                     fvPath/pp.name()
-                  / (meshRef.useSubMesh() ? cellSubsetName : pp.name())
+                  / (meshProxy.useSubMesh() ? cellSubsetName : pp.name())
                   + "_"
                   + timeDesc
                 );
@@ -1302,7 +1302,7 @@ int main(int argc, char *argv[])
 
                 vtk::patchWriter writer
                 (
-                    meshRef.mesh(),
+                    meshProxy.mesh(),
                     outputName,
                     fmtType,
                     nearCellValue,
@@ -1368,8 +1368,8 @@ int main(int argc, char *argv[])
             PtrList<const surfaceScalarField> sScalarFld;
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 sScalarFld
@@ -1379,8 +1379,8 @@ int main(int argc, char *argv[])
             PtrList<const surfaceVectorField> sVectorFld;
             readFields
             (
-                meshRef,
+                meshProxy,
-                meshRef.baseMesh(),
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 sVectorFld
@@ -1394,7 +1394,7 @@ int main(int argc, char *argv[])
                 fileName outputName =
                 (
                     fvPath/fz.name()
-                  / (meshRef.useSubMesh() ? cellSubsetName : fz.name())
+                  / (meshProxy.useSubMesh() ? cellSubsetName : fz.name())
                   + "_"
                   + timeDesc
                 );
@@ -1528,7 +1528,7 @@ int main(int argc, char *argv[])
 
                 vtk::lagrangianWriter writer
                 (
-                    meshRef.mesh(),
+                    meshProxy.mesh(),
                     cloudName,
                     outputName,
                     fmtType
@@ -1553,7 +1553,7 @@ int main(int argc, char *argv[])
             {
                 vtk::lagrangianWriter writer
                 (
-                    meshRef.mesh(),
+                    meshProxy.mesh(),
                     cloudName,
                     outputName,
                     fmtType,

src/dynamicMesh/fvMeshSubset/fvMeshSubset.H

@@ -213,6 +213,9 @@ public:
         //- Original mesh
         inline const fvMesh& baseMesh() const;
 
+        //- Return baseMesh or subMesh, depending on the current state.
+        inline const fvMesh& mesh() const;
+
         //- Have subMesh?
         inline bool hasSubMesh() const;
 
@@ -244,10 +247,10 @@ public:
         void clear();
 
 
-        //- Define the cell subset based on the selectedCells.
+        //- Define cell subset based on the selectedCells.
         //  Create "oldInternalFaces" patch for exposed
         //  internal faces (patchID==-1) or use supplied patch.
-        //  Handles coupled patches by if necessary making coupled patch
+        //  Handles coupled patches if necessary by making coupled patch
         //  face part of patchID (so uncoupled)
         void setCellSubset
         (
@@ -256,7 +259,7 @@ public:
             const bool syncPar = true
         );
 
-        //- Define the cell subset, using the specified cells
+        //- Define cell subset, using the specified cells
         //- to define the selection
         void setCellSubset
         (
@@ -265,7 +268,7 @@ public:
             const bool syncPar = true
         );
 
-        //- Define the cell subset, using the specified cells
+        //- Define cell subset, using the specified cells
         //- labelHashSet to define the selection
         void setCellSubset
         (
@@ -274,8 +277,7 @@ public:
             const bool syncPar = true
         );
 
-        //- Define the cell subset, using the cells for which
+        //- Define cell subset, using the cells for which region == regioni.
-        //  region == regioni.
         void setCellSubset
         (
             const label regioni,
@@ -285,7 +287,7 @@ public:
         );
 
 
-    // Two step subsetting
+    // Two-step subsetting
 
         //- Get labels of exposed faces.
         //  These are
@@ -332,7 +334,7 @@ public:
         );
 
 
-    // Field mapping
+    // Field Mapping
 
         //- Map volume field
         template<class Type>
@@ -341,9 +343,9 @@ public:
         (
             const GeometricField<Type, fvPatchField, volMesh>&,
             const fvMesh& sMesh,
-            const labelList& patchMap,
+            const labelUList& patchMap,
-            const labelList& cellMap,
+            const labelUList& cellMap,
-            const labelList& faceMap
+            const labelUList& faceMap
         );
 
         template<class Type>
@@ -361,9 +363,9 @@ public:
         (
             const GeometricField<Type, fvsPatchField, surfaceMesh>&,
             const fvMesh& sMesh,
-            const labelList& patchMap,
+            const labelUList& patchMap,
-            const labelList& cellMap,
+            const labelUList& cellMap,
-            const labelList& faceMap
+            const labelUList& faceMap
         );
 
         template<class Type>
@@ -380,8 +382,8 @@ public:
         (
             const GeometricField<Type, pointPatchField, pointMesh>&,
             const pointMesh& sMesh,
-            const labelList& patchMap,
+            const labelUList& patchMap,
-            const labelList& pointMap
+            const labelUList& pointMap
         );
 
         template<class Type>
@@ -398,7 +400,7 @@ public:
         (
             const DimensionedField<Type, volMesh>&,
             const fvMesh& sMesh,
-            const labelList& cellMap
+            const labelUList& cellMap
         );
 
         template<class Type>

src/dynamicMesh/fvMeshSubset/fvMeshSubsetI.H

@@ -31,6 +31,12 @@ inline const Foam::fvMesh& Foam::fvMeshSubset::baseMesh() const
 }
 
 
+inline const Foam::fvMesh& Foam::fvMeshSubset::mesh() const
+{
+    return fvMeshSubsetPtr_.valid() ? *fvMeshSubsetPtr_ : baseMesh_;
+}
+
+
 inline bool Foam::fvMeshSubset::hasSubMesh() const
 {
     return fvMeshSubsetPtr_.valid();

src/dynamicMesh/fvMeshSubset/fvMeshSubsetInterpolate.C

@@ -33,20 +33,18 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
 template<class Type>
-tmp<GeometricField<Type, fvPatchField, volMesh>>
+Foam::tmp
-fvMeshSubset::interpolate
+<
+    Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>
+>
+Foam::fvMeshSubset::interpolate
 (
     const GeometricField<Type, fvPatchField, volMesh>& vf,
     const fvMesh& sMesh,
-    const labelList& patchMap,
+    const labelUList& patchMap,
-    const labelList& cellMap,
+    const labelUList& cellMap,
-    const labelList& faceMap
+    const labelUList& faceMap
 )
 {
     // 1. Create the complete field with dummy patch fields
@@ -84,7 +82,7 @@ fvMeshSubset::interpolate
         }
     }
 
-    auto tresF = tmp<GeometricField<Type, fvPatchField, volMesh>>::New
+    auto tresult = tmp<GeometricField<Type, fvPatchField, volMesh>>::New
     (
         IOobject
         (
@@ -99,22 +97,24 @@ fvMeshSubset::interpolate
         Field<Type>(vf.primitiveField(), cellMap),
         patchFields
     );
-    auto& resF = tresF.ref();
+    auto& result = tresult.ref();
-    resF.oriented() = vf.oriented();
+    result.oriented() = vf.oriented();
 
 
     // 2. Change the fvPatchFields to the correct type using a mapper
     //  constructor (with reference to the now correct internal field)
 
-    auto& bf = resF.boundaryFieldRef();
+    auto& bf = result.boundaryFieldRef();
 
     forAll(bf, patchi)
     {
-        if (patchMap[patchi] != -1)
+        const label basePatchId = patchMap[patchi];
+
+        if (basePatchId != -1)
         {
             // Construct addressing
             const fvPatch& subPatch = sMesh.boundary()[patchi];
-            const fvPatch& basePatch = vf.mesh().boundary()[patchMap[patchi]];
+            const fvPatch& basePatch = vf.mesh().boundary()[basePatchId];
             const label baseStart = basePatch.start();
             const label baseSize = basePatch.size();
 
@@ -122,7 +122,7 @@ fvMeshSubset::interpolate
 
             forAll(directAddressing, i)
             {
-                label baseFacei = faceMap[subPatch.start()+i];
+                const label baseFacei = faceMap[subPatch.start()+i];
 
                 if (baseFacei >= baseStart && baseFacei < baseStart+baseSize)
                 {
@@ -141,22 +141,25 @@ fvMeshSubset::interpolate
                 patchi,
                 fvPatchField<Type>::New
                 (
-                    vf.boundaryField()[patchMap[patchi]],
+                    vf.boundaryField()[basePatchId],
                     subPatch,
-                    resF(),
+                    result(),
                     directFvPatchFieldMapper(directAddressing)
                 )
             );
         }
     }
 
-    return tresF;
+    return tresult;
 }
 
 
 template<class Type>
-tmp<GeometricField<Type, fvPatchField, volMesh>>
+Foam::tmp
-fvMeshSubset::interpolate
+<
+    Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>
+>
+Foam::fvMeshSubset::interpolate
 (
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
@@ -173,14 +176,17 @@ fvMeshSubset::interpolate
 
 
 template<class Type>
-tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+Foam::tmp
-fvMeshSubset::interpolate
+<
+    Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>
+>
+Foam::fvMeshSubset::interpolate
 (
     const GeometricField<Type, fvsPatchField, surfaceMesh>& vf,
     const fvMesh& sMesh,
-    const labelList& patchMap,
+    const labelUList& patchMap,
-    const labelList& cellMap,
+    const labelUList& cellMap,
-    const labelList& faceMap
+    const labelUList& faceMap
 )
 {
     // 1. Create the complete field with dummy patch fields
@@ -219,7 +225,7 @@ fvMeshSubset::interpolate
     }
 
     // Create the complete field from the pieces
-    auto tresF = tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>::New
+    auto tresult = tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>::New
     (
         IOobject
         (
@@ -238,14 +244,14 @@ fvMeshSubset::interpolate
         ),
         patchFields
     );
-    auto& resF = tresF.ref();
+    auto& result = tresult.ref();
-    resF.oriented() = vf.oriented();
+    result.oriented() = vf.oriented();
 
 
     // 2. Change the fvsPatchFields to the correct type using a mapper
     //  constructor (with reference to the now correct internal field)
 
-    auto& bf = resF.boundaryFieldRef();
+    auto& bf = result.boundaryFieldRef();
 
     forAll(bf, patchi)
     {
@@ -283,7 +289,7 @@ fvMeshSubset::interpolate
                 (
                     vf.boundaryField()[patchMap[patchi]],
                     subPatch,
-                    resF(),
+                    result(),
                     directFvPatchFieldMapper(directAddressing)
                 )
             );
@@ -327,13 +333,16 @@ fvMeshSubset::interpolate
         }
     }
 
-    return tresF;
+    return tresult;
 }
 
 
 template<class Type>
-tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+Foam::tmp
-fvMeshSubset::interpolate
+<
+    Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>
+>
+Foam::fvMeshSubset::interpolate
 (
     const GeometricField<Type, fvsPatchField, surfaceMesh>& sf
 ) const
@@ -350,13 +359,16 @@ fvMeshSubset::interpolate
 
 
 template<class Type>
-tmp<GeometricField<Type, pointPatchField, pointMesh>>
+Foam::tmp
-fvMeshSubset::interpolate
+<
+    Foam::GeometricField<Type, Foam::pointPatchField, Foam::pointMesh>
+>
+Foam::fvMeshSubset::interpolate
 (
     const GeometricField<Type, pointPatchField, pointMesh>& vf,
     const pointMesh& sMesh,
-    const labelList& patchMap,
+    const labelUList& patchMap,
-    const labelList& pointMap
+    const labelUList& pointMap
 )
 {
     // 1. Create the complete field with dummy patch fields
@@ -395,7 +407,7 @@ fvMeshSubset::interpolate
     }
 
     // Create the complete field from the pieces
-    auto tresF = tmp<GeometricField<Type, pointPatchField, pointMesh>>::New
+    auto tresult = tmp<GeometricField<Type, pointPatchField, pointMesh>>::New
     (
         IOobject
         (
@@ -410,14 +422,14 @@ fvMeshSubset::interpolate
         Field<Type>(vf.primitiveField(), pointMap),
         patchFields
     );
-    auto& resF = tresF.ref();
+    auto& result = tresult.ref();
-    resF.oriented() = vf.oriented();
+    result.oriented() = vf.oriented();
 
 
     // 2. Change the pointPatchFields to the correct type using a mapper
     //  constructor (with reference to the now correct internal field)
 
-    auto& bf = resF.boundaryFieldRef();
+    auto& bf = result.boundaryFieldRef();
 
     forAll(bf, patchi)
     {
@@ -466,20 +478,23 @@ fvMeshSubset::interpolate
                 (
                     vf.boundaryField()[patchMap[patchi]],
                     subPatch,
-                    resF(),
+                    result(),
                     directPointPatchFieldMapper(directAddressing)
                 )
             );
         }
     }
 
-    return tresF;
+    return tresult;
 }
 
 
 template<class Type>
-tmp<GeometricField<Type, pointPatchField, pointMesh>>
+Foam::tmp
-fvMeshSubset::interpolate
+<
+    Foam::GeometricField<Type, Foam::pointPatchField, Foam::pointMesh>
+>
+Foam::fvMeshSubset::interpolate
 (
     const GeometricField<Type, pointPatchField, pointMesh>& sf
 ) const
@@ -495,16 +510,18 @@ fvMeshSubset::interpolate
 
 
 template<class Type>
-tmp<DimensionedField<Type, volMesh>>
+Foam::tmp
-fvMeshSubset::interpolate
+<
+    Foam::DimensionedField<Type, Foam::volMesh>
+>
+Foam::fvMeshSubset::interpolate
 (
     const DimensionedField<Type, volMesh>& df,
     const fvMesh& sMesh,
-    const labelList& cellMap
+    const labelUList& cellMap
 )
 {
-    // Create the complete field from the pieces
+    auto tresult = tmp<DimensionedField<Type, volMesh>>::New
-    auto tresF = tmp<DimensionedField<Type, volMesh>>::New
     (
         IOobject
         (
@@ -518,15 +535,18 @@ fvMeshSubset::interpolate
         df.dimensions(),
         Field<Type>(df, cellMap)
     );
-    tresF.ref().oriented() = df.oriented();
+    tresult.ref().oriented() = df.oriented();
 
-    return tresF;
+    return tresult;
 }
 
 
 template<class Type>
-tmp<DimensionedField<Type, volMesh>>
+Foam::tmp
-fvMeshSubset::interpolate
+<
+    Foam::DimensionedField<Type, Foam::volMesh>
+>
+Foam::fvMeshSubset::interpolate
 (
     const DimensionedField<Type, volMesh>& df
 ) const
@@ -535,8 +555,4 @@ fvMeshSubset::interpolate
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/dynamicMesh/fvMeshSubsetProxy/fvMeshSubsetProxy.C

@@ -30,18 +30,19 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::fvMeshSubsetProxy::fvMeshSubsetProxy
+Foam::fvMeshSubsetProxy::fvMeshSubsetProxy(fvMesh& baseMesh)
-(
-    fvMesh& baseMesh
-)
 :
     baseMesh_(baseMesh),
     subsetter_(baseMesh),
+    exposedPatchId_(-1),
     type_(NONE),
     name_(),
-    exposedPatchId_(-1)
+    names_()
 {
+    if (useSubMesh())
+    {
         correct();
+    }
 }
 
 
@@ -49,17 +50,74 @@ Foam::fvMeshSubsetProxy::fvMeshSubsetProxy
 (
     fvMesh& baseMesh,
     const subsetType type,
-    const word& name,
+    const word& selectionName,
-    const label exposedPatchId
+    label exposedPatchId
 )
 :
     baseMesh_(baseMesh),
     subsetter_(baseMesh),
-    type_(name.empty() ? NONE : type),
+    exposedPatchId_(exposedPatchId),
-    name_(name),
+    type_(selectionName.empty() ? NONE : type),
-    exposedPatchId_(exposedPatchId)
+    name_(),
+    names_()
 {
+    if (type_ == ZONES)
+    {
+        // Populate wordRes for ZONES
+        names_.resize(1);
+        names_.first() = selectionName;
+    }
+    else if (type_ != NONE)
+    {
+        name_ = selectionName;
+    }
+
+    if (useSubMesh())
+    {
         correct();
+    }
+}
+
+
+Foam::fvMeshSubsetProxy::fvMeshSubsetProxy
+(
+    fvMesh& baseMesh,
+    const wordRes& zoneNames,
+    label exposedPatchId
+)
+:
+    baseMesh_(baseMesh),
+    subsetter_(baseMesh),
+    exposedPatchId_(exposedPatchId),
+    type_(ZONES),
+    name_(),
+    names_(zoneNames)
+{
+    if (useSubMesh())
+    {
+        correct();
+    }
+}
+
+
+Foam::fvMeshSubsetProxy::fvMeshSubsetProxy
+(
+    fvMesh& baseMesh,
+    wordRes&& zoneNames,
+    label exposedPatchId
+)
+:
+    baseMesh_(baseMesh),
+    subsetter_(baseMesh),
+    exposedPatchId_(exposedPatchId),
+    type_(ZONES),
+    name_(),
+    names_(std::move(zoneNames))
+{
+    if (useSubMesh())
+    {
+        correct();
+    }
 }
 
 
@@ -67,6 +125,16 @@ Foam::fvMeshSubsetProxy::fvMeshSubsetProxy
 
 void Foam::fvMeshSubsetProxy::correct(bool verbose)
 {
+    if (type_ == NONE)
+    {
+        subsetter_.clear();
+        return;
+    }
+
+    const label nCells = baseMesh_.nCells();
+
+    bitSet selectedCells;
+
     if (type_ == SET)
     {
         if (verbose)
@@ -74,11 +142,13 @@ void Foam::fvMeshSubsetProxy::correct(bool verbose)
             Info<< "Subsetting mesh based on cellSet " << name_ << endl;
         }
 
-        subsetter_.setCellSubset
+        cellSet cset(baseMesh_, name_);
-        (
+
-            cellSet(baseMesh_, name_),
+        selectedCells.resize(nCells);
-            exposedPatchId_
+        for (const label idx : cset)
-        );
+        {
+            selectedCells.set(idx);
+        }
     }
     else if (type_ == ZONE)
     {
@@ -87,12 +157,21 @@ void Foam::fvMeshSubsetProxy::correct(bool verbose)
             Info<< "Subsetting mesh based on cellZone " << name_ << endl;
         }
 
-        subsetter_.setCellSubset
+        selectedCells.resize(nCells);
-        (
+        selectedCells.set(baseMesh_.cellZones()[name_]);
-            baseMesh_.cellZones()[name_],
-            exposedPatchId_
-        );
     }
+    else if (type_ == ZONES)
+    {
+        if (verbose)
+        {
+            Info<< "Subsetting mesh based on cellZones "
+                << flatOutput(names_) << endl;
+        }
+
+        selectedCells = baseMesh_.cellZones().selection(names_);
+    }
+
+    subsetter_.setCellSubset(selectedCells, exposedPatchId_);
 }
 
 

src/dynamicMesh/fvMeshSubsetProxy/fvMeshSubsetProxy.H

@@ -37,6 +37,7 @@ SourceFiles
 #ifndef fvMeshSubsetProxy_H
 #define fvMeshSubsetProxy_H
 
+#include "wordRes.H"
 #include "fvMeshSubset.H"
 #include "zeroGradientFvPatchField.H"
 
@@ -56,33 +57,37 @@ class fvMeshSubsetProxy
 {
 public:
 
-        //- Internal book-keeping for subset type
+        //- Internal bookkeeping for subset type
         enum subsetType
         {
-            NONE = 0,   //<! Not a subset
+            NONE,       //<! No subset
-            SET,        //<! Using a cellSet for the subset
+            SET,        //<! Subset with a cellSet
-            ZONE        //<! Using a cellZone for the subset
+            ZONE,       //<! Subset with a cellZone
+            ZONES       //<! Subset with multiple cellZones
         };
 
 
 private:
 
-    // Private data
+    // Private Data
 
         //- Reference to mesh
         fvMesh& baseMesh_;
 
-        //- Subsetting engine + sub-fvMesh
+        //- Subsetting engine
         fvMeshSubset subsetter_;
 
-        //- Subset type if any.
-        const subsetType type_;
-
-        //- Name of current cellSet/cellZone (or empty)
-        const word name_;
-
         //- Patch ID for exposed internal faces
-        const label exposedPatchId_;
+        label exposedPatchId_;
+
+        //- The subsetting type
+        subsetType type_;
+
+        //- Name of the cellSet/cellZone (or empty)
+        word name_;
+
+        //- Selection for multiple cell zones
+        wordRes names_;
 
 
     // Private Member Functions
@@ -98,16 +103,32 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct a pass-through proxy. No correct() invoked or required.
-        fvMeshSubsetProxy(fvMesh& baseMesh);
+        explicit fvMeshSubsetProxy(fvMesh& baseMesh);
 
         //- Construct from components
         fvMeshSubsetProxy
         (
             fvMesh& baseMesh,
-            const subsetType,
+            const subsetType type,
-            const word& name,
+            const word& selectionName,
-            const label exposedPatchId = -1
+            label exposedPatchId = -1
+        );
+
+        //- Construct from components. The subsetType is ZONES.
+        fvMeshSubsetProxy
+        (
+            fvMesh& baseMesh,
+            const wordRes& zoneNames,
+            label exposedPatchId = -1
+        );
+
+        //- Construct from components. The subsetType is ZONES.
+        fvMeshSubsetProxy
+        (
+            fvMesh& baseMesh,
+            wordRes&& zoneNames,
+            label exposedPatchId = -1
         );
 
 
@@ -153,55 +174,90 @@ public:
 
     // Edit
 
-        //- Update mesh subset
+        //- Update of mesh subset
         void correct(bool verbose = false);
 
         //- Read mesh. Correct on topo-change
         polyMesh::readUpdateState readUpdate();
 
 
+    // Fields
+
         //- Construct volField (with zeroGradient) from an internal field
         template<class Type>
         static tmp<GeometricField<Type, fvPatchField, volMesh>>
         zeroGradientField
         (
-            const typename GeometricField
+            const DimensionedField<Type, volMesh>& df
-            <
+        );
-                Type,
+
-                fvPatchField,
+        //- Convert an internal field to a volume field (with zeroGradient)
-                volMesh
+        template<class Type>
-            >::Internal& fld
+        static tmp<GeometricField<Type, fvPatchField, volMesh>>
+        interpolateInternal
+        (
+            const fvMeshSubset& subsetter,
+            const DimensionedField<Type, volMesh>& df
+        );
+
+        //- Convert an internal field to a volume field (with zeroGradient)
+        //  Currently no proper memory reuse
+        template<class Type>
+        static tmp<GeometricField<Type, fvPatchField, volMesh>>
+        interpolateInternal
+        (
+            const fvMeshSubset& subsetter,
+            const tmp<DimensionedField<Type, volMesh>>& tdf
+        );
+
+        //- Wrapper for field or the subsetted field.
+        //  Pass through or forward to fvMeshSubset::interpolate()
+        template<class GeoField>
+        static tmp<GeoField> interpolate
+        (
+            const fvMeshSubset& subsetter,
+            const GeoField& fld
+        );
+
+        //- Wrapper for field or the subsetted field.
+        //  Pass through or forward to fvMeshSubset::interpolate()
+        template<class GeoField>
+        static tmp<GeoField> interpolate
+        (
+            const fvMeshSubset& subsetter,
+            const tmp<GeoField>& fld
         );
 
 
+    // Fields
+
+        //- Convert an internal field to a volume field (with zeroGradient)
+        template<class Type>
+        tmp<GeometricField<Type, fvPatchField, volMesh>>
+        interpolateInternal
+        (
+            const DimensionedField<Type, volMesh>& df
+        ) const;
+
+        //- Convert an internal field to a volume field (with zeroGradient)
+        //  Currently no proper memory reuse
+        template<class Type>
+        tmp<GeometricField<Type, fvPatchField, volMesh>>
+        interpolateInternal
+        (
+            const tmp<DimensionedField<Type, volMesh>>& tdf
+        ) const;
+
+
         //- Wrapper for field or the subsetted field.
-        //  Map volume field (does in fact do very little interpolation;
+        //  Pass through or forward to fvMeshSubset::interpolate()
-        //  just copied from fvMeshSubset)
-        template<class Type>
-        tmp<GeometricField<Type, fvPatchField, volMesh>>
-        interpolate
-        (
-            const GeometricField<Type, fvPatchField, volMesh>& fld
-        ) const;
-
-
-        //- Convert an internal field to a volume field
-        template<class Type>
-        tmp<GeometricField<Type, fvPatchField, volMesh>>
-        interpolate
-        (
-            const typename GeometricField
-            <
-                Type,
-                fvPatchField,
-                volMesh
-            >::Internal& fld
-        ) const;
-
-
-        //- Map volume field. Just forwards to fvMeshSubset
         template<class GeoField>
         tmp<GeoField> interpolate(const GeoField& fld) const;
+
+        //- Wrapper for field or the subsetted field.
+        //  Pass through or forward to fvMeshSubset::interpolate()
+        template<class GeoField>
+        tmp<GeoField> interpolate(const tmp<GeoField>& fld) const;
 };
 
 

src/dynamicMesh/fvMeshSubsetProxy/fvMeshSubsetProxyTemplates.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 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,18 +26,16 @@ License
 #include "fvMeshSubsetProxy.H"
 #include "volFields.H"
 
-// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
 template<class Type>
-Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
+Foam::tmp
+<
+    Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>
+>
 Foam::fvMeshSubsetProxy::zeroGradientField
 (
-    const typename GeometricField
+    const DimensionedField<Type, volMesh>& df
-    <
-        Type,
-        fvPatchField,
-        volMesh
-    >::Internal& df
 )
 {
     IOobject io(df);
@@ -62,44 +60,55 @@ Foam::fvMeshSubsetProxy::zeroGradientField
 
 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
-Foam::fvMeshSubsetProxy::interpolate
+Foam::fvMeshSubsetProxy::interpolateInternal
 (
-    const GeometricField<Type, fvPatchField, volMesh>& vf
+    const fvMeshSubset& subsetter,
-) const
+    const DimensionedField<Type, volMesh>& df
+)
 {
-    if (subsetter_.hasSubMesh())
+    auto tfield = zeroGradientField<Type>(df);
-    {
-        auto tfield(subsetter_.interpolate(vf));
 
-        tfield.ref().checkOut();
+    if (subsetter.hasSubMesh())
-        tfield.ref().rename(vf.name());
+    {
-        return tfield;
+        return interpolate(subsetter, tfield());
     }
 
-    return vf;
+    return tfield;
 }
 
 
 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
-Foam::fvMeshSubsetProxy::interpolate
+Foam::fvMeshSubsetProxy::interpolateInternal
 (
-    const typename GeometricField
+    const fvMeshSubset& subsetter,
-    <
+    const tmp<DimensionedField<Type, volMesh>>& tdf
-        Type,
+)
-        fvPatchField,
-        volMesh
-    >::Internal& df
-) const
 {
-    auto tfield = zeroGradientField<Type>(df);
+    // TODO - move dimensioned mesh into internal,
+    // but needs different GeometricField constructors
 
-    if (subsetter_.hasSubMesh())
+    if (tdf.valid())
     {
-        return interpolate<Type>(tfield());
+        if (subsetter.hasSubMesh())
+        {
+            auto tproxied = interpolate(subsetter, tdf);
+            auto tfield = zeroGradientField<Type>(tproxied());
+
+            tdf.clear();
+            tproxied.clear();
+            return tfield;
+        }
+        else
+        {
+            auto tfield = zeroGradientField<Type>(tdf());
+
+            tdf.clear();
+            return tfield;
+        }
     }
 
-    return tfield;
+    return nullptr;
 }
 
 
@@ -107,12 +116,14 @@ template<class GeoField>
 Foam::tmp<GeoField>
 Foam::fvMeshSubsetProxy::interpolate
 (
+    const fvMeshSubset& subsetter,
     const GeoField& fld
-) const
+)
 {
-    if (subsetter_.hasSubMesh())
+    if (subsetter.hasSubMesh())
     {
-        tmp<GeoField> tfield = subsetter_.interpolate(fld);
+        auto tfield = subsetter.interpolate(fld);
+
         tfield.ref().checkOut();
         tfield.ref().rename(fld.name());
         return tfield;
@@ -122,4 +133,65 @@ Foam::fvMeshSubsetProxy::interpolate
 }
 
 
+template<class GeoField>
+Foam::tmp<GeoField>
+Foam::fvMeshSubsetProxy::interpolate
+(
+    const fvMeshSubset& subsetter,
+    const tmp<GeoField>& tfield
+)
+{
+    if (tfield.valid() && subsetter.hasSubMesh())
+    {
+        auto tproxied = interpolate(subsetter, tfield());
+        tfield.clear();
+
+        return tproxied;
+    }
+
+    // Nothing to be done
+    return tfield;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
+Foam::fvMeshSubsetProxy::interpolateInternal
+(
+    const DimensionedField<Type, volMesh>& df
+) const
+{
+    return interpolateInternal(subsetter_, df);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
+Foam::fvMeshSubsetProxy::interpolateInternal
+(
+    const tmp<DimensionedField<Type, volMesh>>& tdf
+) const
+{
+    return interpolateInternal(subsetter_, tdf);
+}
+
+
+template<class GeoField>
+Foam::tmp<GeoField>
+Foam::fvMeshSubsetProxy::interpolate(const GeoField& fld) const
+{
+    return interpolate(subsetter_, fld);
+}
+
+
+template<class GeoField>
+Foam::tmp<GeoField>
+Foam::fvMeshSubsetProxy::interpolate(const tmp<GeoField>& tfield) const
+{
+    return interpolate(subsetter_, tfield);
+}
+
+
 // ************************************************************************* //