From ffa939952c6ec02dcf007d4ef9c0b3d6c3cbcbe2 Mon Sep 17 00:00:00 2001
From: mattijs <m.janssens@opencfd.co.uk>
Date: Fri, 23 May 2008 11:06:43 +0100
Subject: [PATCH] subsetting point fields

---
 .../mesh/manipulation/subsetMesh/subsetMesh.C | 161 +++++++++++++++++-
 .../dataConversion/foamToVTK/foamToVTK.C      |  60 +++++--
 .../dataConversion/foamToVTK/readFields.C     |  44 -----
 .../dataConversion/foamToVTK/readFields.H     |  10 --
 .../dataConversion/foamToVTK/vtkMesh.C        |  15 +-
 .../dataConversion/foamToVTK/vtkMesh.H        |  44 +++--
 .../dataConversion/foamToVTK/vtkTopo.H        |  19 ++-
 .../basicSymmetryPointPatchField.C            |   2 +-
 .../calculated/calculatedPointPatchField.C    |   2 +-
 .../basic/coupled/coupledPointPatchField.C    |   2 +-
 .../basic/value/valuePointPatchField.C        |   2 +-
 .../zeroGradientPointPatchField.C             |   2 +-
 .../constraint/cyclic/cyclicPointPatchField.C |   4 +-
 .../constraint/empty/emptyPointPatchField.C   |   2 +-
 .../processor/processorPointPatchField.C      |   4 +-
 .../symmetry/symmetryPointPatchField.C        |   4 +-
 .../constraint/wedge/wedgePointPatchField.C   |   2 +-
 .../derived/global/globalPointPatchField.C    |   4 +-
 .../derived/slip/slipPointPatchField.C        |   4 +-
 .../fvMesh/fvMeshSubset/fvMeshSubset.C        | 111 +++++++-----
 .../fvMesh/fvMeshSubset/fvMeshSubset.H        |  82 ++++++++-
 .../fvMeshSubset/fvMeshSubsetInterpolate.C    | 142 ++++++++++++++-
 ...aceSlipDisplacementPointPatchVectorField.C |   2 +-
 23 files changed, 536 insertions(+), 188 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C b/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C
index 2a9bd022860..7c5e4baf778 100644
--- a/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C
+++ b/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C
@@ -60,7 +60,7 @@ void subsetVolFields
 
         Info<< "Subsetting field " << fieldName << endl;
 
-        GeometricField<Type, fvPatchField, volMesh> volField    
+        GeometricField<Type, fvPatchField, volMesh> fld    
         (
             IOobject
             (
@@ -73,7 +73,7 @@ void subsetVolFields
             baseMesh
         );
 
-        subFields.set(i, subsetter.interpolate(volField));
+        subFields.set(i, subsetter.interpolate(fld));
     }
 }
 
@@ -94,7 +94,7 @@ void subsetSurfaceFields
 
         Info<< "Subsetting field " << fieldName << endl;
 
-        GeometricField<Type, fvsPatchField, surfaceMesh> volField    
+        GeometricField<Type, fvsPatchField, surfaceMesh> fld    
         (
             IOobject
             (
@@ -107,7 +107,42 @@ void subsetSurfaceFields
             baseMesh
         );
 
-        subFields.set(i, subsetter.interpolate(volField));
+        subFields.set(i, subsetter.interpolate(fld));
+    }
+}
+
+
+template<class Type>
+void subsetPointFields
+(
+    const fvMeshSubset& subsetter,
+    const pointMesh& pMesh,
+    const wordList& fieldNames,
+    PtrList<GeometricField<Type, pointPatchField, pointMesh> >& subFields
+)
+{
+    const fvMesh& baseMesh = subsetter.baseMesh();
+
+    forAll(fieldNames, i)
+    {
+        const word& fieldName = fieldNames[i];
+
+        Info<< "Subsetting field " << fieldName << endl;
+
+        GeometricField<Type, pointPatchField, pointMesh> fld    
+        (
+            IOobject
+            (
+                fieldName,
+                baseMesh.time().timeName(),
+                baseMesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            ),
+            pMesh
+        );
+
+        subFields.set(i, subsetter.interpolate(fld));
     }
 }
 
@@ -163,7 +198,9 @@ int main(int argc, char *argv[])
 
     IOobjectList objects(mesh, runTime.timeName());
 
-    // Read vol fields and subset.
+
+    // Read vol fields and subset
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     wordList scalarNames(objects.names(volScalarField::typeName));
     PtrList<volScalarField> scalarFlds(scalarNames.size());
@@ -191,7 +228,9 @@ int main(int argc, char *argv[])
     PtrList<volTensorField> tensorFlds(tensorNames.size());
     subsetVolFields(subsetter, tensorNames, tensorFlds);
 
-    // Read surface fields and subset.
+
+    // Read surface fields and subset
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     wordList surfScalarNames(objects.names(surfaceScalarField::typeName));
     PtrList<surfaceScalarField> surfScalarFlds(surfScalarNames.size());
@@ -231,6 +270,59 @@ int main(int argc, char *argv[])
     subsetSurfaceFields(subsetter, surfTensorNames, surfTensorFlds);
 
 
+    // Read point fields and subset
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    pointMesh pMesh(mesh);
+
+    wordList pointScalarNames(objects.names(pointScalarField::typeName));
+    PtrList<pointScalarField> pointScalarFlds(pointScalarNames.size());
+    subsetPointFields(subsetter, pMesh, pointScalarNames, pointScalarFlds);
+
+    wordList pointVectorNames(objects.names(pointVectorField::typeName));
+    PtrList<pointVectorField> pointVectorFlds(pointVectorNames.size());
+    subsetPointFields(subsetter, pMesh, pointVectorNames, pointVectorFlds);
+
+    wordList pointSphericalTensorNames
+    (
+        objects.names(pointSphericalTensorField::typeName)
+    );
+    PtrList<pointSphericalTensorField> pointSphericalTensorFlds
+    (
+        pointSphericalTensorNames.size()
+    );
+    subsetPointFields
+    (
+        subsetter,
+        pMesh,
+        pointSphericalTensorNames,
+        pointSphericalTensorFlds
+    );
+
+    wordList pointSymmTensorNames
+    (
+        objects.names(pointSymmTensorField::typeName)
+    );
+    PtrList<pointSymmTensorField> pointSymmTensorFlds
+    (
+        pointSymmTensorNames.size()
+    );
+    subsetPointFields
+    (
+        subsetter,
+        pMesh,
+        pointSymmTensorNames,
+        pointSymmTensorFlds
+    );
+
+    wordList pointTensorNames(objects.names(pointTensorField::typeName));
+    PtrList<pointTensorField> pointTensorFlds(pointTensorNames.size());
+    subsetPointFields(subsetter, pMesh, pointTensorNames, pointTensorFlds);
+
+
+
+    // Write mesh and fields to new time
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     runTime++;
 
@@ -252,6 +344,18 @@ int main(int argc, char *argv[])
 
         vectorFlds[i].write();
     }
+    forAll(sphericalTensorFlds, i)
+    {
+        sphericalTensorFlds[i].rename(sphericalTensorNames[i]);
+
+        sphericalTensorFlds[i].write();
+    }
+    forAll(symmTensorFlds, i)
+    {
+        symmTensorFlds[i].rename(symmTensorNames[i]);
+
+        symmTensorFlds[i].write();
+    }
     forAll(tensorFlds, i)
     {
         tensorFlds[i].rename(tensorNames[i]);
@@ -272,6 +376,18 @@ int main(int argc, char *argv[])
 
         surfVectorFlds[i].write();
     }
+    forAll(surfSphericalTensorFlds, i)
+    {
+        surfSphericalTensorFlds[i].rename(surfSphericalTensorNames[i]);
+
+        surfSphericalTensorFlds[i].write();
+    }
+    forAll(surfSymmTensorFlds, i)
+    {
+        surfSymmTensorFlds[i].rename(surfSymmTensorNames[i]);
+
+        surfSymmTensorFlds[i].write();
+    }
     forAll(surfTensorNames, i)
     {
         surfTensorFlds[i].rename(surfTensorNames[i]);
@@ -279,6 +395,39 @@ int main(int argc, char *argv[])
         surfTensorFlds[i].write();
     }
 
+    // Point ones
+    forAll(pointScalarFlds, i)
+    {
+        pointScalarFlds[i].rename(pointScalarNames[i]);
+
+        pointScalarFlds[i].write();
+    }
+    forAll(pointVectorFlds, i)
+    {
+        pointVectorFlds[i].rename(pointVectorNames[i]);
+
+        pointVectorFlds[i].write();
+    }
+    forAll(pointSphericalTensorFlds, i)
+    {
+        pointSphericalTensorFlds[i].rename(pointSphericalTensorNames[i]);
+
+        pointSphericalTensorFlds[i].write();
+    }
+    forAll(pointSymmTensorFlds, i)
+    {
+        pointSymmTensorFlds[i].rename(pointSymmTensorNames[i]);
+
+        pointSymmTensorFlds[i].write();
+    }
+    forAll(pointTensorNames, i)
+    {
+        pointTensorFlds[i].rename(pointTensorNames[i]);
+
+        pointTensorFlds[i].write();
+    }
+
+
     Info << nl << "End" << endl;
 
     return 0;
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
index 85efbfaa103..c9d7409a02a 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
@@ -469,7 +469,6 @@ int main(int argc, char *argv[])
 
         // Construct pointMesh only if nessecary since constructs edge
         // addressing (expensive on polyhedral meshes)
-        autoPtr<pointMesh> pMeshPtr(NULL);
         if (noPointValues)
         {
             Info<< "    pointScalarFields : switched off"
@@ -477,10 +476,6 @@ int main(int argc, char *argv[])
             Info<< "    pointVectorFields : switched off"
                 << " (\"-noPointValues\" option)\n";
         }
-        else
-        {
-            pMeshPtr.reset(new pointMesh(mesh));
-        }
 
         PtrList<pointScalarField> psf;
         PtrList<pointVectorField> pvf;
@@ -488,21 +483,56 @@ int main(int argc, char *argv[])
         PtrList<pointSymmTensorField> pSymmtf;
         PtrList<pointTensorField> ptf;
 
-        if (pMeshPtr.valid())
+        if (!noPointValues)
         {
-            readFields(pMeshPtr(), objects, selectedFields, psf);
+            readFields
+            (
+                vMesh,
+                vMesh.basePointMesh(),
+                objects,
+                selectedFields,
+                psf
+            );
             print("    pointScalarFields          :", Info, psf);
 
-            readFields(pMeshPtr(), objects, selectedFields, pvf);
+            readFields
+            (
+                vMesh,
+                vMesh.basePointMesh(),
+                objects,
+                selectedFields,
+                pvf
+            );
             print("    pointVectorFields          :", Info, pvf);
 
-            readFields(pMeshPtr(), objects, selectedFields, pSpheretf);
+            readFields
+            (
+                vMesh,
+                vMesh.basePointMesh(),
+                objects,
+                selectedFields,
+                pSpheretf
+            );
             print("    pointSphericalTensorFields :", Info, pSpheretf);
 
-            readFields(pMeshPtr(), objects, selectedFields, pSymmtf);
+            readFields
+            (
+                vMesh,
+                vMesh.basePointMesh(),
+                objects,
+                selectedFields,
+                pSymmtf
+            );
             print("    pointSymmTensorFields      :", Info, pSymmtf);
 
-            readFields(pMeshPtr(), objects, selectedFields, ptf);
+            readFields
+            (
+                vMesh,
+                vMesh.basePointMesh(),
+                objects,
+                selectedFields,
+                ptf
+            );
             print("    pointTensorFields          :", Info, ptf);
         }
         Info<< endl;
@@ -550,7 +580,7 @@ int main(int argc, char *argv[])
             writer.write(vSymmtf);
             writer.write(vtf);
 
-            if (pMeshPtr.valid())
+            if (!noPointValues)
             {
                 writeFuns::writePointDataHeader
                 (
@@ -567,7 +597,7 @@ int main(int argc, char *argv[])
                 writer.write(ptf);
 
                 // Interpolated volFields
-                volPointInterpolation pInterp(mesh, pMeshPtr());
+                volPointInterpolation pInterp(mesh, vMesh.pMesh());
                 writer.write(pInterp, vsf);
                 writer.write(pInterp, vvf);
                 writer.write(pInterp, vSpheretf);
@@ -705,7 +735,7 @@ int main(int argc, char *argv[])
             writer.write(vSymmtf);
             writer.write(vtf);
 
-            if (pMeshPtr.valid())
+            if (!noPointValues)
             {
                 writeFuns::writePointDataHeader
                 (
@@ -785,7 +815,7 @@ int main(int argc, char *argv[])
                         writer.write(vSymmtf);
                         writer.write(vtf);
 
-                        if (pMeshPtr.valid())
+                        if (!noPointValues)
                         {
                             writeFuns::writePointDataHeader
                             (
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.C
index 60d09309c30..975124046f7 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.C
@@ -34,50 +34,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
-//template<class GeoField, class Mesh>
-template<class GeoField>
-void readFields
-(
-    const typename GeoField::Mesh& mesh,
-    const IOobjectList& objects,
-    const HashSet<word>& selectedFields,
-    PtrList<GeoField>& fields
-)
-{
-    // Search list of objects for volScalarFields
-    IOobjectList fieldObjects(objects.lookupClass(GeoField::typeName));
-
-    // Construct the vol scalar fields
-    fields.setSize(fieldObjects.size());
-    label nFields = 0;
-
-    for
-    (
-        IOobjectList::iterator iter = fieldObjects.begin();
-        iter != fieldObjects.end();
-        ++iter
-    )
-    {
-        if (!selectedFields.size() || selectedFields.found(iter()->name()))
-        {
-            fields.set
-            (
-                nFields,
-                new GeoField
-                (
-                    *iter(),
-                    mesh
-                )
-            );
-
-            nFields++;
-        }
-    }
-
-    fields.setSize(nFields);
-}
-
-
 template<class GeoField>
 void readFields
 (
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.H
index 577bdf32f0a..db98dd948e9 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.H
@@ -45,16 +45,6 @@ SourceFiles
 namespace Foam
 {
 
-// Read the fields and put on the pointer list
-template<class GeoField>
-void readFields
-(
-    const typename GeoField::Mesh& mesh,
-    const IOobjectList& objects,
-    const HashSet<word>& selectedFields,
-    PtrList<GeoField>& fields
-);
-
 // Read the fields and optionally subset and put on the pointer list
 template<class GeoField>
 void readFields
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.C
index c96cfd4f676..aecd725ddc1 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.C
@@ -45,8 +45,7 @@ Foam::vtkMesh::vtkMesh
 :
     baseMesh_(baseMesh),
     subsetter_(baseMesh),
-    setName_(setName),
-    topoPtr_(NULL)
+    setName_(setName)
 {
     if (setName.size() > 0)
     {
@@ -59,14 +58,6 @@ Foam::vtkMesh::vtkMesh
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::vtkMesh::~vtkMesh()
-{
-    deleteDemandDrivenData(topoPtr_);
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::polyMesh::readUpdateState Foam::vtkMesh::readUpdate()
@@ -78,11 +69,11 @@ Foam::polyMesh::readUpdateState Foam::vtkMesh::readUpdate()
         // Note: since fvMeshSubset has no movePoints() functionality reconstruct
         // the subset even if only movement.
 
-        deleteDemandDrivenData(topoPtr_);
+        topoPtr_.clear();
 
         if (setName_.size() > 0)
         {
-            Pout<< "Subsetting mesh based on cellSet " << setName_ << endl;
+            Info<< "Subsetting mesh based on cellSet " << setName_ << endl;
 
             // Read cellSet using whole mesh
             cellSet currentSet(baseMesh_, setName_);
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.H
index b1ceec83dda..8d4ad36ab47 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.H
@@ -39,6 +39,7 @@ SourceFiles
 
 #include "vtkTopo.H"
 #include "fvMeshSubset.H"
+#include "pointMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,14 +60,17 @@ class vtkMesh
         //- Reference to mesh
         fvMesh& baseMesh_;
 
-        //- Subsetting engine
+        //- Demand driven pointMesh
+        mutable autoPtr<pointMesh> pointMeshPtr_;
+
+        //- Subsetting engine + sub-fvMesh
         fvMeshSubset subsetter_;
 
         //- Current cellSet (or empty)
         const word setName_;
 
         //- Current decomposition of topology
-        mutable vtkTopo* topoPtr_;
+        mutable autoPtr<vtkTopo> topoPtr_;
 
 
 
@@ -87,11 +91,6 @@ public:
         vtkMesh(fvMesh& baseMesh, const word& setName = "");
 
 
-    // Destructor
-
-        ~vtkMesh();
-
-
     // Member Functions
 
         // Access
@@ -102,6 +101,15 @@ public:
                 return baseMesh_;
             }
 
+            const pointMesh& basePointMesh() const
+            {
+                if (!pointMeshPtr_.valid())
+                {
+                    pointMeshPtr_.reset(new pointMesh(baseMesh_));
+                }
+                return pointMeshPtr_();
+            }
+
             const fvMeshSubset& subsetter() const
             {
                 return subsetter_;
@@ -116,11 +124,11 @@ public:
             //- topology
             const vtkTopo& topo() const
             {
-                if (!topoPtr_)
+                if (!topoPtr_.valid())
                 {
-                    topoPtr_ = new vtkTopo(mesh());
+                    topoPtr_.reset(new vtkTopo(mesh()));
                 }
-                return *topoPtr_;
+                return topoPtr_();
             }
 
             //- Access either mesh or submesh
@@ -136,6 +144,19 @@ public:
                 }
             }
 
+            //- Access either pointMesh of base or pointMesh of subset
+            const pointMesh& pMesh() const
+            {
+                if (useSubMesh())
+                {
+                    return subsetter_.subPointMesh();
+                }
+                else
+                {
+                    return basePointMesh();
+                }
+            }
+
             //- Number of field cells
             label nFieldCells() const
             {
@@ -145,7 +166,7 @@ public:
             //- Number of field points
             label nFieldPoints() const
             {
-                return mesh().nPoints() + topo().addPointCellLabels().size();
+                return pMesh().size() + topo().addPointCellLabels().size();
             }
 
 
@@ -171,7 +192,6 @@ public:
                     return fld;
                 }
             }
-
 };
 
 
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkTopo.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkTopo.H
index b45a5f3d2e5..6867eac71c0 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkTopo.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkTopo.H
@@ -81,14 +81,17 @@ public:
    // Public static data
 
         // this must be consistent with the enumeration in "vtkCell.H"
-        static const label VTK_TRIANGLE   = 5;
-        static const label VTK_POLYGON    = 7;
-        static const label VTK_QUAD       = 9;
-
-        static const label VTK_TETRA      = 10;
-        static const label VTK_PYRAMID    = 14;
-        static const label VTK_WEDGE      = 13;
-        static const label VTK_HEXAHEDRON = 12;
+        enum vtkTypes
+        {
+            VTK_TRIANGLE   = 5,
+            VTK_POLYGON    = 7,
+            VTK_QUAD       = 9,
+
+            VTK_TETRA      = 10,
+            VTK_PYRAMID    = 14,
+            VTK_WEDGE      = 13,
+            VTK_HEXAHEDRON = 12,
+        };
 
     // Constructors
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C
index 1dec091e8a2..fe04cb04548 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C
@@ -66,7 +66,7 @@ basicSymmetryPointPatchField<Type>::basicSymmetryPointPatchField
     const pointPatchFieldMapper&
 )
 :
-    pointPatchField<Type>(ptf, iF)
+    pointPatchField<Type>(p, iF)
 {}
 
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C
index f3367bd57ee..b57ae68f80b 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C
@@ -72,7 +72,7 @@ calculatedPointPatchField<Type>::calculatedPointPatchField
     const pointPatchFieldMapper&
 )
 :
-    pointPatchField<Type>(ptf, iF)
+    pointPatchField<Type>(p, iF)
 {}
 
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C
index 80fa72bdeba..6d2a1250d14 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C
@@ -65,7 +65,7 @@ coupledPointPatchField<Type>::coupledPointPatchField
     const pointPatchFieldMapper&
 )
 :
-    pointPatchField<Type>(ptf, iF)
+    pointPatchField<Type>(p, iF)
 {}
 
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C
index 261c42e5d6a..5d95ea2964f 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C
@@ -114,7 +114,7 @@ valuePointPatchField<Type>::valuePointPatchField
     const pointPatchFieldMapper& mapper
 )
 :
-    pointPatchField<Type>(ptf, iF),
+    pointPatchField<Type>(p, iF),
     Field<Type>(ptf, mapper)
 {}
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C
index e71bcfe693d..7e84b13134c 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C
@@ -65,7 +65,7 @@ zeroGradientPointPatchField<Type>::zeroGradientPointPatchField
     const pointPatchFieldMapper&
 )
 :
-    pointPatchField<Type>(ptf, iF)
+    pointPatchField<Type>(p, iF)
 {}
 
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/cyclic/cyclicPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/cyclic/cyclicPointPatchField.C
index fc1dcc205ad..32a68dd977d 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/cyclic/cyclicPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/cyclic/cyclicPointPatchField.C
@@ -81,10 +81,10 @@ cyclicPointPatchField<Type>::cyclicPointPatchField
     const cyclicPointPatchField<Type>& ptf,
     const pointPatch& p,
     const DimensionedField<Type, pointMesh>& iF,
-    const pointPatchFieldMapper&
+    const pointPatchFieldMapper& mapper
 )
 :
-    coupledPointPatchField<Type>(ptf, iF),
+    coupledPointPatchField<Type>(ptf, p, iF, mapper),
     cyclicPatch_(refCast<const cyclicPointPatch>(p))
 {
     if (!isType<cyclicPointPatch>(this->patch()))
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C
index a99e9255a71..e467baf9415 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C
@@ -81,7 +81,7 @@ emptyPointPatchField<Type>::emptyPointPatchField
     const pointPatchFieldMapper&
 )
 :
-    pointPatchField<Type>(ptf, iF)
+    pointPatchField<Type>(p, iF)
 {
     if (!isType<emptyPointPatch>(this->patch()))
     {
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.C
index 525e9b41790..831b023b972 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.C
@@ -66,10 +66,10 @@ processorPointPatchField<Type>::processorPointPatchField
     const processorPointPatchField<Type>& ptf,
     const pointPatch& p,
     const DimensionedField<Type, pointMesh>& iF,
-    const pointPatchFieldMapper&
+    const pointPatchFieldMapper& mapper
 )
 :
-    coupledPointPatchField<Type>(ptf, iF),
+    coupledPointPatchField<Type>(ptf, p, iF, mapper),
     procPatch_(refCast<const processorPointPatch>(ptf.patch()))
 {}
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/symmetry/symmetryPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/symmetry/symmetryPointPatchField.C
index 885a06a7c7d..62eb09b5968 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/symmetry/symmetryPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/symmetry/symmetryPointPatchField.C
@@ -78,10 +78,10 @@ symmetryPointPatchField<Type>::symmetryPointPatchField
     const symmetryPointPatchField<Type>& ptf,
     const pointPatch& p,
     const DimensionedField<Type, pointMesh>& iF,
-    const pointPatchFieldMapper&
+    const pointPatchFieldMapper& mapper
 )
 :
-    basicSymmetryPointPatchField<Type>(ptf, iF)
+    basicSymmetryPointPatchField<Type>(ptf, p, iF, mapper)
 {
     if (!isType<symmetryPointPatch>(this->patch()))
     {
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C
index 33f136da203..2c0b4f721fe 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C
@@ -82,7 +82,7 @@ wedgePointPatchField<Type>::wedgePointPatchField
     const pointPatchFieldMapper&
 )
 :
-    pointPatchField<Type>(ptf, iF)
+    pointPatchField<Type>(p, iF)
 {
     if (!isType<wedgePointPatch>(this->patch()))
     {
diff --git a/src/OpenFOAM/fields/pointPatchFields/derived/global/globalPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/derived/global/globalPointPatchField.C
index eed5b1d0f4e..2dfe6378f91 100644
--- a/src/OpenFOAM/fields/pointPatchFields/derived/global/globalPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/derived/global/globalPointPatchField.C
@@ -82,10 +82,10 @@ globalPointPatchField<Type>::globalPointPatchField
     const globalPointPatchField<Type>& ptf,
     const pointPatch& p,
     const DimensionedField<Type, pointMesh>& iF,
-    const pointPatchFieldMapper&
+    const pointPatchFieldMapper& mapper
 )
 :
-    coupledPointPatchField<Type>(ptf, iF),
+    coupledPointPatchField<Type>(ptf, p, iF, mapper),
     globalPointPatch_(refCast<const globalPointPatch>(ptf.patch()))
 {
     if (!isType<globalPointPatch>(this->patch()))
diff --git a/src/OpenFOAM/fields/pointPatchFields/derived/slip/slipPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/derived/slip/slipPointPatchField.C
index ae1014dcea0..b32cf3e5fe1 100644
--- a/src/OpenFOAM/fields/pointPatchFields/derived/slip/slipPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/derived/slip/slipPointPatchField.C
@@ -62,10 +62,10 @@ slipPointPatchField<Type>::slipPointPatchField
     const slipPointPatchField<Type>& ptf,
     const pointPatch& p,
     const DimensionedField<Type, pointMesh>& iF,
-    const pointPatchFieldMapper&
+    const pointPatchFieldMapper& mapper
 )
 :
-    basicSymmetryPointPatchField<Type>(ptf, iF)
+    basicSymmetryPointPatchField<Type>(ptf, p, iF, mapper)
 {}
 
 
diff --git a/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.C b/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.C
index 26fb6ebbe97..9947aabbe8f 100644
--- a/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.C
+++ b/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.C
@@ -48,7 +48,7 @@ namespace Foam
 
 bool Foam::fvMeshSubset::checkCellSubset() const
 {
-    if (!fvMeshSubsetPtr_)
+    if (!fvMeshSubsetPtr_.valid())
     {
         FatalErrorIn("bool fvMeshSubset::checkCellSubset() const")
             << "Mesh subset not set.  Please set the cell map using "
@@ -267,7 +267,7 @@ void Foam::fvMeshSubset::subsetZones()
             pz.name(),
             subset(baseMesh().nPoints(), pz, pointMap()),
             i,
-            fvMeshSubsetPtr_->pointZones()
+            fvMeshSubsetPtr_().pointZones()
         );
     }
 
@@ -315,7 +315,7 @@ void Foam::fvMeshSubset::subsetZones()
             subAddressing,
             subFlipStatus,
             i,
-            fvMeshSubsetPtr_->faceZones()
+            fvMeshSubsetPtr_().faceZones()
         );
     }
 
@@ -333,13 +333,13 @@ void Foam::fvMeshSubset::subsetZones()
             cz.name(),
             subset(baseMesh().nCells(), cz, cellMap()),
             i,
-            fvMeshSubsetPtr_->cellZones()
+            fvMeshSubsetPtr_().cellZones()
         );
     }
 
 
     // Add the zones
-    fvMeshSubsetPtr_->addZones(pZonePtrs, fZonePtrs, cZonePtrs);
+    fvMeshSubsetPtr_().addZones(pZonePtrs, fZonePtrs, cZonePtrs);
 }
 
 
@@ -357,14 +357,6 @@ Foam::fvMeshSubset::fvMeshSubset(const fvMesh& baseMesh)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::fvMeshSubset::~fvMeshSubset()
-{
-    deleteDemandDrivenData(fvMeshSubsetPtr_);
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::fvMeshSubset::setCellSubset
@@ -671,20 +663,24 @@ void Foam::fvMeshSubset::setCellSubset
 
 
     // Make a new mesh
-    fvMeshSubsetPtr_ = new fvMesh
+    fvMeshSubsetPtr_.reset
     (
-        IOobject
+        new fvMesh
         (
-            baseMesh().name() + "SubSet",
-            baseMesh().time().timeName(),
-            baseMesh().time(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        newPoints,
-        newFaces,
-        newCells
+            IOobject
+            (
+                baseMesh().name() + "SubSet",
+                baseMesh().time().timeName(),
+                baseMesh().time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            newPoints,
+            newFaces,
+            newCells
+        )
     );
+    pointMeshSubsetPtr_.clear();
 
 
     // Add old patches
@@ -700,7 +696,7 @@ void Foam::fvMeshSubset::setCellSubset
             // Patch still exists. Add it
             newBoundary[nNewPatches] = oldPatches[patchI].clone
             (
-                fvMeshSubsetPtr_->boundaryMesh(),
+                fvMeshSubsetPtr_().boundaryMesh(),
                 nNewPatches,
                 boundaryPatchSizes[patchI],
                 patchStart
@@ -723,7 +719,7 @@ void Foam::fvMeshSubset::setCellSubset
                 boundaryPatchSizes[oldInternalPatchID],
                 patchStart,
                 nNewPatches,
-                fvMeshSubsetPtr_->boundaryMesh()
+                fvMeshSubsetPtr_().boundaryMesh()
             );
 
             // The index for the first patch is -1 as it originates from
@@ -738,7 +734,7 @@ void Foam::fvMeshSubset::setCellSubset
     patchMap_.setSize(nNewPatches);
 
     // Add the fvPatches
-    fvMeshSubsetPtr_->addFvPatches(newBoundary);
+    fvMeshSubsetPtr_().addFvPatches(newBoundary);
 
     // Subset and add any zones
     subsetZones();
@@ -1166,22 +1162,25 @@ void Foam::fvMeshSubset::setLargeCellSubset
     // not proper but cannot be avoided since otherwise surfaceInterpolation
     // cannot find its fvSchemes (it will try to read e.g.
     // system/region0SubSet/fvSchemes)
-    fvMeshSubsetPtr_ = new fvMesh
+    fvMeshSubsetPtr_.reset
     (
-        IOobject
+        new fvMesh
         (
-            baseMesh().name(),
-            baseMesh().time().timeName(),
-            baseMesh().time(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        newPoints,
-        newFaces,
-        newCells,
-        syncPar           // parallel synchronisation
+            IOobject
+            (
+                baseMesh().name(),
+                baseMesh().time().timeName(),
+                baseMesh().time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            newPoints,
+            newFaces,
+            newCells,
+            syncPar           // parallel synchronisation
+        )
     );
-
+    pointMeshSubsetPtr_.clear();
 
     // Add old patches
     List<polyPatch*> newBoundary(nbSize);
@@ -1251,7 +1250,7 @@ void Foam::fvMeshSubset::setLargeCellSubset
         // Clone (even if 0 size)
         newBoundary[nNewPatches] = oldPatches[oldPatchI].clone
         (
-            fvMeshSubsetPtr_->boundaryMesh(),
+            fvMeshSubsetPtr_().boundaryMesh(),
             nNewPatches,
             newSize,
             patchStart
@@ -1282,7 +1281,7 @@ void Foam::fvMeshSubset::setLargeCellSubset
                 boundaryPatchSizes[oldInternalPatchID],
                 patchStart,
                 nNewPatches,
-                fvMeshSubsetPtr_->boundaryMesh()
+                fvMeshSubsetPtr_().boundaryMesh()
             );
 
             //Pout<< "    oldInternalFaces : "
@@ -1310,7 +1309,7 @@ void Foam::fvMeshSubset::setLargeCellSubset
         // Patch still exists. Add it
         newBoundary[nNewPatches] = oldPatches[oldPatchI].clone
         (
-            fvMeshSubsetPtr_->boundaryMesh(),
+            fvMeshSubsetPtr_().boundaryMesh(),
             nNewPatches,
             newSize,
             patchStart
@@ -1331,7 +1330,7 @@ void Foam::fvMeshSubset::setLargeCellSubset
 
 
     // Add the fvPatches
-    fvMeshSubsetPtr_->addFvPatches(newBoundary);
+    fvMeshSubsetPtr_().addFvPatches(newBoundary);
 
     // Subset and add any zones
     subsetZones();
@@ -1359,7 +1358,7 @@ const fvMesh& Foam::fvMeshSubset::subMesh() const
 {
     checkCellSubset();
 
-    return *fvMeshSubsetPtr_;
+    return fvMeshSubsetPtr_();
 }
 
 
@@ -1367,7 +1366,27 @@ fvMesh& Foam::fvMeshSubset::subMesh()
 {
     checkCellSubset();
 
-    return *fvMeshSubsetPtr_;
+    return fvMeshSubsetPtr_();
+}
+
+
+const pointMesh& Foam::fvMeshSubset::subPointMesh() const
+{
+    if (!pointMeshSubsetPtr_.valid())
+    {
+        pointMeshSubsetPtr_.reset(new pointMesh(subMesh()));
+    }
+    return pointMeshSubsetPtr_();
+}
+
+
+pointMesh& Foam::fvMeshSubset::subPointMesh()
+{
+    if (!pointMeshSubsetPtr_.valid())
+    {
+        pointMeshSubsetPtr_.reset(new pointMesh(subMesh()));
+    }
+    return pointMeshSubsetPtr_();
 }
 
 
diff --git a/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.H b/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.H
index d5c55900dab..e8428aea07d 100644
--- a/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.H
+++ b/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.H
@@ -56,11 +56,11 @@ SourceFiles
 #define fvMeshSubset_H
 
 #include "fvMesh.H"
+#include "pointMesh.H"
 #include "fvPatchFieldMapper.H"
+#include "pointPatchFieldMapper.H"
 #include "GeometricField.H"
-#include "emptyFvPatchFields.H"
 #include "labelHashSet.H"
-#include "SubField.H"
 #include "surfaceMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -119,6 +119,48 @@ public:
     };
 
 
+    //- Patch-field subset interpolation class
+    class pointPatchFieldSubset
+    :
+        public pointPatchFieldMapper
+    {
+        const labelList& directAddressing_;
+
+    public:
+
+        // Constructors
+
+            //- Construct given addressing
+            pointPatchFieldSubset(const labelList& directAddressing)
+            :
+                directAddressing_(directAddressing)
+            {}
+
+        // Destructor
+
+            virtual ~pointPatchFieldSubset()
+            {}
+
+
+        // Member Functions
+
+            label size() const
+            {
+                return directAddressing_.size();
+            }
+
+            bool direct() const
+            {
+                return true;
+            }
+
+            const unallocLabelList& directAddressing() const
+            {
+                return directAddressing_;
+            }
+    };
+
+
 private:
 
     // Private data
@@ -127,7 +169,9 @@ private:
         const fvMesh& baseMesh_;
 
         //- Subset mesh pointer
-        fvMesh* fvMeshSubsetPtr_;
+        autoPtr<fvMesh> fvMeshSubsetPtr_;
+
+        mutable autoPtr<pointMesh> pointMeshSubsetPtr_;
 
         //- Point mapping array
         labelList pointMap_;
@@ -185,11 +229,6 @@ public:
         explicit fvMeshSubset(const fvMesh&);
 
 
-    // Destructor
-
-        ~fvMeshSubset();
-
-
     // Member Functions
 
         // Edit
@@ -236,8 +275,14 @@ public:
 
             //- Return reference to subset mesh
             const fvMesh& subMesh() const;
+
             fvMesh& subMesh();
 
+            //- Return reference to demand-driven subset pointMesh
+            const pointMesh& subPointMesh() const;
+
+            pointMesh& subPointMesh();
+
             //- Return point map
             const labelList& pointMap() const;
 
@@ -275,7 +320,7 @@ public:
             //- Map surface field
             template<class Type>
             static tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
-             interpolate
+            interpolate
             (
                 const GeometricField<Type, fvsPatchField, surfaceMesh>&,
                 const fvMesh& sMesh,
@@ -289,6 +334,25 @@ public:
             (
                 const GeometricField<Type, fvsPatchField, surfaceMesh>&
             ) const;
+
+            //- Map point field
+            template<class Type>
+            static tmp<GeometricField<Type, pointPatchField, pointMesh> >
+            interpolate
+            (
+                const GeometricField<Type, pointPatchField, pointMesh>&,
+                const pointMesh& sMesh,
+                const objectRegistry& reg,
+                const labelList& patchMap,
+                const labelList& pointMap
+            );
+
+            template<class Type>
+            tmp<GeometricField<Type, pointPatchField, pointMesh> >
+            interpolate
+            (
+                const GeometricField<Type, pointPatchField, pointMesh>&
+            ) const;
 };
 
 
diff --git a/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C b/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
index 40045e89c01..3b59b7b6b9f 100644
--- a/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
+++ b/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
@@ -26,6 +26,8 @@ License
 
 #include "fvMeshSubset.H"
 #include "emptyFvsPatchField.H"
+#include "emptyPointPatchField.H"
+#include "emptyFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -134,6 +136,23 @@ tmp<GeometricField<Type, fvPatchField, volMesh> > fvMeshSubset::interpolate
 }
 
 
+template<class Type>
+tmp<GeometricField<Type, fvPatchField, volMesh> > fvMeshSubset::interpolate
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    return interpolate
+    (
+        vf,
+        subMesh(),
+        patchMap(),
+        cellMap(),
+        faceMap()
+    );
+}
+
+
 template<class Type>
 tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > fvMeshSubset::interpolate
 (
@@ -259,34 +278,141 @@ tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > fvMeshSubset::interpolate
 
 
 template<class Type>
-tmp<GeometricField<Type, fvPatchField, volMesh> > fvMeshSubset::interpolate
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > fvMeshSubset::interpolate
 (
-    const GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sf
 ) const
 {
     return interpolate
     (
-        vf,
+        sf,
         subMesh(),
         patchMap(),
-        cellMap(),
         faceMap()
     );
 }
 
 
 template<class Type>
-tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > fvMeshSubset::interpolate
+tmp<GeometricField<Type, pointPatchField, pointMesh> >
+fvMeshSubset::interpolate
 (
-    const GeometricField<Type, fvsPatchField, surfaceMesh>& sf
+    const GeometricField<Type, pointPatchField, pointMesh>& vf,
+    const pointMesh& sMesh,
+    const objectRegistry& reg,
+    const labelList& patchMap,
+    const labelList& pointMap
+)
+{
+    // Create and map the internal-field values
+    Field<Type> internalField(vf.internalField(), pointMap);
+
+    // Create and map the patch field values
+    PtrList<pointPatchField<Type> > patchFields(patchMap.size());
+
+    forAll (patchFields, patchI)
+    {
+        // Set the first one by hand as it corresponds to the
+        // exposed internal faces.  Additional interpolation can be put here
+        // as necessary.  
+        if (patchMap[patchI] == -1)
+        {
+            patchFields.set
+            (
+                patchI,
+                new emptyPointPatchField<Type>
+                (
+                    sMesh.boundary()[patchI],
+                    DimensionedField<Type, pointMesh>::null()
+                )
+            );
+        }
+        else
+        {
+            // Construct addressing
+            const pointPatch& basePatch =
+                vf.mesh().boundary()[patchMap[patchI]];
+
+            const labelList& meshPoints = basePatch.meshPoints();
+
+            // Make addressing from mesh to patch point
+            Map<label> meshPointMap(2*meshPoints.size());
+            forAll(meshPoints, localI)
+            {
+                meshPointMap.insert(meshPoints[localI], localI);
+            }
+
+            // Find which subpatch points originate from which patch point
+            const pointPatch& subPatch = sMesh.boundary()[patchI];
+            const labelList& subMeshPoints = subPatch.meshPoints();
+
+            // If mapped from outside patch use point 0 for lack of better.
+            labelList directAddressing(subPatch.size(), 0);
+
+            forAll(subMeshPoints, localI)
+            {
+                // Get mesh point on original mesh.
+                label meshPointI = pointMap[subMeshPoints[localI]];
+
+                Map<label>::const_iterator iter = meshPointMap.find(meshPointI);
+
+                if (iter != meshPointMap.end())
+                {
+                    directAddressing[localI] = iter();
+                }
+            }
+
+            patchFields.set
+            (
+                patchI,
+                pointPatchField<Type>::New
+                (
+                    vf.boundaryField()[patchMap[patchI]],
+                    subPatch,
+                    DimensionedField<Type, pointMesh>::null(),
+                    pointPatchFieldSubset(directAddressing)
+                )
+            );
+        }
+    }
+
+    // Create the complete field from the pieces
+    tmp<GeometricField<Type, pointPatchField, pointMesh> > tresF
+    (
+        new GeometricField<Type, pointPatchField, pointMesh>
+        (
+            IOobject
+            (
+                "subset"+vf.name(),
+                vf.time().timeName(),
+                reg,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            sMesh,
+            vf.dimensions(),
+            internalField,
+            patchFields
+        )
+    );
+
+    return tresF;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, pointPatchField, pointMesh> > fvMeshSubset::interpolate
+(
+    const GeometricField<Type, pointPatchField, pointMesh>& sf
 ) const
 {
     return interpolate
     (
         sf,
-        subMesh(),
+        subPointMesh(),     // subsetted point mesh
+        subMesh(),          // registry (pointfields are stored on the polyMesh)
         patchMap(),
-        faceMap()
+        pointMap()
     );
 }
 
diff --git a/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
index 876b1333834..e45c4e470e3 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
@@ -77,7 +77,7 @@ surfaceSlipDisplacementPointPatchVectorField
     const dictionary& dict
 )
 :
-    pointPatchVectorField(p, iF),
+    pointPatchVectorField(p, iF, dict),
     surfaceNames_(dict.lookup("projectSurfaces")),
     projectMode_(followModeNames_.read(dict.lookup("followMode"))),
     projectDir_(dict.lookup("projectDirection")),
-- 
GitLab