diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
index 16432d71653ff9b2d905f82f7da951c3f9ee356f..9493bf39cf4fbd96874e32d0d0f278333e144bec 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
@@ -49,7 +49,6 @@ Description
 
         // Scheme to obtain node values
         // (only used if interpolate=true for the surfaces below)
-
         interpolationScheme cell;
 
         // Output surface format
@@ -74,6 +73,22 @@ Description
     }
     \endverbatim
 
+    Entries:
+    \table
+        Property | Description                             | Required | Default
+        type     | surfaces                                | yes      |
+        surfaces | the list of sample surfaces             | recommended |
+        fields   | word/regex list of fields to sampled    | yes      |
+        sampleScheme | scheme to obtain face centre value  | no       | cell
+        interpolationScheme | scheme to obtain node values | yes      |
+        surfaceFormat | output surface format              | yes      |
+        formatOptions | dictionary of format options       | no       |
+    \endtable
+
+Note
+    The interpolationScheme is only used if interpolate=true is used by any
+    of the surfaces.
+
 See also
     Foam::surfMeshSamplers
 
diff --git a/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.C b/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.C
index 44ca76df3c7b9804307d78e6168e1c3e38674c52..a649d4f03d802ddb0beabffe2e6df20b3f576f9c 100644
--- a/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.C
+++ b/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.C
@@ -95,6 +95,93 @@ void Foam::surfMeshSamplers::checkOutNames
 //     return tmpRegistry().subRegistry(subName, true, false);
 // }
 
+bool Foam::surfMeshSamplers::add_rhoU(const word& derivedName)
+{
+    const objectRegistry& db = mesh_.thisDb();
+
+    const bool existed = db.foundObject<volVectorField>(derivedName);
+
+    if (existed)
+    {
+        return false; // Volume field already existed - not added.
+    }
+
+
+    // rhoU = rho * U
+
+    const auto* rhoPtr = mesh_.lookupObjectPtr<volScalarField>("rho");
+    const volVectorField& U = mesh_.lookupObject<volVectorField>("U");
+
+    tmp<volVectorField> tresult;
+
+    if (rhoPtr)
+    {
+        const auto& rho = *rhoPtr;
+
+        tresult = tmp<volVectorField>::New
+        (
+            derivedName, (rho * U)
+        );
+    }
+    else
+    {
+        const dimensionedScalar rho("rho", dimDensity, rhoRef_);
+
+        tresult = tmp<volVectorField>::New
+        (
+            derivedName, (rho * U)
+        );
+    }
+
+    db.store(tresult.ptr());
+
+    return !existed;
+}
+
+
+bool Foam::surfMeshSamplers::add_pTotal(const word& derivedName)
+{
+    const objectRegistry& db = mesh_.thisDb();
+
+    const bool existed = db.foundObject<volVectorField>(derivedName);
+
+    if (existed)
+    {
+        return false; // Volume field already existed - not added.
+    }
+
+    // pTotal = p + rho * U^2 / 2
+
+    const auto* rhoPtr = mesh_.lookupObjectPtr<volScalarField>("rho");
+    const volScalarField& p = mesh_.lookupObject<volScalarField>("p");
+    const volVectorField& U = mesh_.lookupObject<volVectorField>("U");
+
+    tmp<volScalarField> tresult;
+
+    if (rhoPtr)
+    {
+        const auto& rho = *rhoPtr;
+
+        tresult = tmp<volScalarField>::New
+        (
+            derivedName, (p + 0.5 * rho * magSqr(U))
+        );
+    }
+    else
+    {
+        const dimensionedScalar rho("rho", dimDensity, rhoRef_);
+
+        tresult = tmp<volScalarField>::New
+        (
+            derivedName, (rho * (p + 0.5 * magSqr(U)))
+        );
+    }
+
+    db.store(tresult.ptr());
+
+    return !existed;
+}
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -110,7 +197,8 @@ Foam::surfMeshSamplers::surfMeshSamplers
     mesh_(refCast<const fvMesh>(obr_)),
     fieldSelection_(),
     derivedNames_(),
-    sampleScheme_(word::null)
+    sampleScheme_(word::null),
+    rhoRef_(1.0)
 {
     read(dict);
 }
@@ -128,7 +216,8 @@ Foam::surfMeshSamplers::surfMeshSamplers
     mesh_(refCast<const fvMesh>(obr)),
     fieldSelection_(),
     derivedNames_(),
-    sampleScheme_(word::null)
+    sampleScheme_(word::null),
+    rhoRef_(1.0)
 {
     read(dict);
 }
@@ -157,79 +246,23 @@ bool Foam::surfMeshSamplers::execute()
 
     for (const word& derivedName : derivedNames_)
     {
+        // This is a fairly ugly dispatch mechanism
+
         if (derivedName == "rhoU")
         {
-            added.append(derivedName);
-
-            if (!db.foundObject<volVectorField>(derivedName))
+            if (add_rhoU(derivedName))
             {
                 cleanup.append(derivedName);
-
-                db.store
-                (
-                    new volVectorField
-                    (
-                        derivedName,
-                        // rhoU = rho * U
-                        (
-                            mesh_.lookupObject<volScalarField>("rho")
-                          * mesh_.lookupObject<volVectorField>("U")
-                        )
-                    )
-                );
             }
+            added.append(derivedName);
         }
         else if (derivedName == "pTotal")
         {
-            added.append(derivedName);
-
-            if (!db.foundObject<volScalarField>(derivedName))
+            if (add_pTotal(derivedName))
             {
                 cleanup.append(derivedName);
-
-                const volScalarField& p =
-                    mesh_.lookupObject<volScalarField>("p");
-
-                if (p.dimensions() == dimPressure)
-                {
-                    db.store
-                    (
-                        new volScalarField
-                        (
-                            derivedName,
-                            // pTotal = p + rho U^2 / 2
-                            (
-                                p
-                              + 0.5
-                              * mesh_.lookupObject<volScalarField>("rho")
-                              * magSqr
-                                (
-                                    mesh_.lookupObject<volVectorField>("U")
-                                )
-                            )
-                        )
-                    );
-                }
-                else
-                {
-                    db.store
-                    (
-                        new volScalarField
-                        (
-                            derivedName,
-                            // pTotal = p + U^2 / 2
-                            (
-                                p
-                              + 0.5
-                              * magSqr
-                                (
-                                    mesh_.lookupObject<volVectorField>("U")
-                                )
-                            )
-                        )
-                    );
-                }
             }
+            added.append(derivedName);
         }
         else
         {
@@ -307,6 +340,8 @@ bool Foam::surfMeshSamplers::read(const dictionary& dict)
     fieldSelection_.clear();
     derivedNames_.clear();
 
+    rhoRef_ = dict.lookupOrDefault<scalar>("rhoRef", 1);
+
     const bool createOnRead = dict.lookupOrDefault("createOnRead", false);
 
     if (dict.found("surfaces"))
diff --git a/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.H b/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.H
index 4603240de37ea1996143dc6840c51e43709175ee..f97feb88cd7370c6f0f587dc8417ed2649d2735a 100644
--- a/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.H
+++ b/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.H
@@ -25,7 +25,8 @@ Class
     Foam::surfMeshSamplers
 
 Description
-    Set of surfaces to sample into a surfMesh/surfField.
+    Set of surfaces to sample from a volume field onto surfField that resides
+    on a surfMesh object.
 
     The execute() method is used to sample, and the write() method to write.
     It is fairly common to use for sampling only and have the write disabled.
@@ -53,6 +54,9 @@ Description
         // Optional: pre-defined derived fields to be sampled
         derived         (rhoU pTotal);
 
+        // Reference density for incompressible
+        rhoRef          1.25;
+
         // Optional: create surface immediately on read
         // The default is to create a placeholder without any faces.
         createOnRead    false;
@@ -69,6 +73,23 @@ Description
     }
     \endverbatim
 
+    Entries:
+    \table
+        Property | Description                             | Required | Default
+        type     | surfMeshes                              | yes      |
+        surfaces | the list of sample surfaces             | recommended |
+        fields   | word/regex list of fields to sampled    | yes      |
+        derived  | additional derived fields pTotal/rhoU   | no       |
+        rhoRef | reference density for derived fields (incompressible) | no | 1
+        sampleScheme | scheme to obtain face centre value  | no       | cell
+        createOnRead | Create surface immediately on read  | no       | false;
+    \endtable
+
+Note
+    The default is to create a placeholder surMesh without any faces on
+    construction. This behaviour can be changed by the createOnRead option.
+    For incompressible cases, \c rhoRef can be specified for use in the
+    derived quantities. The default is 1.
 
 See also
     Foam::sampledSurfaces
@@ -93,13 +114,13 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
+// Forward declarations
 class Time;
 class fvMesh;
 class dictionary;
 
 /*---------------------------------------------------------------------------*\
-                       Class surfMeshSamplers Declaration
+                      Class surfMeshSamplers Declaration
 \*---------------------------------------------------------------------------*/
 
 class surfMeshSamplers
@@ -119,7 +140,7 @@ class surfMeshSamplers
         const fvMesh& mesh_;
 
 
-      // Read from dictonary
+    // Read from dictionary
 
         //- Names of fields to sample
         wordRes fieldSelection_;
@@ -130,6 +151,9 @@ class surfMeshSamplers
         //- Sample scheme to obtain face values
         word sampleScheme_;
 
+        //- Reference density (to convert from kinematic to static pressure)
+        scalar rhoRef_;
+
 
     // Private Member Functions
 
@@ -140,6 +164,16 @@ class surfMeshSamplers
             const UList<word>& names
         );
 
+
+        //- Hard-coded derived field (rho * U)
+        //  \return true if field did not previously exist
+        bool add_rhoU(const word& derivedName);
+
+        //- Hard-coded derived field (p + 1/2 * rho * U)
+        //  \return true if field did not previously exist
+        bool add_pTotal(const word& derivedName);
+
+
         //- Access the sampling surfaces
         inline const PtrList<surfMeshSample>& surfaces() const
         {
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fieldTransfer b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fieldTransfer
index c953262b8c514254cd0a7925cc0d00922ab9f8c0..3ccdab376c2100e6834fa925485fabaacd896437 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fieldTransfer
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fieldTransfer
@@ -14,7 +14,7 @@ fieldTransfer
     executeInterval 1;
 
     fields      (p rho U T);
-    derived     (rhoU);
+    derived     (rhoU pTotal);
 
     _plane
     {
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling b/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
index 4fd482ff4fd4e0ee98bef325e50db028a8e8d05c..964e2c6bc02f0946f6cffde35e1ab408f9a9f1f7 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
@@ -16,6 +16,18 @@ massflow
     fields          ( rhoU );
 }
 
+areaAverage
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane1;
+
+    operation       weightedAreaAverage;
+    weightField     rhoU;
+    fields          ( p pTotal );
+}
+
 areaIntegrate
 {
     ${__surfaceFieldValue}
diff --git a/tutorials/incompressible/simpleFoam/squareBend/0/U b/tutorials/incompressible/simpleFoam/squareBend/0/U
new file mode 100644
index 0000000000000000000000000000000000000000..4bae1cccf4a78daced48bb8f552cc31401111eab
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/0/U
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    Default_Boundary_Region
+    {
+        type            noSlip;
+    }
+    inlet
+    {
+        type  fixedValue;
+        value uniform (10 0 0);
+    }
+    outlet
+    {
+        type            inletOutlet;
+        value           $internalField;
+        inletValue      $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/0/epsilon b/tutorials/incompressible/simpleFoam/squareBend/0/epsilon
new file mode 100644
index 0000000000000000000000000000000000000000..4b75d7acc38265d3ac27a7dc58e1b1b74c83a1fa
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/0/epsilon
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 200;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           $internalField;
+    }
+    inlet
+    {
+        type            turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.005;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        value           $internalField;
+        inletValue      $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/0/k b/tutorials/incompressible/simpleFoam/squareBend/0/k
new file mode 100644
index 0000000000000000000000000000000000000000..8a9421b335cfd6c8b34c32c4bfc49ef1dea66553
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/0/k
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            kqRWallFunction;
+        value           $internalField;
+    }
+    inlet
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.05;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        value           $internalField;
+        inletValue      $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/0/nut b/tutorials/incompressible/simpleFoam/squareBend/0/nut
new file mode 100644
index 0000000000000000000000000000000000000000..66f50ded1a9d5e00084613e4212ac800e19fe036
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/0/nut
@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           $internalField;
+    }
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/0/p b/tutorials/incompressible/simpleFoam/squareBend/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..2a620468fbe1610056a1253cf36da54afa45ff6b
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/0/p
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    Default_Boundary_Region
+    {
+        type            zeroGradient;
+    }
+    inlet
+    {
+        type            zeroGradient;
+        refValue        $internalField;
+        refGradient     uniform 0;
+        valueFraction   uniform 0.3;
+    }
+    outlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    ".*"
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/constant/transportProperties b/tutorials/incompressible/simpleFoam/squareBend/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..c4e6700e33209b565645cb3d82471d590265a4b9
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/constant/transportProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              1.5e-05;
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/constant/turbulenceProperties b/tutorials/incompressible/simpleFoam/squareBend/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..8fe9b5ed38d3ed562c1cdf8fcc71b0a6ce10c8bc
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/constant/turbulenceProperties
@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RAS;
+
+RAS
+{
+    RASModel        kEpsilon;
+
+    turbulence      on;
+
+    printCoeffs     on;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/blockMeshDict b/tutorials/incompressible/simpleFoam/squareBend/system/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..4980330b276ee58150f1593f0ba9b6b23c383520
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/blockMeshDict
@@ -0,0 +1,127 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   0.001;
+
+vertices
+(
+    // front-plane: z = +25mm
+    // inlet region
+    (   -50  25   25)           // pt 0
+    (     0  25   25)           // pt 1
+    (   -50  75   25)           // pt 2
+    (     0  75   25)           // pt 3
+    // outlet region
+    (  -500 -75   25)           // pt 4
+    (     0 -75   25)           // pt 5
+    (  -500 -25   25)           // pt 6
+    (     0 -25   25)           // pt 7
+    // bend mid-points
+    (    25   0   25)           // pt 8
+    (    75   0   25)           // pt 9
+    // back-plane: z = -25mm
+    // inlet region
+    (   -50  25   -25)          // pt 0 + 10
+    (     0  25   -25)          // pt 1 + 10
+    (   -50  75   -25)          // pt 2 + 10
+    (     0  75   -25)          // pt 3 + 10
+    // outlet region
+    (  -500 -75   -25)          // pt 4 + 10
+    (     0 -75   -25)          // pt 5 + 10
+    (  -500 -25   -25)          // pt 7 + 10
+    (     0 -25   -25)          // pt 8 + 10
+    // bend mid-points
+    (    25   0   -25)          // pt 8 + 10
+    (    75   0   -25)          // pt 9 + 10
+);
+
+blocks
+(
+    hex (0 1 11 10  2 3 13 12) inlet  ( 20 20 20)  simpleGrading (1 1 1)
+    hex (4 5 15 14  6 7 17 16) outlet (200 20 20)  simpleGrading (1 1 1)
+
+    hex (1 8 18 11  3 9 19 13) bend1  ( 30 20 20)  simpleGrading (1 1 1)
+    hex (5 9 19 15  7 8 18 17) bend2  ( 30 20 20)  simpleGrading (1 1 1)
+);
+
+edges
+(
+   // block 2
+   arc  1  8  ( 17.678  17.678  25)
+   arc 11 18  ( 17.678  17.678 -25)
+   arc  3  9  ( 53.033  53.033  25)
+   arc 13 19  ( 53.033  53.033 -25)
+   // block 3
+   arc  7  8  ( 17.678  -17.678  25)
+   arc 17 18  ( 17.678  -17.678 -25)
+   arc  5  9  ( 53.033  -53.033  25)
+   arc 15 19  ( 53.033  -53.033 -25)
+);
+
+boundary
+(
+        // is there no way of defining all my 'defaultFaces' to be 'wall'?
+    Default_Boundary_Region
+    {
+        type wall;
+        faces
+        (
+            // block0
+            ( 0 1 3 2 )
+            ( 11 10 12 13 )
+            ( 0 10 11 1 )
+            ( 2 3 13 12 )
+            // block1
+            ( 4 5 7 6 )
+            ( 15 14 16 17 )
+            ( 4 14 15 5 )
+            ( 6 7 17 16 )
+            // block2
+            ( 1 8 9 3 )
+            ( 18 11 13 19 )
+            ( 3 9 19 13 )
+            ( 1 11 18 8 )
+            // block3
+            ( 5 9 8 7 )
+            ( 19 15 17 18 )
+            ( 5 15 19 9 )
+            ( 7 8 18 17 )
+        );
+    }
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 2 12 10)
+        );
+    }
+
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (4 6 16 14)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/controlDict b/tutorials/incompressible/simpleFoam/squareBend/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..d39fbeca2e89297a7f337e03f65806caacd11a1a
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/controlDict
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     simpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         500;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   50;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+graphFormat     raw;
+
+runTimeModifiable true;
+
+#include "sampleControls"
+
+functions
+{
+    #include "sampling"
+    #include "samplingDebug"
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/decomposeParDict b/tutorials/incompressible/simpleFoam/squareBend/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..b27efc23b92a517e54a65ca512f095be7728810d
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/decomposeParDict
@@ -0,0 +1,34 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 8;
+
+method          hierarchical;
+
+simpleCoeffs
+{
+    n               (8 1 1);
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               (4 2 1);
+    delta           0.001;
+    order           xyz;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/fieldTransfer b/tutorials/incompressible/simpleFoam/squareBend/system/fieldTransfer
new file mode 100644
index 0000000000000000000000000000000000000000..9c596eb7803c62c4ecf42695222c6c6b4bae3067
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/fieldTransfer
@@ -0,0 +1,65 @@
+// -*- 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);
+    derived     (rhoU pTotal);
+
+    rhoRef      1.25;
+
+    _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/fvSchemes b/tutorials/incompressible/simpleFoam/squareBend/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..c0bf1c9a597be31805763310f2ac8709a53f75ad
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/fvSchemes
@@ -0,0 +1,72 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss limitedLinearV 1;
+    div(phi,k)      Gauss limitedLinear 1;;
+    div(phi,epsilon) Gauss limitedLinear 1;;
+    div(phi,R)      Gauss limitedLinear 1;;
+    div(R)          Gauss limitedLinear 1;;
+    div(phi,omega)  Gauss limitedLinear 1;
+    div(phid,p)     Gauss limitedLinear 1;
+    div(phi,K)      Gauss limitedLinear 1;
+    div(phi,e)      Gauss limitedLinear 1;
+    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+    div((nuEff*dev2(T(grad(U)))))   Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected ;
+}
+
+oversetInterpolation
+{
+    method          inverseDistance;
+}
+
+fluxRequired
+{
+    default         no;
+    pcorr           ;
+    p               ;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/fvSolution b/tutorials/incompressible/simpleFoam/squareBend/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..3cea14579d2a96a9264a047009475dbf7a2d83e7
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/fvSolution
@@ -0,0 +1,86 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          GAMG;
+        smoother        GaussSeidel;
+        tolerance       1e-7;
+        relTol          0.01;
+    }
+
+    Phi
+    {
+        $p;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       1e-8;
+        relTol          0.1;
+        nSweeps         1;
+    }
+
+    k
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       1e-8;
+        relTol          0.1;
+        nSweeps         1;
+    }
+
+    "(epsilon|omega)"
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       1e-8;
+        relTol          0.1;
+        nSweeps         1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    consistent yes;
+}
+
+potentialFlow
+{
+    nNonOrthogonalCorrectors 10;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        U               0.9;
+        k               0.7;
+        "(epsilon|omega)" 0.7;
+    }
+}
+
+cache
+{
+    grad(U);
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/sampleControls b/tutorials/incompressible/simpleFoam/squareBend/system/sampleControls
new file mode 100644
index 0000000000000000000000000000000000000000..23e1cf4246d806e2ebed072dd335b6c0928ecd36
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/sampleControls
@@ -0,0 +1,40 @@
+// -*- 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
new file mode 100644
index 0000000000000000000000000000000000000000..9feabd4078954e047a06184d3dbd6b19906d0383
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/sampling
@@ -0,0 +1,110 @@
+// -*- C++ -*-
+
+// ************************************************************************* //
+
+#include "fieldTransfer"
+
+massflow
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane1;
+
+    operation       areaNormalIntegrate;
+
+    fields          ( rhoU );
+}
+
+areaAverage
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane1;
+
+    operation       weightedAreaAverage;
+    weightField     rhoU;
+    fields          ( p pTotal );
+}
+
+areaIntegrate
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane1;
+
+    operation       weightedAreaIntegrate;
+    weightField     rhoU;
+    fields          ( p pTotal );
+}
+
+// Inflow uniformity
+UI1
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane1;
+
+    operation       uniformity;
+    fields          ( U p );
+}
+
+
+// Uniformity after the bend
+UI2
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane2;
+
+    operation       uniformity;
+    fields          ( U p );
+}
+
+
+// Uniformity on sampled surface
+UI3
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane3;
+
+    operation       uniformity;
+    fields          ( U p );
+}
+
+
+// rhoU-weighted uniformity, including weighting U too (weird but possible)
+rhoU_UI1
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane1;
+
+    operation       weightedUniformity;
+    weightField     rhoU;
+    fields          ( p U rhoU );
+}
+
+
+// rhoU-weighted uniformity
+rhoU_UI2
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane2;
+
+    operation       weightedUniformity;
+    weightField     rhoU;
+    fields          ( p U rhoU );
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/samplingDebug b/tutorials/incompressible/simpleFoam/squareBend/system/samplingDebug
new file mode 100644
index 0000000000000000000000000000000000000000..52a15b82f6ce1f238acb66119f7fcbabec4ce0cc
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/samplingDebug
@@ -0,0 +1,53 @@
+// -*- C++ -*-
+
+// ************************************************************************* //
+
+debug
+{
+    type    surfaces;
+    libs    ("libsampling.so");
+    log     true;
+    writeControl  timeStep;
+    writeInterval 1;
+
+    fields  (U);
+
+    sampleScheme  cellPoint;
+    interpolationScheme cellPoint;
+    // surfaceFormat ensight;
+    surfaceFormat vtk;
+
+    formatOptions
+    {
+        ensight
+        {
+            collateTimes  true;
+            // collateTimes  false;
+        }
+    }
+
+    surfaces
+    (
+        angledPlane
+        {
+            type        distanceSurface;
+            distance    0;
+            signed      true;
+            regularise  true;
+            surfaceType triSurfaceMesh;
+            surfaceName angledPlane.obj;
+        }
+
+        iso
+        {
+            type        isoSurface;
+            isoField    p;
+            isoValue    1e5;
+            regularise  true;
+            interpolate true;
+        }
+    );
+}
+
+
+// ************************************************************************* //