From 610127213398c82d3c07ca6607de031d2c5424f7 Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Tue, 22 Jun 2021 08:46:28 -0700
Subject: [PATCH 1/2] ENH: Adding subMesh option to momentumError and div FOs

1) Adding subMesh capabilities to momentumError and div FOs.
	- A subMesh is created from cellZones.
	- The operators (div, etc) are only calculated in the subMesh.

2) Optionally, halo cells can be added to the cellZones.

3) New helper class to handle the subMesh creation and field mapping.
---
 .../Make/options                              |   1 +
 .../regionFunctionObject.H                    |  12 +-
 .../regionFunctionObjectTemplates.C           |  43 +++-
 src/dynamicMesh/Make/files                    |   2 +
 src/dynamicMesh/zoneSubSet/zoneSubSet.C       | 139 +++++++++++
 src/dynamicMesh/zoneSubSet/zoneSubSet.H       | 204 ++++++++++++++++
 .../zoneSubSet/zoneSubSetTemplates.C          |  70 ++++++
 src/functionObjects/field/Make/options        |   2 +
 src/functionObjects/field/div/div.C           |  24 +-
 src/functionObjects/field/div/div.H           |  23 +-
 src/functionObjects/field/div/divTemplates.C  |  43 +++-
 .../field/fieldExpression/fieldExpression.C   |  28 ++-
 .../field/fieldExpression/fieldExpression.H   |  24 +-
 .../field/momentumError/momentumError.C       | 221 +++++++++++++-----
 .../field/momentumError/momentumError.H       |  17 +-
 src/functionObjects/tools/Make/files          |   3 +
 src/functionObjects/tools/Make/options        |  13 ++
 .../functionObjects/Make/options              |   6 +-
 .../rhoSimpleFoam/aerofoilNACA0012/Allrun     |  14 +-
 .../aerofoilNACA0012/Allrun-parallel          |  18 ++
 .../rhoSimpleFoam/aerofoilNACA0012/Allrun.pre |  20 ++
 .../aerofoilNACA0012/system/controlDict       |  33 ++-
 .../aerofoilNACA0012/system/decomposeParDict  |  27 +++
 .../aerofoilNACA0012/system/fvSchemes         |   1 +
 .../aerofoilNACA0012/system/topoSetDict       |  37 +++
 25 files changed, 914 insertions(+), 111 deletions(-)
 create mode 100644 src/dynamicMesh/zoneSubSet/zoneSubSet.C
 create mode 100644 src/dynamicMesh/zoneSubSet/zoneSubSet.H
 create mode 100644 src/dynamicMesh/zoneSubSet/zoneSubSetTemplates.C
 create mode 100644 src/functionObjects/tools/Make/files
 create mode 100644 src/functionObjects/tools/Make/options
 create mode 100755 tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/Allrun-parallel
 create mode 100755 tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/Allrun.pre
 create mode 100644 tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/decomposeParDict
 create mode 100644 tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/topoSetDict

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/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
+    }
+);
+
+
+// ************************************************************************* //

From acd8e533102adae2c4798c37f79dfb067b5df8c4 Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Tue, 22 Jun 2021 11:07:26 -0700
Subject: [PATCH 2/2] BUG: Managing topoBitSet as source in haloToCell. See
 issue #2137

---
 .../sets/cellSources/haloToCell/haloToCell.C         | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

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