diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
index 8161e4ae15..c397cdcbc1 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
+++ b/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
-    fvMeshSubsetProxy meshRef(mesh, fvMeshSubsetProxy::ZONE, cellZoneName);
+    // Ignored (unproxied) if cellZoneName is empty
+    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
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
index faf3c7301f..e253be1f19 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
+++ b/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,
-                meshRef.baseMesh(),
+                meshProxy,
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 vScalarFld
@@ -728,8 +728,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
-                meshRef.baseMesh(),
+                meshProxy,
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 vVectorFld
@@ -738,8 +738,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
-                meshRef.baseMesh(),
+                meshProxy,
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 vSphTensorf
@@ -748,8 +748,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
-                meshRef.baseMesh(),
+                meshProxy,
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 vSymTensorFld
@@ -758,8 +758,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
-                meshRef.baseMesh(),
+                meshProxy,
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 vTensorFld
@@ -797,8 +797,8 @@ int main(int argc, char *argv[])
         {
             readFields
             (
-                meshRef,
-                meshRef.baseMesh(),
+                meshProxy,
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 dScalarFld
@@ -807,8 +807,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
-                meshRef.baseMesh(),
+                meshProxy,
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 dVectorFld
@@ -817,8 +817,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
-                meshRef.baseMesh(),
+                meshProxy,
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 dSphTensorFld
@@ -827,8 +827,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
-                meshRef.baseMesh(),
+                meshProxy,
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 dSymTensorFld
@@ -837,8 +837,8 @@ int main(int argc, char *argv[])
 
             readFields
             (
-                meshRef,
-                meshRef.baseMesh(),
+                meshProxy,
+                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,
-                meshRef.baseMesh(),
+                meshProxy,
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 sScalarFld
@@ -1158,8 +1158,8 @@ int main(int argc, char *argv[])
             PtrList<const surfaceVectorField> sVectorFld;
             readFields
             (
-                meshRef,
-                meshRef.baseMesh(),
+                meshProxy,
+                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,
-                meshRef.baseMesh(),
+                meshProxy,
+                meshProxy.baseMesh(),
                 objects,
                 selectedFields,
                 sScalarFld
@@ -1379,8 +1379,8 @@ int main(int argc, char *argv[])
             PtrList<const surfaceVectorField> sVectorFld;
             readFields
             (
-                meshRef,
-                meshRef.baseMesh(),
+                meshProxy,
+                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,
diff --git a/src/dynamicMesh/fvMeshSubset/fvMeshSubset.H b/src/dynamicMesh/fvMeshSubset/fvMeshSubset.H
index c6e83c361f..fd06ccee09 100644
--- a/src/dynamicMesh/fvMeshSubset/fvMeshSubset.H
+++ b/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
-        //  region == regioni.
+        //- Define cell subset, using the cells for which 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 labelList& cellMap,
-            const labelList& faceMap
+            const labelUList& patchMap,
+            const labelUList& cellMap,
+            const labelUList& faceMap
         );
 
         template<class Type>
@@ -361,9 +363,9 @@ public:
         (
             const GeometricField<Type, fvsPatchField, surfaceMesh>&,
             const fvMesh& sMesh,
-            const labelList& patchMap,
-            const labelList& cellMap,
-            const labelList& faceMap
+            const labelUList& patchMap,
+            const labelUList& cellMap,
+            const labelUList& faceMap
         );
 
         template<class Type>
@@ -380,8 +382,8 @@ public:
         (
             const GeometricField<Type, pointPatchField, pointMesh>&,
             const pointMesh& sMesh,
-            const labelList& patchMap,
-            const labelList& pointMap
+            const labelUList& patchMap,
+            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>
diff --git a/src/dynamicMesh/fvMeshSubset/fvMeshSubsetI.H b/src/dynamicMesh/fvMeshSubset/fvMeshSubsetI.H
index a58486f1e0..f19ecf6779 100644
--- a/src/dynamicMesh/fvMeshSubset/fvMeshSubsetI.H
+++ b/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();
diff --git a/src/dynamicMesh/fvMeshSubset/fvMeshSubsetInterpolate.C b/src/dynamicMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
index 0bb4b87cc2..c4920b330e 100644
--- a/src/dynamicMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
+++ b/src/dynamicMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
@@ -33,20 +33,18 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
 template<class Type>
-tmp<GeometricField<Type, fvPatchField, volMesh>>
-fvMeshSubset::interpolate
+Foam::tmp
+<
+    Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>
+>
+Foam::fvMeshSubset::interpolate
 (
     const GeometricField<Type, fvPatchField, volMesh>& vf,
     const fvMesh& sMesh,
-    const labelList& patchMap,
-    const labelList& cellMap,
-    const labelList& faceMap
+    const labelUList& patchMap,
+    const labelUList& cellMap,
+    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();
-    resF.oriented() = vf.oriented();
+    auto& result = tresult.ref();
+    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>>
-fvMeshSubset::interpolate
+Foam::tmp
+<
+    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>>
-fvMeshSubset::interpolate
+Foam::tmp
+<
+    Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>
+>
+Foam::fvMeshSubset::interpolate
 (
     const GeometricField<Type, fvsPatchField, surfaceMesh>& vf,
     const fvMesh& sMesh,
-    const labelList& patchMap,
-    const labelList& cellMap,
-    const labelList& faceMap
+    const labelUList& patchMap,
+    const labelUList& cellMap,
+    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();
-    resF.oriented() = vf.oriented();
+    auto& result = tresult.ref();
+    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>>
-fvMeshSubset::interpolate
+Foam::tmp
+<
+    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>>
-fvMeshSubset::interpolate
+Foam::tmp
+<
+    Foam::GeometricField<Type, Foam::pointPatchField, Foam::pointMesh>
+>
+Foam::fvMeshSubset::interpolate
 (
     const GeometricField<Type, pointPatchField, pointMesh>& vf,
     const pointMesh& sMesh,
-    const labelList& patchMap,
-    const labelList& pointMap
+    const labelUList& patchMap,
+    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();
-    resF.oriented() = vf.oriented();
+    auto& result = tresult.ref();
+    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>>
-fvMeshSubset::interpolate
+Foam::tmp
+<
+    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>>
-fvMeshSubset::interpolate
+Foam::tmp
+<
+    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 tresF = tmp<DimensionedField<Type, volMesh>>::New
+    auto tresult = 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>>
-fvMeshSubset::interpolate
+Foam::tmp
+<
+    Foam::DimensionedField<Type, Foam::volMesh>
+>
+Foam::fvMeshSubset::interpolate
 (
     const DimensionedField<Type, volMesh>& df
 ) const
@@ -535,8 +555,4 @@ fvMeshSubset::interpolate
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/dynamicMesh/fvMeshSubsetProxy/fvMeshSubsetProxy.C b/src/dynamicMesh/fvMeshSubsetProxy/fvMeshSubsetProxy.C
index b638746ed8..490cbceba3 100644
--- a/src/dynamicMesh/fvMeshSubsetProxy/fvMeshSubsetProxy.C
+++ b/src/dynamicMesh/fvMeshSubsetProxy/fvMeshSubsetProxy.C
@@ -30,18 +30,19 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::fvMeshSubsetProxy::fvMeshSubsetProxy
-(
-    fvMesh& baseMesh
-)
+Foam::fvMeshSubsetProxy::fvMeshSubsetProxy(fvMesh& baseMesh)
 :
     baseMesh_(baseMesh),
     subsetter_(baseMesh),
+    exposedPatchId_(-1),
     type_(NONE),
     name_(),
-    exposedPatchId_(-1)
+    names_()
 {
-    correct();
+    if (useSubMesh())
+    {
+        correct();
+    }
 }
 
 
@@ -49,17 +50,74 @@ Foam::fvMeshSubsetProxy::fvMeshSubsetProxy
 (
     fvMesh& baseMesh,
     const subsetType type,
-    const word& name,
-    const label exposedPatchId
+    const word& selectionName,
+    label exposedPatchId
 )
 :
     baseMesh_(baseMesh),
     subsetter_(baseMesh),
-    type_(name.empty() ? NONE : type),
-    name_(name),
-    exposedPatchId_(exposedPatchId)
+    exposedPatchId_(exposedPatchId),
+    type_(selectionName.empty() ? NONE : type),
+    name_(),
+    names_()
 {
-    correct();
+    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(baseMesh_, name_),
-            exposedPatchId_
-        );
+        cellSet cset(baseMesh_, name_);
+
+        selectedCells.resize(nCells);
+        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
-        (
-            baseMesh_.cellZones()[name_],
-            exposedPatchId_
-        );
+        selectedCells.resize(nCells);
+        selectedCells.set(baseMesh_.cellZones()[name_]);
     }
+    else if (type_ == ZONES)
+    {
+        if (verbose)
+        {
+            Info<< "Subsetting mesh based on cellZones "
+                << flatOutput(names_) << endl;
+        }
+
+        selectedCells = baseMesh_.cellZones().selection(names_);
+    }
+
+    subsetter_.setCellSubset(selectedCells, exposedPatchId_);
 }
 
 
diff --git a/src/dynamicMesh/fvMeshSubsetProxy/fvMeshSubsetProxy.H b/src/dynamicMesh/fvMeshSubsetProxy/fvMeshSubsetProxy.H
index b6b36916be..c29fde7cf8 100644
--- a/src/dynamicMesh/fvMeshSubsetProxy/fvMeshSubsetProxy.H
+++ b/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
-            SET,        //<! Using a cellSet for the subset
-            ZONE        //<! Using a cellZone for the subset
+            NONE,       //<! No subset
+            SET,        //<! Subset with a cellSet
+            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
-        fvMeshSubsetProxy(fvMesh& baseMesh);
+        //- Construct a pass-through proxy. No correct() invoked or required.
+        explicit fvMeshSubsetProxy(fvMesh& baseMesh);
 
         //- Construct from components
         fvMeshSubsetProxy
         (
             fvMesh& baseMesh,
-            const subsetType,
-            const word& name,
-            const label exposedPatchId = -1
+            const subsetType type,
+            const word& selectionName,
+            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
-            <
-                Type,
-                fvPatchField,
-                volMesh
-            >::Internal& fld
+            const DimensionedField<Type, volMesh>& df
+        );
+
+        //- Convert an internal field to a volume field (with zeroGradient)
+        template<class Type>
+        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;
-        //  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
+        //  Pass through or forward to fvMeshSubset::interpolate()
         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;
 };
 
 
diff --git a/src/dynamicMesh/fvMeshSubsetProxy/fvMeshSubsetProxyTemplates.C b/src/dynamicMesh/fvMeshSubsetProxy/fvMeshSubsetProxyTemplates.C
index 67dc2a7f95..85ad1e79f7 100644
--- a/src/dynamicMesh/fvMeshSubsetProxy/fvMeshSubsetProxyTemplates.C
+++ b/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
-    <
-        Type,
-        fvPatchField,
-        volMesh
-    >::Internal& df
+    const DimensionedField<Type, volMesh>& 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
+    const fvMeshSubset& subsetter,
+    const DimensionedField<Type, volMesh>& df
+)
 {
-    if (subsetter_.hasSubMesh())
-    {
-        auto tfield(subsetter_.interpolate(vf));
+    auto tfield = zeroGradientField<Type>(df);
 
-        tfield.ref().checkOut();
-        tfield.ref().rename(vf.name());
-        return tfield;
+    if (subsetter.hasSubMesh())
+    {
+        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
-    <
-        Type,
-        fvPatchField,
-        volMesh
-    >::Internal& df
-) const
+    const fvMeshSubset& subsetter,
+    const tmp<DimensionedField<Type, volMesh>>& tdf
+)
 {
-    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);
+}
+
+
 // ************************************************************************* //