diff --git a/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/Make/options b/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/Make/options
index 3c22d1310b..2a75a610c1 100644
--- a/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/Make/options
+++ b/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/Make/options
@@ -6,4 +6,5 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
+    -ldynamicMesh \
     -lfieldFunctionObjects
diff --git a/src/OpenFOAM/db/functionObjects/regionFunctionObject/regionFunctionObject.H b/src/OpenFOAM/db/functionObjects/regionFunctionObject/regionFunctionObject.H
index d9fed7b56b..061b2300a7 100644
--- a/src/OpenFOAM/db/functionObjects/regionFunctionObject/regionFunctionObject.H
+++ b/src/OpenFOAM/db/functionObjects/regionFunctionObject/regionFunctionObject.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -74,7 +74,6 @@ class regionFunctionObject
 :
     public stateFunctionObject
 {
-
 protected:
 
     // Protected Member Data
@@ -147,6 +146,15 @@ protected:
             bool cacheable = false
         );
 
+        //- Store the field in an optional objectRegistry under the given name
+        template<class ObjectType>
+        bool storeInDb
+        (
+            const word& fieldName,
+            const tmp<ObjectType>& tfield,
+            const objectRegistry& obr
+        );
+
         //- Write field if present in the (sub) objectRegistry
         bool writeObject(const word& fieldName);
 
diff --git a/src/OpenFOAM/db/functionObjects/regionFunctionObject/regionFunctionObjectTemplates.C b/src/OpenFOAM/db/functionObjects/regionFunctionObject/regionFunctionObjectTemplates.C
index f253d83e05..c1c0f740b6 100644
--- a/src/OpenFOAM/db/functionObjects/regionFunctionObject/regionFunctionObjectTemplates.C
+++ b/src/OpenFOAM/db/functionObjects/regionFunctionObject/regionFunctionObjectTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016 OpenFOAM Foundation
-    Copyright (C) 2016-2018 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -122,16 +122,20 @@ bool Foam::functionObjects::regionFunctionObject::store
         return false;
     }
 
-    if (fieldName.size() && foundObject<ObjectType>(fieldName))
-    {
-        const ObjectType& field = lookupObject<ObjectType>(fieldName);
+    ObjectType* fieldptr;
 
+    if
+    (
+        !fieldName.empty()
+     && (fieldptr = getObjectPtr<ObjectType>(fieldName)) != nullptr
+    )
+    {
         // If there is a result field already registered, assign to the new
         // result field. Otherwise transfer ownership of the new result field to
         // the object registry
-        if (&field != &tfield())
+        if (fieldptr != &tfield())
         {
-            const_cast<ObjectType&>(field) = tfield;
+            (*fieldptr) = tfield;
         }
         else
         {
@@ -156,4 +160,31 @@ bool Foam::functionObjects::regionFunctionObject::store
 }
 
 
+template<class ObjectType>
+bool Foam::functionObjects::regionFunctionObject::storeInDb
+(
+    const word& fieldName,
+    const tmp<ObjectType>& tfield,
+    const objectRegistry& obr
+)
+{
+    ObjectType* fieldptr;
+    if
+    (
+        !fieldName.empty()
+     && (fieldptr = obr.getObjectPtr<ObjectType>(fieldName)) != nullptr
+    )
+    {
+        (*fieldptr) = tfield;
+    }
+    else
+    {
+        tfield.ref().rename(fieldName);
+        obr.objectRegistry::store(tfield.ptr());
+    }
+
+    return true;
+}
+
+
 // ************************************************************************* //
diff --git a/src/dynamicMesh/Make/files b/src/dynamicMesh/Make/files
index 08723a0583..eedcd7fe8c 100644
--- a/src/dynamicMesh/Make/files
+++ b/src/dynamicMesh/Make/files
@@ -138,4 +138,6 @@ polyMeshFilter/polyMeshFilter.C
 pointPatchDist/externalPointEdgePoint.C
 pointPatchDist/pointPatchDist.C
 
+zoneSubSet/zoneSubSet.C
+
 LIB = $(FOAM_LIBBIN)/libdynamicMesh
diff --git a/src/dynamicMesh/zoneSubSet/zoneSubSet.C b/src/dynamicMesh/zoneSubSet/zoneSubSet.C
new file mode 100644
index 0000000000..34b60e6c2d
--- /dev/null
+++ b/src/dynamicMesh/zoneSubSet/zoneSubSet.C
@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "zoneSubSet.H"
+#include "cellBitSet.H"
+#include "haloToCell.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Detail
+{
+    defineTypeNameAndDebug(zoneSubSet, 0);
+}
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::Detail::zoneSubSet::correct()
+{
+    subsetter_.clear();
+    haloCells_.clearStorage();
+
+    if (zoneMatcher_.empty())
+    {
+        return false;
+    }
+
+    // Select named zones
+    cellBitSet selectedCells
+    (
+        subsetter_.baseMesh(),
+        subsetter_.baseMesh().cellZones().selection(zoneMatcher_)
+    );
+
+    if (debug)
+    {
+        Pout<< "Subsetting "
+            << selectedCells.addressing().count()
+            << " cells based on cellZones "
+            << flatOutput(zoneMatcher_) << endl;
+    }
+
+    if (nLayers_ > 0)
+    {
+        // Add halo layer(s)
+        haloToCell haloSource(subsetter_.baseMesh(), nLayers_);
+        haloSource.verbose(false);
+
+        // Before adding halo cells
+        haloCells_ = selectedCells.addressing();
+
+        haloSource.applyToSet(topoSetSource::ADD, selectedCells);
+
+        // Halo cells: anything new, not in the original set
+        haloCells_ ^= selectedCells.addressing();
+    }
+
+    if (debug)
+    {
+        const label nHalo = haloCells_.count();
+        const label nSubCell = selectedCells.addressing().count();
+
+        Info<< "    overall "
+            << returnReduce(nSubCell, sumOp<label>())
+            << " cells after adding " << nLayers_ << " layers with "
+            << returnReduce(nHalo, sumOp<label>())
+            << " halo cells"
+            << endl;
+    }
+
+    subsetter_.setCellSubset(selectedCells.addressing());
+
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::Detail::zoneSubSet::zoneSubSet
+(
+    const fvMesh& mesh,
+    const wordRes& zoneSelector,
+    const label nZoneLayers
+)
+:
+    subsetter_(mesh),
+    zoneMatcher_(zoneSelector),
+    nLayers_(nZoneLayers),
+    haloCells_()
+{
+    correct();
+}
+
+
+Foam::Detail::zoneSubSet::zoneSubSet
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    subsetter_(mesh),
+    zoneMatcher_(),
+    nLayers_(dict.getOrDefault<label>("nLayers", 0)),
+    haloCells_()
+{
+    dict.readIfPresent("cellZones", zoneMatcher_);
+
+    correct();
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/zoneSubSet/zoneSubSet.H b/src/dynamicMesh/zoneSubSet/zoneSubSet.H
new file mode 100644
index 0000000000..6a87af2025
--- /dev/null
+++ b/src/dynamicMesh/zoneSubSet/zoneSubSet.H
@@ -0,0 +1,204 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::zoneSubSet
+
+Description
+    Intermediate tool for handling \c cellZones
+    for function objects (e.g. \c momentumError)
+    wherein the sub-mesh option is available.
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    <functionObject>
+    {
+        // Mandatory and optional entries
+        ...
+
+        // Inherited entries
+        cellZones         ( <cellZone1> <cellZone2> ... <cellZoneN> );
+        nLayers           <label>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                             | Type  | Reqd | Deflt
+      cellZones | Names of cellZones                    | wordRes | no   | -
+      nLayers   | Number of cell layers around cellZones  | label | no   | 0
+    \endtable
+
+See also
+  - Foam::functionObjects::momentumError
+  - Foam::functionObjects::div
+
+SourceFiles
+    zoneSubSet.C
+    zoneSubSetTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef zoneSubSet_H
+#define zoneSubSet_H
+
+#include "fvMeshSubset.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Detail
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class zoneSubSet Declaration
+\*---------------------------------------------------------------------------*/
+
+class zoneSubSet
+{
+    // Private Data
+
+        //- Subsetting engine
+        fvMeshSubset subsetter_;
+
+        //- Matcher for zones
+        wordRes zoneMatcher_;
+
+        //- Number of layers around zones
+        label nLayers_;
+
+        //- The halo cells
+        bitSet haloCells_;
+
+
+    // Private Members Functions
+
+        //- Initialise the sub-mesh
+        bool correct();
+
+        //- No copy construct
+        zoneSubSet(const zoneSubSet&) = delete;
+
+        //- No copy assignment
+        void operator=(const zoneSubSet&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("zoneSubSet");
+
+
+    // Constructors
+
+        //- Construct from components
+        zoneSubSet
+        (
+            const fvMesh& mesh,
+            const wordRes& zoneSelector,
+            const label nZoneLayers = 0
+        );
+
+        //- Construct from mesh and dictionary
+        zoneSubSet
+        (
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~zoneSubSet() = default;
+
+
+    // Member Functions
+
+    // Access
+
+        //- The entire base mesh
+        const fvMesh& baseMesh() const noexcept
+        {
+            return subsetter_.baseMesh();
+        }
+
+        //- The mesh subsetter
+        const fvMeshSubset& subsetter() const
+        {
+            return subsetter_;
+        }
+
+        //- The mesh subsetter
+        fvMeshSubset& subsetter()
+        {
+            return subsetter_;
+        }
+
+        //- True if sub-mesh is expected to be used
+        bool useSubMesh() const noexcept
+        {
+            return !zoneMatcher_.empty();
+        }
+
+        //- Return the current zones selector
+        const wordRes& zones() const noexcept
+        {
+            return zoneMatcher_;
+        }
+
+        //- Return number of layers around cell zones
+        label nLayers() const noexcept
+        {
+            return nLayers_;
+        }
+
+
+    // Fields
+
+        //- Map from the sub-mesh to original cell zones
+        template<class Type>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> mapToZone
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& subVolField
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Detail
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "zoneSubSetTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#endif
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/zoneSubSet/zoneSubSetTemplates.C b/src/dynamicMesh/zoneSubSet/zoneSubSetTemplates.C
new file mode 100644
index 0000000000..9cfb015d97
--- /dev/null
+++ b/src/dynamicMesh/zoneSubSet/zoneSubSetTemplates.C
@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "volFieldsFwd.H"
+
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
+Foam::Detail::zoneSubSet::mapToZone
+(
+    const GeometricField<Type, fvPatchField, volMesh>& subVolField
+) const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> volFieldType;
+
+    // The full-mesh field, with zero for non-zoned areas
+    auto tcellZonesField = volFieldType::New
+    (
+        subVolField.name(),
+        subsetter_.baseMesh(),
+        dimensioned<Type>(subVolField.dimensions())
+    );
+    auto& cellZonesField = tcellZonesField.ref();
+
+
+    // Map field in global mesh on original zones
+    // - without any halo cells
+
+    const labelList& cellMap = subsetter_.cellMap();
+
+    forAll(cellMap, subCelli)
+    {
+        const label celli = cellMap[subCelli];
+
+        if (!haloCells_.test(celli))
+        {
+            cellZonesField[celli] = subVolField[subCelli];
+        }
+    }
+
+    return tcellZonesField;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/Make/options b/src/functionObjects/field/Make/options
index 93f082e386..a99cf3db84 100644
--- a/src/functionObjects/field/Make/options
+++ b/src/functionObjects/field/Make/options
@@ -4,6 +4,7 @@ EXE_INC = \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
@@ -27,6 +28,7 @@ LIB_LIBS = \
     -lfileFormats \
     -lsurfMesh \
     -lmeshTools \
+    -ldynamicMesh \
     -lsampling \
     -llagrangian \
     -ldistributionModels \
diff --git a/src/functionObjects/field/div/div.C b/src/functionObjects/field/div/div.C
index 51fc391155..47f084deb0 100644
--- a/src/functionObjects/field/div/div.C
+++ b/src/functionObjects/field/div/div.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,12 +47,26 @@ namespace functionObjects
 
 bool Foam::functionObjects::div::calc()
 {
-    bool processed = false;
+    return
+    (
+        calcDiv<surfaceScalarField>()
+     || calcDiv<volVectorField>()
+    );
+}
 
-    processed = processed || calcDiv<surfaceScalarField>();
-    processed = processed || calcDiv<volVectorField>();
 
-    return processed;
+bool Foam::functionObjects::div::write()
+{
+    if (zoneSubSetPtr_)
+    {
+        return
+        (
+            writeField<scalar>()
+         || writeField<vector>()
+        );
+    }
+
+    return writeObject(resultName_);
 }
 
 
diff --git a/src/functionObjects/field/div/div.H b/src/functionObjects/field/div/div.H
index 74dcbe4052..0f63d750e2 100644
--- a/src/functionObjects/field/div/div.H
+++ b/src/functionObjects/field/div/div.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -53,15 +53,15 @@ Usage
         // Mandatory (inherited) entry (runtime modifiable)
         field           \<field\>;
 
-        // Optional (inherited) entries
+        // Inherited entries
         ...
     }
     \endverbatim
 
     where the entries mean:
     \table
-      Property     | Description                        | Type | Req'd | Dflt
-      type         | Type name: add                     | word |  yes  | -
+      Property     | Description                        | Type | Reqd  | Deflt
+      type         | Type name: div                     | word |  yes  | -
       libs         | Library name: fieldFunctionObjects | word |  yes  | -
       field        | Name of the operand field          | word |  yes  | -
     \endtable
@@ -69,6 +69,7 @@ Usage
     The inherited entries are elaborated in:
      - \link functionObject.H \endlink
      - \link fieldExpression.H \endlink
+     - \link zoneSubSet.H \endlink
 
     Minimal example by using the \c postProcess utility:
     \verbatim
@@ -76,8 +77,10 @@ Usage
     \endverbatim
 
 Note
-    To execute \c div function object on an input \<field\>, a numerical scheme
+  - To execute \c div function object on an input \<field\>, a numerical scheme
     should be defined for \c div(\<field\>) in \c system/fvSchemes.divSchemes.
+  - Optionally the user can specify \c cellZones to create a sub-mesh for the
+    \c div calculation.
 
 See also
     - Foam::functionObject
@@ -121,6 +124,10 @@ class div
         //- Calculate the divergence field and return true if successful
         virtual bool calc();
 
+        //- Helper function for writing submesh fields
+        template<class Type>
+        bool writeField();
+
 
 public:
 
@@ -147,6 +154,12 @@ public:
 
     //- Destructor
     virtual ~div() = default;
+
+
+    // Member Functions
+
+        //- Write the result field
+        virtual bool write();
 };
 
 
diff --git a/src/functionObjects/field/div/divTemplates.C b/src/functionObjects/field/div/divTemplates.C
index 40f486bf81..cf50acfb62 100644
--- a/src/functionObjects/field/div/divTemplates.C
+++ b/src/functionObjects/field/div/divTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,15 +33,50 @@ License
 template<class FieldType>
 bool Foam::functionObjects::div::calcDiv()
 {
-    if (foundObject<FieldType>(fieldName_, false))
+    const auto* fieldptr = cfindObject<FieldType>(fieldName_);
+
+    if (!fieldptr)
     {
-        return store
+        return false;
+    }
+
+    if (zoneSubSetPtr_)
+    {
+        const fvMeshSubset& subsetter = zoneSubSetPtr_->subsetter();
+
+        return storeInDb
         (
             resultName_,
-            fvc::div(lookupObject<FieldType>(fieldName_))
+            fvc::div(subsetter.interpolate(*fieldptr, false)),
+            subsetter.subMesh().thisDb()
         );
     }
 
+
+    return store
+    (
+        resultName_,
+        fvc::div(*fieldptr)
+    );
+}
+
+
+template<class Type>
+bool Foam::functionObjects::div::writeField()
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> volFieldType;
+
+    const fvMesh& subMesh = zoneSubSetPtr_->subsetter().subMesh();
+    const auto* fieldptr = subMesh.findObject<volFieldType>(resultName_);
+
+    if (fieldptr)
+    {
+        tmp<volFieldType> tfield = zoneSubSetPtr_->mapToZone<Type>(*fieldptr);
+        tfield().write();
+
+        return true;
+    }
+
     return false;
 }
 
diff --git a/src/functionObjects/field/fieldExpression/fieldExpression.C b/src/functionObjects/field/fieldExpression/fieldExpression.C
index df4284db04..53b22aba1e 100644
--- a/src/functionObjects/field/fieldExpression/fieldExpression.C
+++ b/src/functionObjects/field/fieldExpression/fieldExpression.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -90,16 +90,28 @@ Foam::functionObjects::fieldExpression::fieldExpression
 
 bool Foam::functionObjects::fieldExpression::read(const dictionary& dict)
 {
-    fvMeshFunctionObject::read(dict);
-
-    if (fieldName_.empty() || dict.found("field"))
+    if (fvMeshFunctionObject::read(dict))
     {
-        dict.readEntry("field", fieldName_);
+        if (fieldName_.empty() || dict.found("field"))
+        {
+            dict.readEntry("field", fieldName_);
+        }
+
+        dict.readIfPresent("result", resultName_);
+
+        if (dict.found("cellZones"))
+        {
+            zoneSubSetPtr_.reset(new Detail::zoneSubSet(mesh_, dict));
+        }
+        else
+        {
+            zoneSubSetPtr_.reset(nullptr);
+        }
+
+        return true;
     }
 
-    dict.readIfPresent("result", resultName_);
-
-    return true;
+    return false;
 }
 
 
diff --git a/src/functionObjects/field/fieldExpression/fieldExpression.H b/src/functionObjects/field/fieldExpression/fieldExpression.H
index ae6df7703a..ed9977ea0a 100644
--- a/src/functionObjects/field/fieldExpression/fieldExpression.H
+++ b/src/functionObjects/field/fieldExpression/fieldExpression.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,22 +37,26 @@ Description
 Usage
     Minimal example by using \c system/controlDict.functions:
     \verbatim
-    <userDefinedSubDictName1>
+    <functionObject>
     {
-        // Mandatory and other optional entries
+        // Mandatory and optional entries
         ...
 
-        // Optional (inherited) entries (runtime modifiable)
+        // Inherited entries
         field             <field>;
         result            <fieldResult>;
+        cellZones         ( <cellZone1> <cellZone2> ... <cellZoneN> );
+        nLayers           <label>;
     }
     \endverbatim
 
     where the entries mean:
     \table
-      Property  | Description               | Type | Req'd | Dflt
-      field     | Name of the operand field | word | yes   | -
-      result    | Name of the output field  | word | no    | \<FO\>(\<field\>)
+      Property  | Description                             | Type  | Reqd | Deflt
+      field     | Name of the operand field               | word  | yes  | -
+      result    | Name of the output field  | word | no   | \<FO\>(\<field\>)
+      cellZones | Names of cellZones                    | wordRes | no   | -
+      nLayers   | Number of cell layers around cellZones  | label | no   | 0
     \endtable
 
 See also
@@ -71,6 +75,7 @@ SourceFiles
 
 #include "fvMeshFunctionObject.H"
 #include "volFieldsFwd.H"
+#include "zoneSubSet.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -89,7 +94,7 @@ class fieldExpression
 {
 protected:
 
-    // Protected Member Data
+    // Protected Data
 
         //- Name of field to process
         word fieldName_;
@@ -97,6 +102,9 @@ protected:
         //- Name of result field
         word resultName_;
 
+        //- Sub-set mesh
+        autoPtr<Detail::zoneSubSet> zoneSubSetPtr_;
+
 
     // Protected Member Functions
 
diff --git a/src/functionObjects/field/momentumError/momentumError.C b/src/functionObjects/field/momentumError/momentumError.C
index 00920f545c..0761c4a105 100644
--- a/src/functionObjects/field/momentumError/momentumError.C
+++ b/src/functionObjects/field/momentumError/momentumError.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,6 +54,9 @@ Foam::functionObjects::momentumError::divDevRhoReff()
     typedef compressible::turbulenceModel cmpTurbModel;
     typedef incompressible::turbulenceModel icoTurbModel;
 
+    const auto& U = lookupObject<volVectorField>(UName_);
+    tmp<volVectorField> tU(U);
+
     {
         auto* turb = findObject<cmpTurbModel>
         (
@@ -62,19 +65,31 @@ Foam::functionObjects::momentumError::divDevRhoReff()
 
         if (turb)
         {
+            tmp<volScalarField> trho = turb->rho();
+            tmp<volScalarField> tnuEff = turb->nuEff();
+
+            if (zoneSubSetPtr_)
+            {
+                const fvMeshSubset& subsetter = zoneSubSetPtr_->subsetter();
+
+                tU = subsetter.interpolate(U, false);
+                trho = subsetter.interpolate(turb->rho(), false);
+                tnuEff = subsetter.interpolate(turb->nuEff()(), false);
+            }
+
             return tmp<volVectorField>::New
             (
                 "divDevRhoReff",
-              - fvc::div
+                - fvc::div
                 (
-                    (turb->rho()*turb->nuEff())
-                   *dev2(T(fvc::grad(turb->U()))),
-                   "div(((rho*nuEff)*dev2(T(grad(U)))))"
+                    (trho()*tnuEff())
+                    *dev2(T(fvc::grad(tU()))),
+                    "div(((rho*nuEff)*dev2(T(grad(U)))))"
                 )
-              - fvc::laplacian
+                - fvc::laplacian
                 (
-                    turb->rho()*turb->nuEff(),
-                    turb->U(),
+                    trho()*tnuEff(),
+                    tU(),
                     "laplacian(nuEff,U)"
                 )
             );
@@ -89,18 +104,25 @@ Foam::functionObjects::momentumError::divDevRhoReff()
 
         if (turb)
         {
+            tmp<volScalarField> tnuEff = turb->nuEff();
+
+            if (zoneSubSetPtr_)
+            {
+                const fvMeshSubset& subsetter = zoneSubSetPtr_->subsetter();
+
+                tU = subsetter.interpolate(U, false);
+                tnuEff = subsetter.interpolate(turb->nuEff()(), false);
+            }
+
             return tmp<volVectorField>::New
             (
-                "divDevReff",
-              - fvc::div
+                "divDevRhoReff",
+                - fvc::div
                 (
-                    (turb->nuEff())*dev2(T(fvc::grad(turb->U()))),
-                    "div((nuEff*dev2(T(grad(U)))))"
-                )
-              - fvc::laplacian
-                (
-                    turb->nuEff(), turb->U(), "laplacian(nuEff,U)"
+                    tnuEff()*dev2(T(fvc::grad(tU()))),
+                    "div(((nuEff)*dev2(T(grad(U)))))"
                 )
+                - fvc::laplacian(tnuEff(), tU(), "laplacian(nuEff,U)")
             );
         }
      }
@@ -125,69 +147,139 @@ Foam::functionObjects::momentumError::momentumError
 {
     read(dict);
 
-    const surfaceScalarField& phi =
-        lookupObject<surfaceScalarField>(phiName_);
+    const auto& phi =lookupObject<surfaceScalarField>(phiName_);
 
-    volVectorField* momentPtr
+    const dimensionSet momDims
     (
-        new volVectorField
+        phi.dimensions()*dimVelocity/dimVolume
+    );
+
+
+    volVectorField* momentPtr = nullptr;
+
+    word momName(scopedName("momentError"));
+
+    if (zoneSubSetPtr_)
+    {
+        const fvMesh& subMesh = zoneSubSetPtr_->subsetter().subMesh();
+
+        // momentErrorMap
+
+        momentPtr = new volVectorField
         (
             IOobject
             (
-                "momentError",
-                time_.timeName(),
-                mesh_,
+                scopedName("momentErrorMap"),
+                subMesh.time().timeName(),
+                subMesh,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
+            subMesh,
+            dimensionedVector(momDims)
+        );
+
+        subMesh.objectRegistry::store(momentPtr);
+
+        momName = scopedName("momentErrorZone");
+    }
+
+    momentPtr = new volVectorField
+    (
+        IOobject
+        (
+            momName,
+            time_.timeName(),
             mesh_,
-            dimensionedVector(phi.dimensions()*dimVelocity/dimVolume, Zero)
-        )
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector(momDims)
     );
 
     mesh_.objectRegistry::store(momentPtr);
 }
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::functionObjects::momentumError::read(const dictionary& dict)
 {
-    fvMeshFunctionObject::read(dict);
-
-    Info<< type() << " " << name() << ":" << nl;
-
-    // Optional field name entries
-    if (dict.readIfPresent<word>("p", pName_))
+    if (fvMeshFunctionObject::read(dict))
     {
-        Info<< "    p: " << pName_ << endl;
-    }
-    if (dict.readIfPresent<word>("U", UName_))
-    {
-        Info<< "    U: " << UName_ << endl;
+        Info<< type() << ' ' << name() << ':' << nl;
+
+        // Optional field name entries
+        if (dict.readIfPresent<word>("p", pName_))
+        {
+            Info<< "    p: " << pName_ << endl;
+        }
+        if (dict.readIfPresent<word>("U", UName_))
+        {
+            Info<< "    U: " << UName_ << endl;
+        }
+
+        if (dict.readIfPresent<word>("phi", phiName_))
+        {
+            Info<< "    phi: " << phiName_ << endl;
+        }
+
+        if (dict.found("cellZones"))
+        {
+            zoneSubSetPtr_.reset(new Detail::zoneSubSet(mesh_, dict));
+        }
+        else
+        {
+            zoneSubSetPtr_.reset(nullptr);
+        }
+
+        return true;
     }
 
-    if (dict.readIfPresent<word>("phi", phiName_))
-    {
-        Info<< "    phi: " << phiName_ << endl;
-    }
-
-    return true;
+    return false;
 }
 
 
 void Foam::functionObjects::momentumError::calcMomentError()
 {
+    const auto& p = lookupObject<volScalarField>(pName_);
+    const auto& U = lookupObject<volVectorField>(UName_);
+    const auto& phi = lookupObject<surfaceScalarField>(phiName_);
 
-    volVectorField& momentErr =
-        lookupObjectRef<volVectorField>("momentError");
+    if (zoneSubSetPtr_)
+    {
+        const fvMeshSubset& subsetter = zoneSubSetPtr_->subsetter();
 
-    const volScalarField& p = lookupObject<volScalarField>(pName_);
-    const volVectorField& U = lookupObject<volVectorField>(UName_);
-    const surfaceScalarField& phi =
-        lookupObject<surfaceScalarField>(phiName_);
+        fvMesh& subMesh = zoneSubSetPtr_->subsetter().subMesh();
 
-    momentErr = divDevRhoReff() + fvc::div(phi, U) + fvc::grad(p);
+        subMesh.fvSchemes::readOpt() = mesh_.fvSchemes::readOpt();
+        subMesh.fvSchemes::read();
 
+        auto& momentErrMap =
+            subMesh.lookupObjectRef<volVectorField>
+            (
+                scopedName("momentErrorMap")
+            );
+
+        tmp<volScalarField> tp = subsetter.interpolate(p, false);
+        tmp<volVectorField> tU = subsetter.interpolate(U, false);
+        tmp<surfaceScalarField> tphi = subsetter.interpolate(phi, false);
+
+         momentErrMap =
+         (
+            divDevRhoReff()
+          + fvc::div(tphi, tU, "div(phi,U)")
+          + fvc::grad(tp, "grad(p)")
+         );
+    }
+    else
+    {
+        auto& momentErr =
+            lookupObjectRef<volVectorField>(scopedName("momentError"));
+
+        momentErr =  fvc::div(phi, U) + fvc::grad(p) + divDevRhoReff();
+    }
 }
 
 
@@ -201,10 +293,35 @@ bool Foam::functionObjects::momentumError::execute()
 
 bool Foam::functionObjects::momentumError::write()
 {
-    const volVectorField& momentErr =
-        lookupObjectRef<volVectorField>("momentError");
+    Log << "    functionObjects::" << type() << " " << name();
 
-    momentErr.write();
+    if (!zoneSubSetPtr_)
+    {
+        Log << " writing field: " << scopedName("momentError") << endl;
+
+        const auto& momentErr =
+            lookupObjectRef<volVectorField>(scopedName("momentError"));
+
+        momentErr.write();
+    }
+    else
+    {
+        Log << " writing field: " << scopedName("momentErrorMap") << endl;
+
+        const fvMeshSubset& subsetter = zoneSubSetPtr_->subsetter();
+        const fvMesh& subMesh = subsetter.subMesh();
+
+        const auto& momentErrMap =
+            subMesh.lookupObject<volVectorField>
+            (
+                scopedName("momentErrorMap")
+            );
+
+        tmp<volVectorField> mapMomErr =
+            zoneSubSetPtr_->mapToZone<vector>(momentErrMap);
+
+        mapMomErr().write();
+    }
 
     return true;
 }
diff --git a/src/functionObjects/field/momentumError/momentumError.H b/src/functionObjects/field/momentumError/momentumError.H
index 63746c52dd..b8558c4a47 100644
--- a/src/functionObjects/field/momentumError/momentumError.H
+++ b/src/functionObjects/field/momentumError/momentumError.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,14 +54,14 @@ Usage
         U           <UName>;
         phi         <phiName>;
 
-        // Optional (inherited) entries
+        // Inherited entries
         ...
     }
     \endverbatim
 
     where the entries mean:
     \table
-      Property     | Description                        | Type | Req'd | Dflt
+      Property     | Description                        | Type | Reqd  | Deflt
       type         | Type name: momentumError           | word |  yes  | -
       libs         | Library name: fieldFunctionObjects | word |  yes  | -
       p            | Name of pressure field             | word |  no   | p
@@ -71,9 +71,14 @@ Usage
 
     The inherited entries are elaborated in:
      - \link functionObject.H \endlink
+     - \link zoneSubSet.H \endlink
 
     Usage by the \c postProcess utility is not available.
 
+Note
+  - Optionally the user can specify \c cellZones to create a sub-mesh for the
+    \c momentumError calculation.
+
 See also
     - Foam::functionObject
     - Foam::functionObjects::fvMeshFunctionObject
@@ -89,6 +94,7 @@ SourceFiles
 
 #include "fvMeshFunctionObject.H"
 #include "volFieldsFwd.H"
+#include "zoneSubSet.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -120,10 +126,13 @@ protected:
             //- Name of flux field
             word phiName_;
 
+            //- Sub-set mesh
+            autoPtr<Detail::zoneSubSet> zoneSubSetPtr_;
+
 
     // Protected Member Functions
 
-         //- Return the effective viscous stress (laminar + turbulent).
+        //- Return the effective viscous stress (laminar + turbulent).
         tmp<volVectorField> divDevRhoReff();
 
 
diff --git a/src/functionObjects/tools/Make/files b/src/functionObjects/tools/Make/files
new file mode 100644
index 0000000000..7733b4ddab
--- /dev/null
+++ b/src/functionObjects/tools/Make/files
@@ -0,0 +1,3 @@
+zoneSubSet/zoneSubSet.C
+
+LIB = $(FOAM_LIBBIN)/libfunctionObjectsTools
diff --git a/src/functionObjects/tools/Make/options b/src/functionObjects/tools/Make/options
new file mode 100644
index 0000000000..213ece0a3d
--- /dev/null
+++ b/src/functionObjects/tools/Make/options
@@ -0,0 +1,13 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -lfileFormats \
+    -lsurfMesh \
+    -lmeshTools \
+    -ldynamicMesh
diff --git a/src/meshTools/sets/cellSources/haloToCell/haloToCell.C b/src/meshTools/sets/cellSources/haloToCell/haloToCell.C
index c5bc2ff9eb..60c2dc436a 100644
--- a/src/meshTools/sets/cellSources/haloToCell/haloToCell.C
+++ b/src/meshTools/sets/cellSources/haloToCell/haloToCell.C
@@ -28,6 +28,7 @@ License
 #include "haloToCell.H"
 #include "polyMesh.H"
 #include "cellSet.H"
+#include "topoBitSet.H"
 #include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -82,9 +83,16 @@ void Foam::haloToCell::combine(topoSet& set, const bool add) const
     // The starting set of cells
     bitSet current(cells.size());
 
-    for (const label celli : set)
+    if (isA<topoBitSet>(set))
     {
-        current.set(celli);
+        current |= refCast<const topoBitSet>(set).addressing();
+    }
+    else
+    {
+        for (const label celli : set)
+        {
+            current.set(celli);
+        }
     }
 
     // The perimeter faces of the cell set
diff --git a/src/phaseSystemModels/reactingEuler/functionObjects/Make/options b/src/phaseSystemModels/reactingEuler/functionObjects/Make/options
index 0d77a3a0b3..b50eca09b5 100644
--- a/src/phaseSystemModels/reactingEuler/functionObjects/Make/options
+++ b/src/phaseSystemModels/reactingEuler/functionObjects/Make/options
@@ -5,9 +5,11 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/functionObjects/field/lnInclude
+    -I$(LIB_SRC)/functionObjects/field/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lfieldFunctionObjects \
-    -lreactingMultiphaseSystem
+    -lreactingMultiphaseSystem \
+    -ldynamicMesh
diff --git a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/Allrun b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/Allrun
index 97c3b0d0c8..315ebd9f1c 100755
--- a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/Allrun
+++ b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/Allrun
@@ -3,19 +3,9 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
 #------------------------------------------------------------------------------
 
-mkdir -p constant/geometry
+./Allrun.pre
 
-cp -f \
-    "$FOAM_TUTORIALS"/resources/geometry/NACA0012.obj.gz \
-    constant/geometry
-
-restore0Dir
-
-runApplication blockMesh
-
-runApplication transformPoints -scale "(1 0 1)"
-
-runApplication extrudeMesh
+runApplication topoSet
 
 runApplication $(getApplication)
 
diff --git a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/Allrun-parallel b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/Allrun-parallel
new file mode 100755
index 0000000000..2e9d377e12
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/Allrun-parallel
@@ -0,0 +1,18 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+./Allrun.pre
+
+runParallel -s decompose \
+    redistributePar -decompose -overwrite -withZero
+
+runParallel topoSet
+
+runParallel $(getApplication)
+
+runParallel -s reconstruct \
+    redistributePar -reconstruct -overwrite
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/Allrun.pre b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/Allrun.pre
new file mode 100755
index 0000000000..102703e969
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/Allrun.pre
@@ -0,0 +1,20 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+mkdir -p constant/geometry
+
+cp -f \
+    "$FOAM_TUTORIALS"/resources/geometry/NACA0012.obj.gz \
+    constant/geometry
+
+restore0Dir
+
+runApplication blockMesh
+
+runApplication transformPoints -scale "(1 0 1)"
+
+runApplication extrudeMesh
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/controlDict b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/controlDict
index 6a8840cf64..cc4e0f4994 100644
--- a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/controlDict
+++ b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/controlDict
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
+|  \\    /   O peration     | Version:  v2106                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -16,7 +16,7 @@ FoamFile
 
 application     rhoSimpleFoam;
 
-startFrom       latestTime;
+startFrom       startTime;
 
 startTime       0;
 
@@ -34,7 +34,7 @@ purgeWrite      0;
 
 writeFormat     ascii;
 
-writePrecision   8;
+writePrecision  8;
 
 writeCompression off;
 
@@ -48,6 +48,33 @@ functions
 {
     #includeFunc MachNo
     #includeFunc solverInfo
+
+    readFields1
+    {
+        type        readFields;
+        libs        (fieldFunctionObjects);
+        fields      (phi);
+        readOnStart true;
+    }
+
+    contErr
+    {
+        type            div;
+        libs            (fieldFunctionObjects);
+        field           phi;
+        executeControl  writeTime;
+        writeControl    writeTime;
+    }
+
+    momErr
+    {
+        type            momentumError;
+        libs            (fieldFunctionObjects);
+        cellZones       (z1);
+        nLayers         2;
+        executeControl  writeTime;
+        writeControl    writeTime;
+    }
 }
 
 
diff --git a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/decomposeParDict b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/decomposeParDict
new file mode 100644
index 0000000000..c191ee4fef
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/decomposeParDict
@@ -0,0 +1,27 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2106                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 8;
+
+method          hierarchical;
+
+coeffs
+{
+    n           (4 1 2);
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/fvSchemes b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/fvSchemes
index 4eb8843990..da2dd4700e 100644
--- a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/fvSchemes
+++ b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/fvSchemes
@@ -27,6 +27,7 @@ gradSchemes
     grad(U)         $limited;
     grad(k)         $limited;
     grad(omega)     $limited;
+    grad(subsetU)   $limited;
 }
 
 divSchemes
diff --git a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/topoSetDict b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/topoSetDict
new file mode 100644
index 0000000000..d7919c2ee8
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/topoSetDict
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2106                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      topoSetDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    {
+        name    c1;
+        type    cellSet;
+        action  new;
+        source  boxToCell;
+        box     (-0.2 -0.1 -0.2) (1.2 0.1 0.2);
+    }
+
+    {
+        name    z1;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        set     c1;     // cellSet
+    }
+);
+
+
+// ************************************************************************* //