diff --git a/src/functionObjects/field/fieldAverage/fieldAverage.H b/src/functionObjects/field/fieldAverage/fieldAverage.H
index 6985a21670..6d6a748d34 100644
--- a/src/functionObjects/field/fieldAverage/fieldAverage.H
+++ b/src/functionObjects/field/fieldAverage/fieldAverage.H
@@ -32,7 +32,9 @@ Group
 Description
     Calculates average quantities for a user-specified selection of volumetric
     and surface fields.
-    With the %subRegion option, also supports fields on a surfMesh.
+    With the %subRegion option, also supports fields on functionObject
+    surface output (eg, sampledSurfaces).
+
 
     Fields are entered as a list of sub-dictionaries, which indicate the type of
     averages to perform, and can be updated during the calculation.  The current
@@ -104,15 +106,15 @@ Usage
 
     Where the entries comprise:
     \table
-        Property          | Description               | Required | Default
-        type              | type name: fieldAverage              | yes |
-        restartOnRestart  | Restart the averaging on restart     | no  | no
-        restartOnOutput   | Restart the averaging on output      | no  | no
-        periodicRestart   | Periodically restart the averaging   | no  | no
-        restartPeriod     | Periodic restart period              | conditional |
-        restartTime       | One-shot reset of the averaging | no | great
-        fields            | list of fields and averaging options | yes |
-        subRegion         | name of a sub-region such as a surface name | no |
+        Property        | Description                       | Required | Default
+        type            | Type name: fieldAverage           | yes |
+        restartOnRestart| Restart the averaging on restart  | no  | no
+        restartOnOutput | Restart the averaging on output   | no  | no
+        periodicRestart | Periodically restart the averaging | no  | no
+        restartPeriod   | Periodic restart period            | conditional |
+        restartTime     | One-shot reset of the averaging    | no | great
+        fields          | list of fields and averaging options | yes |
+        subRegion       | Name for alternative objectRegistry  | no  | ""
     \endtable
 
 
diff --git a/src/functionObjects/field/fieldAverage/fieldAverageTemplates.C b/src/functionObjects/field/fieldAverage/fieldAverageTemplates.C
index 6b7fe6e1be..edfa9628e1 100644
--- a/src/functionObjects/field/fieldAverage/fieldAverageTemplates.C
+++ b/src/functionObjects/field/fieldAverage/fieldAverageTemplates.C
@@ -28,7 +28,7 @@ License
 #include "fieldAverageItem.H"
 #include "volFields.H"
 #include "surfaceFields.H"
-#include "surfFields.H"
+#include "polySurfaceFields.H"
 #include "OFstream.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -99,7 +99,7 @@ void Foam::functionObjects::fieldAverage::addMeanField
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
-    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+    typedef DimensionedField<Type, polySurfaceGeoMesh> SurfFieldType;
 
     if (item.mean())
     {
@@ -169,7 +169,7 @@ void Foam::functionObjects::fieldAverage::restoreWindowFields
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
-    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+    typedef DimensionedField<Type, polySurfaceGeoMesh> SurfFieldType;
 
     if (item.window() > 0)
     {
@@ -243,11 +243,11 @@ void Foam::functionObjects::fieldAverage::addPrime2MeanField
 {
     typedef GeometricField<Type1, fvPatchField, volMesh> VolFieldType1;
     typedef GeometricField<Type1, fvsPatchField, surfaceMesh> SurfaceFieldType1;
-    typedef DimensionedField<Type1, surfGeoMesh> SurfFieldType1;
+    typedef DimensionedField<Type1, polySurfaceGeoMesh> SurfFieldType1;
 
     typedef GeometricField<Type2, fvPatchField, volMesh> VolFieldType2;
     typedef GeometricField<Type2, fvsPatchField, surfaceMesh> SurfaceFieldType2;
-    typedef DimensionedField<Type2, surfGeoMesh> SurfFieldType2;
+    typedef DimensionedField<Type2, polySurfaceGeoMesh> SurfFieldType2;
 
     if (item.prime2Mean())
     {
@@ -312,7 +312,7 @@ void Foam::functionObjects::fieldAverage::storeWindowFields()
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
-    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+    typedef DimensionedField<Type, polySurfaceGeoMesh> SurfFieldType;
 
     for (fieldAverageItem& item : faItems_)
     {
@@ -331,7 +331,7 @@ void Foam::functionObjects::fieldAverage::calculateMeanFields() const
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
-    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+    typedef DimensionedField<Type, polySurfaceGeoMesh> SurfFieldType;
 
     for (const fieldAverageItem& item : faItems_)
     {
@@ -347,11 +347,11 @@ void Foam::functionObjects::fieldAverage::calculatePrime2MeanFields() const
 {
     typedef GeometricField<Type1, fvPatchField, volMesh> VolFieldType1;
     typedef GeometricField<Type1, fvsPatchField, surfaceMesh> SurfaceFieldType1;
-    typedef DimensionedField<Type1, surfGeoMesh> SurfFieldType1;
+    typedef DimensionedField<Type1, polySurfaceGeoMesh> SurfFieldType1;
 
     typedef GeometricField<Type2, fvPatchField, volMesh> VolFieldType2;
     typedef GeometricField<Type2, fvsPatchField, surfaceMesh> SurfaceFieldType2;
-    typedef DimensionedField<Type2, surfGeoMesh> SurfFieldType2;
+    typedef DimensionedField<Type2, polySurfaceGeoMesh> SurfFieldType2;
 
     for (const fieldAverageItem& item : faItems_)
     {
@@ -391,11 +391,11 @@ void Foam::functionObjects::fieldAverage::addMeanSqrToPrime2Mean() const
 {
     typedef GeometricField<Type1, fvPatchField, volMesh> VolFieldType1;
     typedef GeometricField<Type1, fvsPatchField, surfaceMesh> SurfaceFieldType1;
-    typedef DimensionedField<Type1, surfGeoMesh> SurfFieldType1;
+    typedef DimensionedField<Type1, polySurfaceGeoMesh> SurfFieldType1;
 
     typedef GeometricField<Type2, fvPatchField, volMesh> VolFieldType2;
     typedef GeometricField<Type2, fvsPatchField, surfaceMesh> SurfaceFieldType2;
-    typedef DimensionedField<Type2, surfGeoMesh> SurfFieldType2;
+    typedef DimensionedField<Type2, polySurfaceGeoMesh> SurfFieldType2;
 
     for (const fieldAverageItem& item : faItems_)
     {
@@ -431,7 +431,7 @@ void Foam::functionObjects::fieldAverage::writeFields() const
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
-    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+    typedef DimensionedField<Type, polySurfaceGeoMesh> SurfFieldType;
 
     for (const fieldAverageItem& item : faItems_)
     {
diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
index 056b00a0c5..65225a544a 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
@@ -60,7 +60,7 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::regionTypeNames_
 ({
     { regionTypes::stFaceZone, "faceZone" },
     { regionTypes::stPatch, "patch" },
-    { regionTypes::stSurface, "surface" },
+    { regionTypes::stObject, "functionObjectSurface" },
     { regionTypes::stSampled, "sampledSurface" },
 });
 
@@ -118,9 +118,9 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::postOperationTypeNames_
 const Foam::objectRegistry&
 Foam::functionObjects::fieldValues::surfaceFieldValue::obr() const
 {
-    if (stSurface == regionType_)
+    if (stObject == regionType_)
     {
-        return mesh_.lookupObject<objectRegistry>(regionName_);
+        return storedObjects().lookupObject<polySurface>(regionName_);
     }
 
     return mesh_;
@@ -368,9 +368,9 @@ combineSurfaceGeometry
     pointField& points
 ) const
 {
-    if (stSurface == regionType_)
+    if (stObject == regionType_)
     {
-        const surfMesh& s = dynamicCast<const surfMesh>(obr());
+        const polySurface& s = dynamicCast<const polySurface>(obr());
 
         if (Pstream::parRun())
         {
@@ -436,9 +436,9 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::totalArea() const
 {
     scalar totalArea = 0;
 
-    if (stSurface == regionType_)
+    if (stObject == regionType_)
     {
-        const surfMesh& s = dynamicCast<const surfMesh>(obr());
+        const polySurface& s = dynamicCast<const polySurface>(obr());
 
         totalArea = gSum(s.magSf());
     }
@@ -505,9 +505,9 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::update()
             setPatchFaces();
             break;
         }
-        case stSurface:
+        case stObject:
         {
-            const surfMesh& s = dynamicCast<const surfMesh>(obr());
+            const polySurface& s = dynamicCast<const polySurface>(obr());
             nFaces_ = returnReduce(s.size(), sumOp<label>());
             break;
         }
@@ -1007,9 +1007,9 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::write()
     vectorField Sf;
     if (usesSf())
     {
-        if (stSurface == regionType_)
+        if (stObject == regionType_)
         {
-            const surfMesh& s = dynamicCast<const surfMesh>(obr());
+            const polySurface& s = dynamicCast<const polySurface>(obr());
             Sf = s.Sf();
         }
         else if (sampledPtr_.valid())
diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
index 211c9ae5ca..7d64f9cc70 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
@@ -98,7 +98,7 @@ Usage
     \plaintable
         faceZone     | The \b name entry to specify the faceZone
         patch        | The \b name entry to specify the patch
-        surface      | The \b name entry to specify the surfMesh
+        functionObjectSurface | The \b name entry to specify a polySurface
         sampledSurface | A \b sampledSurfaceDict sub-dictionary and \b name
     \endplaintable
 
@@ -134,13 +134,14 @@ Note
     - The values reported by the areaNormalAverage and areaNormalIntegrate
       operations are written as the first component of a field with the same
       rank as the input field.
-    - faces on empty patches get ignored
-    - if the field is a volField the \c faceZone can only consist of boundary
+    - Faces on empty patches get ignored
+    - If the field is a volField the \c faceZone can only consist of boundary
       faces
-    - Using \c surface:
+    - Using \c functionObjectSurface:
       - The keyword %subRegion should not be used to select surfaces.
-        Specify instead the regionType 'surface' and provide the surface name.
-    - using \c sampledSurface:
+        Instead specify the regionType 'functionObjectSurface' and provide
+        the name.
+    - Using \c sampledSurface:
         - not available for surface fields
         - if interpolate=true they use \c interpolationCellPoint
           otherwise they use cell values
@@ -180,11 +181,11 @@ SourceFiles
 
 #include "fieldValue.H"
 #include "Enum.H"
-#include "meshedSurf.H"
 #include "surfaceMesh.H"
+#include "polySurface.H"
 #include "fvsPatchField.H"
 #include "volFieldsFwd.H"
-#include "surfFieldsFwd.H"
+#include "polySurfaceFieldsFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -217,7 +218,7 @@ public:
         {
             stFaceZone = 0x01,   //!< Calculate on a faceZone
             stPatch    = 0x02,   //!< Calculate on a patch
-            stSurface  = 0x11,   //!< Calculate with fields on a surfMesh
+            stObject   = 0x11,   //!< Calculate with function object surface
             stSampled  = 0x12    //!< Sample onto surface and calculate
         };
 
diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueI.H b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueI.H
index 7975fd66a2..a5fa1b4353 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueI.H
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueI.H
@@ -26,13 +26,19 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "Time.H"
+#include "sampledSurface.H"
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 inline bool Foam::functionObjects::fieldValues::surfaceFieldValue::
 withSurfaceFields() const
 {
-    return (stFaceZone == regionType_ || stPatch == regionType_);
+    return
+    (
+        stFaceZone == regionType_
+     || stPatch == regionType_
+     || (sampledPtr_.valid() && sampledPtr_->withSurfaceFields())
+    );
 }
 
 
diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C
index d665e74d58..d50bc1dee9 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C
@@ -27,7 +27,7 @@ License
 
 #include "surfaceFieldValue.H"
 #include "surfaceFields.H"
-#include "surfFields.H"
+#include "polySurfaceFields.H"
 #include "volFields.H"
 #include "sampledSurface.H"
 #include "surfaceWriter.H"
@@ -58,7 +58,7 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::validField
 {
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> sf;
     typedef GeometricField<Type, fvPatchField, volMesh> vf;
-    typedef DimensionedField<Type, surfGeoMesh> smt;
+    typedef DimensionedField<Type, polySurfaceGeoMesh> smt;
 
     return
     (
@@ -79,7 +79,7 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::getFieldValues
 {
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> sf;
     typedef GeometricField<Type, fvPatchField, volMesh> vf;
-    typedef DimensionedField<Type, surfGeoMesh> smt;
+    typedef DimensionedField<Type, polySurfaceGeoMesh> smt;
 
     if (foundObject<smt>(fieldName))
     {
diff --git a/src/functionObjects/field/fieldsExpression/fieldsExpressionTemplates.C b/src/functionObjects/field/fieldsExpression/fieldsExpressionTemplates.C
index 6704dd75f0..17a11a9b2d 100644
--- a/src/functionObjects/field/fieldsExpression/fieldsExpressionTemplates.C
+++ b/src/functionObjects/field/fieldsExpression/fieldsExpressionTemplates.C
@@ -27,7 +27,7 @@ License
 
 #include "volFields.H"
 #include "surfaceFields.H"
-#include "surfFields.H"
+#include "polySurfaceFields.H"
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
@@ -36,7 +36,7 @@ bool Foam::functionObjects::fieldsExpression::calcFieldTypes(FOType& fo)
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
-    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+    typedef DimensionedField<Type, polySurfaceGeoMesh> SurfFieldType;
 
     if (foundObject<VolFieldType>(fieldNames_[0]))
     {
diff --git a/src/functionObjects/field/fluxSummary/fluxSummary.C b/src/functionObjects/field/fluxSummary/fluxSummary.C
index 9aff420a03..b31467f38c 100644
--- a/src/functionObjects/field/fluxSummary/fluxSummary.C
+++ b/src/functionObjects/field/fluxSummary/fluxSummary.C
@@ -27,8 +27,7 @@ License
 
 #include "fluxSummary.H"
 #include "surfaceFields.H"
-#include "surfFields.H"
-#include "surfMesh.H"
+#include "polySurfaceFields.H"
 #include "dictionary.H"
 #include "Time.H"
 #include "syncTools.H"
@@ -65,6 +64,7 @@ Foam::functionObjects::fluxSummary::modeTypeNames_
     { modeType::mdFaceZone , "faceZone" },
     { modeType::mdFaceZoneAndDirection, "faceZoneAndDirection" },
     { modeType::mdCellZoneAndDirection, "cellZoneAndDirection" },
+    { modeType::mdSurface, "functionObjectSurface" },
     { modeType::mdSurface, "surface" },
     { modeType::mdSurfaceAndDirection, "surfaceAndDirection" },
 });
@@ -116,14 +116,15 @@ void Foam::functionObjects::fluxSummary::initialiseSurface
     DynamicList<boolList>& faceFlip
 ) const
 {
-    const surfMesh* surfptr = mesh_.cfindObject<surfMesh>(surfName);
+    const polySurface* surfptr =
+        storedObjects().cfindObject<polySurface>(surfName);
 
     if (!surfptr)
     {
         FatalErrorInFunction
             << "Unable to find surface " << surfName
             << ".  Valid surfaces: "
-            << mesh_.sortedNames<surfMesh>() << nl
+            << storedObjects().sortedNames<polySurface>() << nl
             << exit(FatalError);
     }
 
@@ -142,14 +143,15 @@ void Foam::functionObjects::fluxSummary::initialiseSurfaceAndDirection
     DynamicList<boolList>& faceFlip
 ) const
 {
-    const surfMesh* surfptr = mesh_.cfindObject<surfMesh>(surfName);
+    const polySurface* surfptr =
+        storedObjects().cfindObject<polySurface>(surfName);
 
     if (!surfptr)
     {
         FatalErrorInFunction
             << "Unable to find surface " << surfName
             << ".  Valid surfaces: "
-            << mesh_.sortedNames<surfMesh>() << nl
+            << storedObjects().sortedNames<polySurface>() << nl
             << exit(FatalError);
     }
 
@@ -661,7 +663,9 @@ Foam::scalar Foam::functionObjects::fluxSummary::totalArea
 
     if (isSurfaceMode())
     {
-        const surfMesh& s = mesh_.lookupObject<surfMesh>(zoneNames_[idx]);
+        const polySurface& s =
+            storedObjects().lookupObject<polySurface>(zoneNames_[idx]);
+
         sumMagSf = sum(s.magSf());
     }
     else
@@ -693,11 +697,12 @@ Foam::scalar Foam::functionObjects::fluxSummary::totalArea
 
 bool Foam::functionObjects::fluxSummary::surfaceModeWrite()
 {
-    if (zoneNames_.size())
+    for (const word& surfName : zoneNames_)
     {
-        const label surfi = 0;
-        const surfMesh& s = mesh_.lookupObject<surfMesh>(zoneNames_[surfi]);
-        const surfVectorField& phi = s.lookupObject<surfVectorField>(phiName_);
+        const polySurface& s =
+            storedObjects().lookupObject<polySurface>(surfName);
+
+        const auto& phi = s.lookupObject<polySurfaceVectorField>(phiName_);
 
         Log << type() << ' ' << name() << ' '
             << checkFlowType(phi.dimensions(), phi.name()) << " write:" << nl;
@@ -706,8 +711,10 @@ bool Foam::functionObjects::fluxSummary::surfaceModeWrite()
 
     forAll(zoneNames_, surfi)
     {
-        const surfMesh& s = mesh_.lookupObject<surfMesh>(zoneNames_[surfi]);
-        const surfVectorField& phi = s.lookupObject<surfVectorField>(phiName_);
+        const polySurface& s =
+            storedObjects().lookupObject<polySurface>(zoneNames_[surfi]);
+
+        const auto& phi = s.lookupObject<polySurfaceVectorField>(phiName_);
 
         checkFlowType(phi.dimensions(), phi.name());
 
diff --git a/src/functionObjects/field/fluxSummary/fluxSummary.H b/src/functionObjects/field/fluxSummary/fluxSummary.H
index 0e76c0c49a..f1e7f81deb 100644
--- a/src/functionObjects/field/fluxSummary/fluxSummary.H
+++ b/src/functionObjects/field/fluxSummary/fluxSummary.H
@@ -120,14 +120,14 @@ public:
         //- Face mode type
         enum modeType
         {
-            mdFaceZone,                 //!< face zone
-            mdFaceZoneAndDirection,     //!< face zone with prescribed direction
-            mdCellZoneAndDirection,     //!< cell zone with prescribed direction
-            mdSurface,                  //!< surfMesh
-            mdSurfaceAndDirection       //!< surfMesh with prescribed direction
+            mdFaceZone,             //!< Face zone
+            mdFaceZoneAndDirection, //!< Face zone with prescribed direction
+            mdCellZoneAndDirection, //!< Cell zone with prescribed direction
+            mdSurface,              //!< A functionObject surface
+            mdSurfaceAndDirection   //!< A surface with prescribed direction
         };
 
-        //- Mode type names
+        //- Face mode names
         static const Enum<modeType> modeTypeNames_;
 
 
diff --git a/src/functionObjects/field/mag/mag.H b/src/functionObjects/field/mag/mag.H
index e8b09cf855..cff6e244ae 100644
--- a/src/functionObjects/field/mag/mag.H
+++ b/src/functionObjects/field/mag/mag.H
@@ -34,7 +34,8 @@ Description
 
     The operation can be applied to any volume or surface fields generating a
     volume or surface scalar field.
-    With the %subRegion option, also supports fields on a surfMesh.
+    With the %subRegion option, also supports fields on functionObject
+    surface output (eg, sampledSurfaces).
 
 See also
     Foam::functionObjects::fvMeshFunctionObject
diff --git a/src/functionObjects/field/mag/magTemplates.C b/src/functionObjects/field/mag/magTemplates.C
index 82f59f601e..3a7c840668 100644
--- a/src/functionObjects/field/mag/magTemplates.C
+++ b/src/functionObjects/field/mag/magTemplates.C
@@ -27,7 +27,7 @@ License
 
 #include "volFields.H"
 #include "surfaceFields.H"
-#include "surfFields.H"
+#include "polySurfaceFields.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -36,7 +36,7 @@ bool Foam::functionObjects::mag::calcMag()
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
-    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+    typedef DimensionedField<Type, polySurfaceGeoMesh> SurfFieldType;
 
     if (foundObject<VolFieldType>(fieldName_, false))
     {
diff --git a/src/functionObjects/field/magSqr/magSqr.H b/src/functionObjects/field/magSqr/magSqr.H
index 5e92bbbd26..146c1b2c44 100644
--- a/src/functionObjects/field/magSqr/magSqr.H
+++ b/src/functionObjects/field/magSqr/magSqr.H
@@ -34,7 +34,8 @@ Description
 
     The operation can be applied to any volume or surface field generating a
     volume or surface scalar field.
-    With the %subRegion option, also supports fields on a surfMesh.
+    With the %subRegion option, also supports fields on functionObject
+    surface output (eg, sampledSurfaces).
 
 See also
     Foam::functionObjects::fvMeshFunctionObject
diff --git a/src/functionObjects/field/magSqr/magSqrTemplates.C b/src/functionObjects/field/magSqr/magSqrTemplates.C
index 1f636b2496..2a38479282 100644
--- a/src/functionObjects/field/magSqr/magSqrTemplates.C
+++ b/src/functionObjects/field/magSqr/magSqrTemplates.C
@@ -27,7 +27,7 @@ License
 
 #include "volFields.H"
 #include "surfaceFields.H"
-#include "surfFields.H"
+#include "polySurfaceFields.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -36,7 +36,7 @@ bool Foam::functionObjects::magSqr::calcMagSqr()
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
-    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+    typedef DimensionedField<Type, polySurfaceGeoMesh> SurfFieldType;
 
     if (foundObject<VolFieldType>(fieldName_, false))
     {
diff --git a/src/functionObjects/field/readFields/readFieldsTemplates.C b/src/functionObjects/field/readFields/readFieldsTemplates.C
index 81e2963089..662e5ad6ff 100644
--- a/src/functionObjects/field/readFields/readFieldsTemplates.C
+++ b/src/functionObjects/field/readFields/readFieldsTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2015-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -28,7 +28,6 @@ License
 #include "readFields.H"
 #include "volFields.H"
 #include "surfaceFields.H"
-#include "surfFields.H"
 #include "Time.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -38,7 +37,7 @@ bool Foam::functionObjects::readFields::loadField(const word& fieldName)
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
-    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+    /// typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
 
     if (foundObject<VolFieldType>(fieldName))
     {
@@ -54,13 +53,13 @@ bool Foam::functionObjects::readFields::loadField(const word& fieldName)
             << " " << fieldName << " already exists in database"
             << " already in database" << endl;
     }
-    else if (foundObject<SurfFieldType>(fieldName))
-    {
-        DebugInfo
-            << "readFields: " << SurfFieldType::typeName
-            << " " << fieldName << " already exists in database"
-            << " already in database" << endl;
-    }
+    /// else if (foundObject<SurfFieldType>(fieldName))
+    /// {
+    ///     DebugInfo
+    ///         << "readFields: " << SurfFieldType::typeName
+    ///         << " " << fieldName << " already exists in database"
+    ///         << " already in database" << endl;
+    /// }
     else
     {
         IOobject fieldHeader
@@ -76,35 +75,32 @@ bool Foam::functionObjects::readFields::loadField(const word& fieldName)
         {
             // Store field on mesh database
             Log << "    Reading " << fieldName << endl;
-            VolFieldType* vfPtr(new VolFieldType(fieldHeader, mesh_));
-            mesh_.objectRegistry::store(vfPtr);
+            VolFieldType* fldPtr(new VolFieldType(fieldHeader, mesh_));
+            mesh_.objectRegistry::store(fldPtr);
             return true;
         }
         else if (fieldHeader.typeHeaderOk<SurfaceFieldType>(true, true, false))
         {
             // Store field on mesh database
             Log << "    Reading " << fieldName << endl;
-            SurfaceFieldType* sfPtr(new SurfaceFieldType(fieldHeader, mesh_));
-            mesh_.objectRegistry::store(sfPtr);
+            SurfaceFieldType* fldPtr(new SurfaceFieldType(fieldHeader, mesh_));
+            mesh_.objectRegistry::store(fldPtr);
             return true;
         }
-        else if (fieldHeader.typeHeaderOk<SurfFieldType>(true, true, false))
-        {
-            if (isA<surfMesh>(obr()))
-            {
-                const surfMesh& s = dynamicCast<const surfMesh>(obr());
-
-                // Store field on surfMesh database
-                Log << "    Reading " << fieldName << endl;
-                SurfFieldType* sfPtr(new SurfFieldType(fieldHeader, s));
-                s.store(sfPtr);
-                return true;
-            }
-            else
-            {
-                return false;
-            }
-        }
+        /// else if (fieldHeader.typeHeaderOk<SurfFieldType>(true, true, false))
+        /// {
+        ///     const surfMesh* surfptr = isA<surfMesh>(obr());
+        ///     if (surfptr)
+        ///     {
+        ///         const surfMesh& s = surfptr;
+        ///
+        ///         // Store field on surfMesh database
+        ///         Log << "    Reading " << fieldName << endl;
+        ///         SurfFieldType* fldPtr(new SurfFieldType(fieldHeader, s));
+        ///         s.store(fldPtr);
+        ///         return true;
+        ///     }
+        /// }
     }
 
     return false;
diff --git a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/avg-tracer0 b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/avg-tracer0
index b0063ac36e..860916f58a 100644
--- a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/avg-tracer0
+++ b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/avg-tracer0
@@ -3,20 +3,20 @@
 // Averaging of volume fields or fields on an externally loaded surface.
 avg-tracer0-0.25
 {
-    ${settings_avg}
-    subRegion   plane-0.25;
+    ${__settings_avg}
+    subRegion   sampled.plane-0.25;
 }
 
 avg-tracer0-0.45
 {
-    ${settings_avg}
-    subRegion   plane-0.45;
+    ${__settings_avg}
+    subRegion   sampled.plane-0.45;
 }
 
 avg-tracer0-0.55
 {
-    ${settings_avg}
-    subRegion   plane-0.55;
+    ${__settings_avg}
+    subRegion   sampled.plane-0.55;
 }
 
 
diff --git a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/controlDict b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/controlDict
index f5918a37aa..46851bcc18 100644
--- a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/controlDict
+++ b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/controlDict
@@ -55,8 +55,6 @@ endTime         1200;
 deltaT          1;
 writeInterval   100;
 
-#include "sampleControls"
-
 functions
 {
     #include "abort"
diff --git a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fieldTransfer b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fieldTransfer
index b1eae09bb3..acbca5e483 100644
--- a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fieldTransfer
+++ b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fieldTransfer
@@ -1,17 +1,19 @@
 // -*- C++ -*-
 
-// ************************************************************************* //
-
-// Transcribe volume fields to surfaces.
-fieldTransfer
+// Sample volume fields to surfaces and hold on registry.
+sampled
 {
-    type    surfMeshes;
+    type    surfaces;
     libs    ("libsampling.so");
     log     true;
-    writeControl    none;
-    createOnRead    true;
+    verbose true;
+
     executeControl  timeStep;
     executeInterval 1;
+    writeControl    none;
+    surfaceFormat   none;
+
+    sampleOnExecute true;
 
     // Includes a rhoU derived field
     fields  (rho U tracer0 rhoU);
@@ -19,31 +21,31 @@ fieldTransfer
     // Remove derived fields we created prior
     removeFieldsOnExecute   (pTotal);
 
-
-    baseCfg
+    _sampleMesh
     {
-        type        sampledTriSurfaceMesh;
-        source      cells;
+        type    sampledTriSurfaceMesh;
+        source  cells;
+        store   true;
     }
 
     surfaces
-    (
+    {
         plane-0.25
         {
-            $baseCfg
+            ${_sampleMesh}
             surface     plane-0.25.stl;
         }
         plane-0.45
         {
-            $baseCfg
+            ${_sampleMesh}
             surface     plane-0.45.stl;
         }
         plane-0.55
         {
-            $baseCfg
+            ${_sampleMesh}
             surface     plane-0.55.stl;
         }
-    );
+    }
 }
 
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fluxSummary b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fluxSummary
new file mode 100644
index 0000000000..7bd6592ede
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fluxSummary
@@ -0,0 +1,29 @@
+// -*- C++ -*-
+//
+// Summary of fluxes
+
+flux
+{
+    type            fluxSummary;
+    libs            ("libfieldFunctionObjects.so");
+    log             on;
+    enabled         true;
+
+    executeControl  none;
+    writeControl    timeStep;
+    writeInterval   4;
+
+    mode            functionObjectSurface;
+
+    phi             rhoU;
+
+    surfaces
+    (
+        sampled.plane-0.25
+        sampled.plane-0.45
+        sampled.plane-0.55
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sampling b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sampling
index 233c5e68ca..e2547f7cb8 100644
--- a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sampling
+++ b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sampling
@@ -4,12 +4,20 @@
 
 #include "derivedFields"
 #include "fieldTransfer"
+#include "fluxSummary"
+
+#include "surfaceFieldValue.cfg"
+
 #include "avg-tracer0"
 #include "sum-tracer0"
 // #include "sum-tracer0Mean"
 
+
 // #include "averaging.1612"            // generate time-averaged fields
 // #include "spatialDeviation.1612"     // spatial deviation compared to a mean field
 // #include "statistics.1612"           // transcribe volume fields -> surface
 
+// Cleanup
+#remove "__.*"
+
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sum-tracer0 b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sum-tracer0
index 0e19fb6392..82d306231c 100644
--- a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sum-tracer0
+++ b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sum-tracer0
@@ -4,10 +4,10 @@
 
 sum-tracer0-0.25
 {
-    $settings_stats
+    ${__settings_stats}
 
-    regionType      surface;
-    name            plane-0.25;
+    regionType      functionObjectSurface;
+    name            sampled.plane-0.25;
 
     operation       weightedAreaIntegrate;
     // postOperation   sqrt;
@@ -19,10 +19,10 @@ sum-tracer0-0.25
 
 sum-tracer0-0.45
 {
-    $settings_stats
+    ${__settings_stats}
 
-    regionType      surface;
-    name            plane-0.45;
+    regionType      functionObjectSurface;
+    name            sampled.plane-0.45;
 
     operation       weightedAreaIntegrate;
     // postOperation   sqrt;
@@ -33,10 +33,10 @@ sum-tracer0-0.45
 
 sum-tracer0-0.55
 {
-    $settings_stats
+    ${__settings_stats}
 
-    regionType      surface;
-    name            plane-0.55;
+    regionType      functionObjectSurface;
+    name            sampled.plane-0.55;
 
     operation       weightedAreaIntegrate;
     // postOperation   sqrt;
diff --git a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sampleControls b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/surfaceFieldValue.cfg
similarity index 98%
rename from tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sampleControls
rename to tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/surfaceFieldValue.cfg
index f1fba81312..9c8c0ed2f0 100644
--- a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sampleControls
+++ b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/surfaceFieldValue.cfg
@@ -16,7 +16,7 @@
 //   restartTime:
 //     - a 'one-shot' reset at a particular time
 //
-settings_avg
+__settings_avg
 {
     type            fieldAverage;
     libs            ("libfieldFunctionObjects.so");
@@ -77,7 +77,7 @@ settings_avg
 //     The sqrt operation is useful when determining RMS values.
 //
 // The 'output/write' control triggers the calculation.
-settings_stats
+__settings_stats
 {
     type            surfaceFieldValue;
     libs            ("libfieldFunctionObjects.so");
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/controlDict b/tutorials/compressible/rhoSimpleFoam/squareBend/system/controlDict
index d5dff1e013..441cddd999 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/controlDict
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/controlDict
@@ -47,10 +47,9 @@ graphFormat     raw;
 
 runTimeModifiable true;
 
-#include "sampleControls"
-
 functions
 {
+    #include "derivedFields"
     #include "sampling"
     // #include "samplingDebug"
     // #include "sampleCellCentres"
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fieldTransfer b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fieldTransfer
deleted file mode 100644
index 3733f5c8c8..0000000000
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fieldTransfer
+++ /dev/null
@@ -1,65 +0,0 @@
-// -*- C++ -*-
-
-// ************************************************************************* //
-
-// Transcribe volume fields to surfaces.
-fieldTransfer
-{
-    type    surfMeshes;
-    libs    ("libsampling.so");
-    log     true;
-    writeControl    none;
-    createOnRead    true;
-    executeControl  timeStep;
-    executeInterval 1;
-
-    fields      (p rho U T rhoU pTotal);
-
-    // Remove derived fields we created prior
-    removeFieldsOnExecute   (rhoU pTotal);
-
-    _plane
-    {
-        type    plane;
-        source  cells;
-
-        planeType   pointAndNormal;
-
-        pointAndNormalDict
-        {
-            normal (-1 0 0);
-            point  (-0.04 0 0);
-        }
-    }
-
-    surfaces
-    (
-        // Top channel
-        plane1
-        {
-            ${_plane}
-            bounds (-1 0 -1) (0 1 1);
-        }
-
-        // Bottom channel
-        plane2
-        {
-            ${_plane}
-            bounds (-1 -1 -1) (0 0 1);
-        }
-
-        // Angled plane - for general testing
-        plane3
-        {
-            type        distanceSurface;
-            distance    0;
-            signed      true;
-
-            surfaceType triSurfaceMesh;
-            surfaceName angledPlane.obj;
-        }
-    );
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampleControls b/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampleControls
deleted file mode 100644
index 23e1cf4246..0000000000
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampleControls
+++ /dev/null
@@ -1,40 +0,0 @@
-// -*- C++ -*-
-
-//   restartTime:
-//     - a 'one-shot' reset at a particular time
-//
-// fields [required]
-//     Pairs of fields to use for calculating the deviation.
-//     The fields must already exist on the surfaces.
-//
-// weightField [optional]
-//     A scalar or vector field for weighting.
-//
-// postOperation [optional]
-//     Modify the results by particular operations.
-//     (none | sqrt)
-//     The sqrt operation is useful when determining RMS values.
-//
-// The 'output/write' control triggers the calculation.
-__surfaceFieldValue
-{
-    type            surfaceFieldValue;
-    libs            ("libfieldFunctionObjects.so");
-    log             on;
-    enabled         true;
-
-    writeControl    timeStep;
-    writeInterval   1;
-
-    writeFields     false;
-    surfaceFormat   vtk;
-    // writeArea       true;
-
-    // resetOnStartUp  true;
-    // resetOnOutput   false;
-    // periodicRestart true;
-    // restartPeriod   0.0005;
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling b/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
index eeab502a7d..3f5bff95f6 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
@@ -1,16 +1,134 @@
-// -*- C++ -*-
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1812                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      sampling;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// ************************************************************************* //
+// Setup for surfaceFieldValue
+//
+//   restartTime:
+//     - a 'one-shot' reset at a particular time
+//
+// fields [required]
+//     Pairs of fields to use for calculating the deviation.
+//     The fields must already exist on the surfaces.
+//
+// weightField [optional]
+//     A scalar or vector field for weighting.
+//
+// postOperation [optional]
+//     Modify the results by particular operations.
+//     (none | sqrt)
+//     The sqrt operation is useful when determining RMS values.
+//
+// The 'output/write' control triggers the calculation.
+__surfaceFieldValue
+{
+    type            surfaceFieldValue;
+    libs            ("libfieldFunctionObjects.so");
+    log             on;
+    enabled         true;
 
-#include "derivedFields"
-#include "fieldTransfer"
+    writeControl    timeStep;
+    writeInterval   1;
+
+    writeFields     false;
+    surfaceFormat   vtk;
+    // writeArea       true;
+
+    // resetOnStartUp  true;
+    // resetOnOutput   false;
+    // periodicRestart true;
+    // restartPeriod   0.0005;
+}
+
+
+// * * * * * * * * * * * * * * * * Sampling  * * * * * * * * * * * * * * * * //
+
+// Sample volume fields to surfaces and hold on registry.
+sampled
+{
+    type    surfaces;
+    libs    ("libsampling.so");
+    log     true;
+
+    executeControl  timeStep;
+    executeInterval 1;
+    writeControl    none;
+    sampleOnExecute true;
+
+    surfaceFormat   none;
+
+    fields      (p rho U T rhoU pTotal);
+
+    // Remove derived fields we created prior
+    removeFieldsOnExecute   (rhoU pTotal);
+
+    _plane
+    {
+        type    plane;
+        source  cells;
+        store   true;
+
+        planeType   pointAndNormal;
+
+        pointAndNormalDict
+        {
+            normal (-1 0 0);
+            point  (-0.04 0 0);
+        }
+    }
+
+    surfaces
+    {
+        // Top channel
+        plane1
+        {
+            ${_plane}
+            bounds (-1 0 -1) (0 1 1);
+        }
+
+        // Bottom channel
+        plane2
+        {
+            ${_plane}
+            bounds (-1 -1 -1) (0 0 1);
+        }
+
+        // Angled plane - for general testing
+        plane3
+        {
+            type        distanceSurface;
+            distance    0;
+            signed      true;
+            store       true;
+
+            surfaceType triSurfaceMesh;
+            surfaceName angledPlane.obj;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Calculations  * * * * * * * * * * * * * * * //
 
 massflow
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane1;
+    regionType      functionObjectSurface;
+    name            sampled.plane1;
 
     operation       areaNormalIntegrate;
 
@@ -21,8 +139,8 @@ areaAverage
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane1;
+    regionType      functionObjectSurface;
+    name            sampled.plane1;
 
     operation       weightedAreaAverage;
     weightField     rhoU;
@@ -33,8 +151,8 @@ areaIntegrate
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane1;
+    regionType      functionObjectSurface;
+    name            sampled.plane1;
 
     operation       weightedAreaIntegrate;
     weightField     rhoU;
@@ -46,8 +164,8 @@ UI1
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane1;
+    regionType      functionObjectSurface;
+    name            sampled.plane1;
 
     operation       uniformity;
     fields          ( U T );
@@ -59,8 +177,8 @@ UI2
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane2;
+    regionType      functionObjectSurface;
+    name            sampled.plane2;
 
     operation       uniformity;
     fields          ( U T );
@@ -72,8 +190,8 @@ UI3
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane3;
+    regionType      functionObjectSurface;
+    name            sampled.plane3;
 
     operation       uniformity;
     fields          ( U T );
@@ -86,8 +204,8 @@ T_UI1
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane1;
+    regionType      functionObjectSurface;
+    name            sampled.plane1;
 
     operation       weightedUniformity;
     weightField     T;
@@ -100,8 +218,8 @@ rhoU_UI1
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane1;
+    regionType      functionObjectSurface;
+    name            sampled.plane1;
 
     operation       weightedUniformity;
     weightField     rhoU;
@@ -114,8 +232,8 @@ rhoU_UI2
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane2;
+    regionType      functionObjectSurface;
+    name            sampled.plane2;
 
     operation       weightedUniformity;
     weightField     rhoU;
@@ -123,4 +241,9 @@ rhoU_UI2
 }
 
 
+// * * * * * * * * * * * * * * * * * Cleanup * * * * * * * * * * * * * * * * //
+
+#remove "__.*"
+
+
 // ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/controlDict b/tutorials/incompressible/simpleFoam/squareBend/system/controlDict
index a4f3061f0d..4227f18a54 100644
--- a/tutorials/incompressible/simpleFoam/squareBend/system/controlDict
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/controlDict
@@ -47,10 +47,9 @@ graphFormat     raw;
 
 runTimeModifiable true;
 
-#include "sampleControls"
-
 functions
 {
+    #include "derivedFields"
     #include "sampling"
     // #include "samplingDebug"
 }
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/fieldTransfer b/tutorials/incompressible/simpleFoam/squareBend/system/fieldTransfer
deleted file mode 100644
index db33a6546f..0000000000
--- a/tutorials/incompressible/simpleFoam/squareBend/system/fieldTransfer
+++ /dev/null
@@ -1,65 +0,0 @@
-// -*- C++ -*-
-
-// ************************************************************************* //
-
-// Transcribe volume fields to surfaces.
-fieldTransfer
-{
-    type    surfMeshes;
-    libs    ("libsampling.so");
-    log     true;
-    writeControl    none;
-    createOnRead    true;
-    executeControl  timeStep;
-    executeInterval 1;
-
-    fields      (p U rhoU pTotal);
-
-    // Remove derived fields we created prior
-    removeFieldsOnExecute   (rhoU pTotal);
-
-    _plane
-    {
-        type    plane;
-        source  cells;
-
-        planeType   pointAndNormal;
-
-        pointAndNormalDict
-        {
-            normal (-1 0 0);
-            point  (-0.04 0 0);
-        }
-    }
-
-    surfaces
-    (
-        // Top channel
-        plane1
-        {
-            ${_plane}
-            bounds (-1 0 -1) (0 1 1);
-        }
-
-        // Bottom channel
-        plane2
-        {
-            ${_plane}
-            bounds (-1 -1 -1) (0 0 1);
-        }
-
-        // Angled plane - for general testing
-        plane3
-        {
-            type        distanceSurface;
-            distance    0;
-            signed      true;
-
-            surfaceType triSurfaceMesh;
-            surfaceName angledPlane.obj;
-        }
-    );
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/sampleControls b/tutorials/incompressible/simpleFoam/squareBend/system/sampleControls
deleted file mode 100644
index 23e1cf4246..0000000000
--- a/tutorials/incompressible/simpleFoam/squareBend/system/sampleControls
+++ /dev/null
@@ -1,40 +0,0 @@
-// -*- C++ -*-
-
-//   restartTime:
-//     - a 'one-shot' reset at a particular time
-//
-// fields [required]
-//     Pairs of fields to use for calculating the deviation.
-//     The fields must already exist on the surfaces.
-//
-// weightField [optional]
-//     A scalar or vector field for weighting.
-//
-// postOperation [optional]
-//     Modify the results by particular operations.
-//     (none | sqrt)
-//     The sqrt operation is useful when determining RMS values.
-//
-// The 'output/write' control triggers the calculation.
-__surfaceFieldValue
-{
-    type            surfaceFieldValue;
-    libs            ("libfieldFunctionObjects.so");
-    log             on;
-    enabled         true;
-
-    writeControl    timeStep;
-    writeInterval   1;
-
-    writeFields     false;
-    surfaceFormat   vtk;
-    // writeArea       true;
-
-    // resetOnStartUp  true;
-    // resetOnOutput   false;
-    // periodicRestart true;
-    // restartPeriod   0.0005;
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/sampling b/tutorials/incompressible/simpleFoam/squareBend/system/sampling
index eac6e6ce23..2e9c094497 100644
--- a/tutorials/incompressible/simpleFoam/squareBend/system/sampling
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/sampling
@@ -1,16 +1,134 @@
-// -*- C++ -*-
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1812                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      sampling;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// ************************************************************************* //
+// Setup for surfaceFieldValue
+//
+//   restartTime:
+//     - a 'one-shot' reset at a particular time
+//
+// fields [required]
+//     Pairs of fields to use for calculating the deviation.
+//     The fields must already exist on the surfaces.
+//
+// weightField [optional]
+//     A scalar or vector field for weighting.
+//
+// postOperation [optional]
+//     Modify the results by particular operations.
+//     (none | sqrt)
+//     The sqrt operation is useful when determining RMS values.
+//
+// The 'output/write' control triggers the calculation.
+__surfaceFieldValue
+{
+    type            surfaceFieldValue;
+    libs            ("libfieldFunctionObjects.so");
+    log             on;
+    enabled         true;
 
-#include "derivedFields"
-#include "fieldTransfer"
+    writeControl    timeStep;
+    writeInterval   1;
+
+    writeFields     false;
+    surfaceFormat   vtk;
+    // writeArea       true;
+
+    // resetOnStartUp  true;
+    // resetOnOutput   false;
+    // periodicRestart true;
+    // restartPeriod   0.0005;
+}
+
+
+// * * * * * * * * * * * * * * * * Sampling  * * * * * * * * * * * * * * * * //
+
+// Sample volume fields to surfaces and hold on registry.
+sampled
+{
+    type    surfaces;
+    libs    ("libsampling.so");
+    log     true;
+
+    executeControl  timeStep;
+    executeInterval 1;
+    writeControl    none;
+    sampleOnExecute true;
+
+    surfaceFormat   none;
+
+    fields      (p U rhoU pTotal);
+
+    // Remove derived fields we created prior
+    removeFieldsOnExecute   (rhoU pTotal);
+
+    _plane
+    {
+        type    plane;
+        source  cells;
+        store   true;
+
+        planeType   pointAndNormal;
+
+        pointAndNormalDict
+        {
+            normal (-1 0 0);
+            point  (-0.04 0 0);
+        }
+    }
+
+    surfaces
+    (
+        // Top channel
+        plane1
+        {
+            ${_plane}
+            bounds (-1 0 -1) (0 1 1);
+        }
+
+        // Bottom channel
+        plane2
+        {
+            ${_plane}
+            bounds (-1 -1 -1) (0 0 1);
+        }
+
+        // Angled plane - for general testing
+        plane3
+        {
+            type        distanceSurface;
+            distance    0;
+            signed      true;
+            store       true;
+
+            surfaceType triSurfaceMesh;
+            surfaceName angledPlane.obj;
+        }
+    );
+}
+
+
+// * * * * * * * * * * * * * * * Calculations  * * * * * * * * * * * * * * * //
 
 massflow
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane1;
+    regionType      functionObjectSurface;
+    name            sampled.plane1;
 
     operation       areaNormalIntegrate;
 
@@ -21,8 +139,8 @@ areaAverage
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane1;
+    regionType      functionObjectSurface;
+    name            sampled.plane1;
 
     operation       weightedAreaAverage;
     weightField     rhoU;
@@ -33,8 +151,8 @@ areaIntegrate
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane1;
+    regionType      functionObjectSurface;
+    name            sampled.plane1;
 
     operation       weightedAreaIntegrate;
     weightField     rhoU;
@@ -46,8 +164,8 @@ UI1
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane1;
+    regionType      functionObjectSurface;
+    name            sampled.plane1;
 
     operation       uniformity;
     fields          ( U p );
@@ -59,8 +177,8 @@ UI2
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane2;
+    regionType      functionObjectSurface;
+    name            sampled.plane2;
 
     operation       uniformity;
     fields          ( U p );
@@ -72,8 +190,8 @@ UI3
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane3;
+    regionType      functionObjectSurface;
+    name            sampled.plane3;
 
     operation       uniformity;
     fields          ( U p );
@@ -85,8 +203,8 @@ rhoU_UI1
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane1;
+    regionType      functionObjectSurface;
+    name            sampled.plane1;
 
     operation       weightedUniformity;
     weightField     rhoU;
@@ -99,8 +217,8 @@ rhoU_UI2
 {
     ${__surfaceFieldValue}
 
-    regionType      surface;
-    name            plane2;
+    regionType      functionObjectSurface;
+    name            sampled.plane2;
 
     operation       weightedUniformity;
     weightField     rhoU;
@@ -108,4 +226,9 @@ rhoU_UI2
 }
 
 
+// * * * * * * * * * * * * * * * * * Cleanup * * * * * * * * * * * * * * * * //
+
+#remove "__.*"
+
+
 // ************************************************************************* //