diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index 18ec9b28dbf1d0b667be96ac98e4c59eef182636..9e1d2eec8704dac016ae73b1039771222f8687ae 100644
--- a/src/functionObjects/field/Make/files
+++ b/src/functionObjects/field/Make/files
@@ -4,6 +4,8 @@ columnAverage/columnAverage.C
 
 continuityError/continuityError.C
 
+derivedFields/derivedFields.C
+
 fieldAverage/fieldAverage.C
 fieldAverage/fieldAverageItem/fieldAverageItem.C
 fieldAverage/fieldAverageItem/fieldAverageItemIO.C
diff --git a/src/functionObjects/field/derivedFields/derivedFields.C b/src/functionObjects/field/derivedFields/derivedFields.C
new file mode 100644
index 0000000000000000000000000000000000000000..d3c998fd0b5e9e6798e1383cd59dbe7c680ca981
--- /dev/null
+++ b/src/functionObjects/field/derivedFields/derivedFields.C
@@ -0,0 +1,348 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "derivedFields.H"
+#include "volFields.H"
+#include "dictionary.H"
+#include "Time.H"
+#include "mapPolyMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(derivedFields, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        derivedFields,
+        dictionary
+    );
+}
+}
+
+
+const Foam::Enum
+<
+    Foam::functionObjects::derivedFields::derivedType
+>
+Foam::functionObjects::derivedFields::knownNames
+({
+    { derivedType::NONE , "none" },
+    { derivedType::MASS_FLUX , "rhoU" },
+    { derivedType::TOTAL_PRESSURE , "pTotal" },
+});
+
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+static bool calc_rhoU
+(
+    const fvMesh& mesh,
+    const word& derivedName,
+    const scalar rhoRef
+)
+{
+    // rhoU = rho * U
+
+    const auto* rhoPtr = mesh.findObject<volScalarField>("rho");
+    const volVectorField& U = mesh.lookupObject<volVectorField>("U");
+
+    volVectorField* result = mesh.getObjectPtr<volVectorField>(derivedName);
+
+    const bool isNew = !result;
+
+    if (!result)
+    {
+        result = new volVectorField
+        (
+            IOobject
+            (
+                derivedName,
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                true
+            ),
+            mesh,
+            (dimDensity * dimVelocity)
+        );
+
+        result->store();
+    }
+
+    if (rhoPtr)
+    {
+        const auto& rho = *rhoPtr;
+
+        *result = (rho * U);
+    }
+    else
+    {
+        const dimensionedScalar rho("rho", dimDensity, rhoRef);
+
+        *result = (rho * U);
+    }
+
+    return isNew;
+}
+
+
+static bool calc_pTotal
+(
+    const fvMesh& mesh,
+    const word& derivedName,
+    const scalar rhoRef
+)
+{
+    // pTotal = p + rho * U^2 / 2
+
+    const auto* rhoPtr = mesh.findObject<volScalarField>("rho");
+    const volScalarField& p = mesh.lookupObject<volScalarField>("p");
+    const volVectorField& U = mesh.lookupObject<volVectorField>("U");
+
+    volScalarField* result = mesh.getObjectPtr<volScalarField>(derivedName);
+
+    const bool isNew = !result;
+
+    if (!result)
+    {
+        result = new volScalarField
+        (
+            IOobject
+            (
+                derivedName,
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                true
+            ),
+            mesh,
+            dimPressure
+        );
+
+        result->store();
+    }
+
+    if (rhoPtr)
+    {
+        const auto& rho = *rhoPtr;
+
+        *result = (p + 0.5 * rho * magSqr(U));
+    }
+    else
+    {
+        const dimensionedScalar rho("rho", dimDensity, rhoRef);
+
+        *result = (rho * (p + 0.5 * magSqr(U)));
+    }
+
+    return isNew;
+}
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::derivedFields::derivedFields
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    derivedTypes_(),
+    rhoRef_(1.0)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::derivedFields::read(const dictionary& dict)
+{
+    fvMeshFunctionObject::read(dict),
+
+    rhoRef_ = dict.lookupOrDefault<scalar>("rhoRef", 1);
+
+    wordList derivedNames(dict.get<wordList>("derived"));
+
+    derivedTypes_.resize(derivedNames.size());
+
+    label nbad = 0, ngood = 0;
+
+    for (const word& key : derivedNames)
+    {
+        derivedTypes_[ngood] = knownNames.get(key, derivedType::UNKNOWN);
+
+        switch (derivedTypes_[ngood])
+        {
+            case derivedType::NONE:
+                break;
+
+            case derivedType::UNKNOWN:
+            {
+                derivedNames[nbad++] = key;
+                break;
+            }
+
+            default:
+            {
+                ++ngood;
+                break;
+            }
+        }
+    }
+
+    if (nbad)
+    {
+        WarningInFunction
+            << "Ignoring unknown derived names: "
+            << SubList<word>(derivedNames, nbad) << nl;
+    }
+
+    derivedTypes_.resize(ngood);
+
+    // Output the good names
+    forAll(derivedTypes_, i)
+    {
+        derivedNames[i] = knownNames[derivedTypes_[i]];
+    }
+
+    Info<< type() << " derived: "
+        << flatOutput(SubList<word>(derivedNames, ngood)) << nl;
+
+    return true;
+}
+
+
+bool Foam::functionObjects::derivedFields::execute()
+{
+    Log << type() << " calculating:";
+
+    for (const derivedType category : derivedTypes_)
+    {
+        bool isNew = false;
+
+        switch (category)
+        {
+            case derivedType::MASS_FLUX:
+            {
+                isNew = calc_rhoU(mesh_, knownNames[category], rhoRef_);
+
+                Log << "  " << knownNames[category];
+                if (isNew) Log << " (new)";
+                break;
+            }
+
+            case derivedType::TOTAL_PRESSURE:
+            {
+                isNew = calc_pTotal(mesh_, knownNames[category], rhoRef_);
+
+                Log << "  " << knownNames[category];
+                if (isNew) Log << " (new)";
+                break;
+            }
+
+            default:
+                break;
+        }
+    }
+
+    Log << nl;
+
+    return true;
+}
+
+
+bool Foam::functionObjects::derivedFields::write()
+{
+    for (const derivedType category : derivedTypes_)
+    {
+        switch (category)
+        {
+            case derivedType::NONE:
+            case derivedType::UNKNOWN:
+                break;
+
+            default:
+            {
+                const auto* ioptr =
+                    mesh_.cfindObject<regIOobject>(knownNames[category]);
+
+                if (ioptr)
+                {
+                    Log << type() << " " << name() << " write:" << nl
+                        << "    writing field " << ioptr->name() << endl;
+
+                    ioptr->write();
+                }
+                break;
+            }
+        }
+    }
+
+    return true;
+}
+
+
+void Foam::functionObjects::derivedFields::removeDerivedFields()
+{
+    for (const derivedType category : derivedTypes_)
+    {
+        mesh_.thisDb().checkOut(knownNames[category]);
+    }
+}
+
+
+void Foam::functionObjects::derivedFields::updateMesh(const mapPolyMesh& mpm)
+{
+    if (&mpm.mesh() == &mesh_)
+    {
+        removeDerivedFields();
+    }
+}
+
+
+void Foam::functionObjects::derivedFields::movePoints(const polyMesh& m)
+{
+    if (&m == &mesh_)
+    {
+        removeDerivedFields();
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/derivedFields/derivedFields.H b/src/functionObjects/field/derivedFields/derivedFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..0bf3b7e23d4dec7c4b5091802163e063c28dab44
--- /dev/null
+++ b/src/functionObjects/field/derivedFields/derivedFields.H
@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::derivedFields
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    A limited set of predefined derived fields ("rhoU", "pTotal").
+
+    \verbatim
+    derived
+    {
+        type    derivedFields;
+        libs    ("libfieldFunctionObjects.so");
+
+        fields  (rhoU pTotal);
+
+        // Optional: reference density for incompressible
+        rhoRef  1.25;
+    }
+    \endverbatim
+
+    Entries:
+    \table
+        Property | Description                              | Required | Default
+        type     | derivedFields                            | yes |
+        derived  | Derived fields (pTotal/rhoU)             | yes |
+        rhoRef   | Reference density (incompressible)       | no  | 1
+    \endtable
+
+    The known derived fields
+    \plaintable
+       rhoU          | (rho * U)
+       pTotal        | (p + 1/2 * rho * U)
+    \endplaintable
+
+SourceFiles
+    derivedFields.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_derivedFields_H
+#define functionObjects_derivedFields_H
+
+#include "fvMeshFunctionObject.H"
+#include "Enum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class derivedFields Declaration
+\*---------------------------------------------------------------------------*/
+
+class derivedFields
+:
+    public fvMeshFunctionObject
+{
+public:
+
+    // Public Enumerations
+
+        //- Derived/calculated field type
+        enum derivedType
+        {
+            NONE = 0,           //!< "none"
+            MASS_FLUX,          //!< "rhoU"
+            TOTAL_PRESSURE,     //!< "pTotal"
+            UNKNOWN
+        };
+
+        //- Known derived field types
+        static const Enum<derivedType> knownNames;
+
+
+protected:
+
+    // Read from dictionary
+
+        //- List of derived field (types) to create
+        List<derivedType> derivedTypes_;
+
+        //- Reference density (to convert from kinematic to static pressure)
+        scalar rhoRef_;
+
+
+    // Private Member Functions
+
+        //- 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);
+
+
+public:
+
+    //- Run-time type information
+    TypeName("derivedFields");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        derivedFields
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~derivedFields() = default;
+
+
+    // Member Functions
+
+        //- Remove (checkOut) derived fields from the object registry
+        void removeDerivedFields();
+
+        //- Read the data
+        virtual bool read(const dictionary& dict);
+
+        //- Calculate the derived fields
+        virtual bool execute();
+
+        //- Write derived fields
+        virtual bool write();
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh& mpm);
+
+        //- Update for mesh point-motion
+        virtual void movePoints(const polyMesh& m);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.C b/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.C
index fb7019f0980972770ab6479dc10f434c7b8389bb..eb8b8fde32279ea5c6449f50e7d055f2f6196bfa 100644
--- a/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.C
+++ b/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.C
@@ -77,93 +77,6 @@ bool Foam::surfMeshSamplers::verbose_ = false;
 //     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_.findObject<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_.findObject<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  * * * * * * * * * * * * * * //
 
@@ -177,9 +90,7 @@ Foam::surfMeshSamplers::surfMeshSamplers
     functionObjects::fvMeshFunctionObject(name, runTime, dict),
     PtrList<surfMeshSample>(),
     fieldSelection_(),
-    derivedNames_(),
-    sampleFaceScheme_(),
-    rhoRef_(1.0)
+    sampleFaceScheme_()
 {
     read(dict);
 }
@@ -195,9 +106,7 @@ Foam::surfMeshSamplers::surfMeshSamplers
     functionObjects::fvMeshFunctionObject(name, obr, dict),
     PtrList<surfMeshSample>(),
     fieldSelection_(),
-    derivedNames_(),
-    sampleFaceScheme_(),
-    rhoRef_(1.0)
+    sampleFaceScheme_()
 {
     read(dict);
 }
@@ -218,44 +127,8 @@ bool Foam::surfMeshSamplers::execute()
         return true;
     }
 
-    const objectRegistry& db = mesh_.thisDb();
-
-    // Manage derived names
-    DynamicList<word> added(derivedNames_.size());
-    DynamicList<word> cleanup(derivedNames_.size());
-
-    for (const word& derivedName : derivedNames_)
-    {
-        // This is a fairly ugly dispatch mechanism
-
-        if (derivedName == "rhoU")
-        {
-            if (add_rhoU(derivedName))
-            {
-                cleanup.append(derivedName);
-            }
-            added.append(derivedName);
-        }
-        else if (derivedName == "pTotal")
-        {
-            if (add_pTotal(derivedName))
-            {
-                cleanup.append(derivedName);
-            }
-            added.append(derivedName);
-        }
-        else
-        {
-            WarningInFunction
-                << "unknown derived name: " << derivedName << nl
-                << "Use one of 'rhoU', 'pTotal'" << nl
-                << endl;
-        }
-    }
-
-
     // The acceptable fields
-    wordHashSet acceptable(added);
+    wordHashSet acceptable;
     acceptable.insert(acceptType<scalar>());
     acceptable.insert(acceptType<vector>());
     acceptable.insert(acceptType<sphericalTensor>());
@@ -277,11 +150,7 @@ bool Foam::surfMeshSamplers::execute()
         }
     }
 
-    // Cleanup any locally introduced names
-    for (const word& fieldName : cleanup)
-    {
-        db.checkOut(fieldName);
-    }
+    clearObjects(removeFieldsOnExecute_);
 
     return true;
 }
@@ -294,27 +163,13 @@ bool Foam::surfMeshSamplers::write()
     // Doesn't bother checking which fields have been generated here
     // or elsewhere
 
-    // This could be more efficient
-    wordRes select(fieldSelection_.size() + derivedNames_.size());
-
-    label nElem = 0;
-    for (const wordRe& item : fieldSelection_)
-    {
-        select[nElem++] = item;
-    }
-    for (const auto& derivedName : derivedNames_)
-    {
-        select[nElem++] = derivedName;
-    }
-
-    // Avoid duplicate entries
-    select.uniq();
-
     for (const surfMeshSample& s : surfaces())
     {
-        s.write(select);
+        s.write(fieldSelection_);
     }
 
+    clearObjects(removeFieldsOnWrite_);
+
     return true;
 }
 
@@ -325,9 +180,11 @@ bool Foam::surfMeshSamplers::read(const dictionary& dict)
 
     PtrList<surfMeshSample>::clear();
     fieldSelection_.clear();
-    derivedNames_.clear();
+    removeFieldsOnExecute_.clear();
+    removeFieldsOnWrite_.clear();
 
-    rhoRef_ = dict.lookupOrDefault<scalar>("rhoRef", 1);
+    dict.readIfPresent("removeFieldsOnExecute", removeFieldsOnExecute_);
+    dict.readIfPresent("removeFieldsOnWrite", removeFieldsOnWrite_);
 
     const bool createOnRead = dict.lookupOrDefault("createOnRead", false);
 
@@ -393,11 +250,6 @@ bool Foam::surfMeshSamplers::read(const dictionary& dict)
         fieldSelection_.uniq();
 
         Info<< type() << " fields:  " << flatOutput(fieldSelection_) << nl;
-
-        if (dict.readIfPresent("derived", derivedNames_))
-        {
-            Info<< type() << " derived: " << flatOutput(derivedNames_) << nl;
-        }
         Info<< nl;
 
         // Ensure all surfaces and merge information are expired
@@ -431,6 +283,20 @@ bool Foam::surfMeshSamplers::read(const dictionary& dict)
         Info<< nl;
     }
 
+    if (removeFieldsOnExecute_.size())
+    {
+        Info<< type() << " Remove fields on-execute : "
+            << removeFieldsOnExecute_ << nl;
+    }
+    if (removeFieldsOnWrite_.size())
+    {
+        Info<< type() << " Remove fields on-write   :"
+            << removeFieldsOnWrite_ << nl;
+    }
+
+    // Ensure all surfaces are expired (unsampled)
+    expire();
+
     return true;
 }
 
diff --git a/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.H b/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.H
index 82426d404285135e7dfea6d1f30f0f5f979b69b9..9e44f2999aa697934cc85ce21a673a73d9c84576 100644
--- a/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.H
+++ b/src/sampling/surfMeshSample/surfMeshSamplers/surfMeshSamplers.H
@@ -51,16 +51,13 @@ Description
         // Scheme to obtain face centre value
         sampleScheme    cell;
 
-        // 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;
 
+        // Optional: remove derived fields created prior
+        removeFieldsOnExecute ( someDerivedField );
+
         surfaces
         (
             f0surf
@@ -83,13 +80,13 @@ Description
         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
+        removeFieldsOnExecute | List of fields to remove as "consumed" | no  |
+        removeFieldsOnWrite   | List of fields to remove as "consumed" | no  |
     \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
@@ -134,27 +131,18 @@ class surfMeshSamplers
         //- Names of fields to sample
         wordRes fieldSelection_;
 
-        //- Names of derived fields to create and sample
-        wordList derivedNames_;
+        //- List of fields to remove as "consumed"
+        wordList removeFieldsOnExecute_;
+
+        //- List of fields to remove as "consumed"
+        wordList removeFieldsOnWrite_;
 
         //- Sample scheme to obtain face values
         word sampleFaceScheme_;
 
-        //- Reference density (to convert from kinematic to static pressure)
-        scalar rhoRef_;
-
 
     // Private Member Functions
 
-        //- 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/gasMixing/injectorPipe/system/derivedFields b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/derivedFields
new file mode 100644
index 0000000000000000000000000000000000000000..68ae53f8678d3b29d5e574e03b788a3704d3f7b8
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/derivedFields
@@ -0,0 +1,19 @@
+// -*- C++ -*-
+
+// ************************************************************************* //
+
+// Create additional volume fields (for sampling)
+derivedFields
+{
+    type    derivedFields;
+    libs    ("libfieldFunctionObjects.so");
+    log     true;
+
+    writeControl    none;
+    executeControl  timeStep;
+    executeInterval 1;
+
+    derived (rhoU pTotal);
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fieldTransfer b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fieldTransfer
index 82e8a625a1039b5493aade63ddf678f6b788dccb..b1eae09bb3142b9925149a635e67fafe4dc09284 100644
--- a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fieldTransfer
+++ b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fieldTransfer
@@ -13,8 +13,12 @@ fieldTransfer
     executeControl  timeStep;
     executeInterval 1;
 
-    fields  (rho U tracer0);
-    derived (rhoU);
+    // Includes a rhoU derived field
+    fields  (rho U tracer0 rhoU);
+
+    // Remove derived fields we created prior
+    removeFieldsOnExecute   (pTotal);
+
 
     baseCfg
     {
diff --git a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sampling b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sampling
index a14f5bdec2c75798c9bc6460d90868f1e3bb37e3..233c5e68caf1a4d397690b16d51009fd60692080 100644
--- a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sampling
+++ b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/sampling
@@ -2,6 +2,7 @@
 
 // ************************************************************************* //
 
+#include "derivedFields"
 #include "fieldTransfer"
 #include "avg-tracer0"
 #include "sum-tracer0"
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/derivedFields b/tutorials/compressible/rhoSimpleFoam/squareBend/system/derivedFields
new file mode 100644
index 0000000000000000000000000000000000000000..e28a2a9f8c2d15ca0d3c4890bedfa84ca5a75278
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/derivedFields
@@ -0,0 +1,21 @@
+// -*- C++ -*-
+
+// ************************************************************************* //
+
+// Create additional volume fields (for sampling)
+derivedFields
+{
+    type    derivedFields;
+    libs    ("libfieldFunctionObjects.so");
+    log     true;
+
+    writeControl    none;
+    executeControl  timeStep;
+    executeInterval 1;
+
+    derived     (rhoU pTotal);
+
+    rhoRef      1.25;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fieldTransfer b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fieldTransfer
index 3ccdab376c2100e6834fa925485fabaacd896437..12f4e4276b341785d5e410a1e0ae4e2ff72fd3d9 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fieldTransfer
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fieldTransfer
@@ -14,7 +14,9 @@ fieldTransfer
     executeInterval 1;
 
     fields      (p rho U T);
-    derived     (rhoU pTotal);
+
+    // Remove derived fields we created prior
+    removeFieldsOnExecute   (rhoU pTotal);
 
     _plane
     {
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling b/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
index 964e2c6bc02f0946f6cffde35e1ab408f9a9f1f7..eeab502a7dfa6b6c0d5f12d6633646782e244f1a 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
@@ -2,6 +2,7 @@
 
 // ************************************************************************* //
 
+#include "derivedFields"
 #include "fieldTransfer"
 
 massflow
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/fieldTransfer b/tutorials/incompressible/simpleFoam/squareBend/system/fieldTransfer
index 9c596eb7803c62c4ecf42695222c6c6b4bae3067..649eccae185567de3c1d63082de2c87db0c78563 100644
--- a/tutorials/incompressible/simpleFoam/squareBend/system/fieldTransfer
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/fieldTransfer
@@ -14,9 +14,9 @@ fieldTransfer
     executeInterval 1;
 
     fields      (p U);
-    derived     (rhoU pTotal);
 
-    rhoRef      1.25;
+    // Remove derived fields we created prior
+    removeFieldsOnExecute   (rhoU pTotal);
 
     _plane
     {
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/sampling b/tutorials/incompressible/simpleFoam/squareBend/system/sampling
index 9feabd4078954e047a06184d3dbd6b19906d0383..eac6e6ce233e158662800bd144f4f5432ab5c2fd 100644
--- a/tutorials/incompressible/simpleFoam/squareBend/system/sampling
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/sampling
@@ -2,6 +2,7 @@
 
 // ************************************************************************* //
 
+#include "derivedFields"
 #include "fieldTransfer"
 
 massflow