From 9dada5f3f204cfb2d2041644b851c14916fd3c1d Mon Sep 17 00:00:00 2001
From: HenningScheufler <henning.scheufler@dlr.de>
Date: Wed, 12 May 2021 14:32:53 +0200
Subject: [PATCH] ENH: geoVoF module has the capability to run AMR with load
 balancing

code style and quality improvements
renamed recon::centre to interfaceCentre.{groupName}
ranmed recon::normal to interfaceNormal.{groupName}
centre and normal field are not written by default
---
 .../fvMesh/zoneDistribute/zoneDistribute.C    |  45 +++----
 .../fvMesh/zoneDistribute/zoneDistribute.H    |  25 ++--
 .../fvMesh/zoneDistribute/zoneDistributeI.H   |  13 +-
 .../zoneStencils/zoneCPCStencil.C             |  29 +++--
 .../zoneStencils/zoneCPCStencil.H             |  18 ++-
 .../zoneStencils/zoneCellStencils.C           |  54 +++-----
 .../zoneStencils/zoneCellStencils.H           |  38 +++---
 src/functionObjects/field/setFlow/setFlow.C   |  60 ++++-----
 .../isoAdvection/isoAdvection.C               | 115 +++++++++++-------
 .../isoAdvection/isoAdvection.H               |  24 ++--
 .../isoAdvection/isoAdvectionTemplates.C      |  60 +++++----
 .../geometricVoF/cellCuts/cutCell/cutCell.C   |  32 ++++-
 .../geometricVoF/cellCuts/cutCell/cutCell.H   |   2 +-
 .../cellCuts/cutCell/cutCellIso.C             |  46 +------
 .../cellCuts/cutCell/cutCellIso.H             |  41 +++++--
 .../cellCuts/cutCell/cutCellPLIC.C            |  46 +------
 .../cellCuts/cutCell/cutCellPLIC.H            |  35 ++++--
 .../geometricVoF/cellCuts/cutFace/cutFace.C   |  17 +--
 .../cellCuts/cutFace/cutFaceAdvect.C          |  28 +----
 .../cellCuts/cutFace/cutFaceAdvect.H          |  23 +++-
 .../cellCuts/cutFace/cutFaceIso.C             |  24 ----
 .../cellCuts/cutFace/cutFaceIso.H             |  20 ++-
 .../cellCuts/cutFace/cutFacePLIC.C            |  24 ----
 .../cellCuts/cutFace/cutFacePLIC.H            |  20 ++-
 .../reconstructedDistanceFunction.C           |   5 +-
 .../reconstructedDistanceFunction.H           |   4 +-
 .../plicSchemes/gradAlpha/gradAlpha.C         |  15 +--
 .../plicSchemes/gradAlpha/gradAlpha.H         |   4 -
 .../plicSchemes/plicRDF/plicRDF.C             |  58 ++++-----
 .../plicSchemes/plicRDF/plicRDF.H             |   3 -
 .../reconstructionSchemes.C                   |  12 +-
 .../reconstructionSchemes.H                   |  24 +++-
 .../interIsoFoam/damBreak/system/fvSolution   |   1 +
 33 files changed, 481 insertions(+), 484 deletions(-)

diff --git a/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.C b/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.C
index 3cafff4ff9..aad80178c0 100644
--- a/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.C
+++ b/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.C
@@ -91,10 +91,11 @@ Foam::zoneDistribute::coupledFacesPatch() const
 
 Foam::zoneDistribute::zoneDistribute(const fvMesh& mesh)
 :
-    MeshObject<fvMesh, Foam::UpdateableMeshObject, zoneDistribute>(mesh),
-    stencil_(mesh),
+    MeshObject<fvMesh, Foam::TopologicalMeshObject, zoneDistribute>(mesh),
     coupledBoundaryPoints_(coupledFacesPatch()().meshPoints()),
-    send_(Pstream::nProcs())
+    send_(Pstream::nProcs()),
+    stencil_(zoneCPCStencil::New(mesh)),
+    gblIdx_(stencil_.globalNumbering())
 {
 }
 
@@ -103,15 +104,11 @@ Foam::zoneDistribute::zoneDistribute(const fvMesh& mesh)
 
 Foam::zoneDistribute& Foam::zoneDistribute::New(const fvMesh& mesh)
 {
-    zoneDistribute* ptr = mesh.thisDb().getObjectPtr<zoneDistribute>
-    (
-        zoneDistribute::typeName
-    );
+    auto* ptr = mesh.thisDb().getObjectPtr<zoneDistribute>("zoneDistribute");
 
     if (!ptr)
     {
         ptr = new zoneDistribute(mesh);
-
         regIOobject::store(ptr);
     }
 
@@ -123,24 +120,22 @@ Foam::zoneDistribute& Foam::zoneDistribute::New(const fvMesh& mesh)
 
 void Foam::zoneDistribute::updateStencil(const boolList& zone)
 {
-    stencil_.updateStencil(zone);
+    zoneCPCStencil::New(mesh_).updateStencil(zone);
 }
 
 
-void Foam::zoneDistribute::setUpCommforZone
-(
-    const boolList& zone,
-    bool updateStencil
-)
+void Foam::zoneDistribute::setUpCommforZone(const boolList& zone,bool updateStencil)
 {
+    zoneCPCStencil& stencil = zoneCPCStencil::New(mesh_);
+
     if (updateStencil)
     {
-        stencil_.updateStencil(zone);
+        stencil.updateStencil(zone);
     }
 
-    const labelHashSet comms = stencil_.needsComm();
+    const labelHashSet comms = stencil.needsComm();
 
-    List<labelHashSet> needed_(Pstream::nProcs());
+    List<labelHashSet> needed(Pstream::nProcs());
 
     if (Pstream::parRun())
     {
@@ -150,11 +145,10 @@ void Foam::zoneDistribute::setUpCommforZone
             {
                 for (const label gblIdx : stencil_[celli])
                 {
-                    if (!globalNumbering().isLocal(gblIdx))
+                    if (!gblIdx_.isLocal(gblIdx))
                     {
-                        const label procID =
-                            globalNumbering().whichProcID(gblIdx);
-                        needed_[procID].insert(gblIdx);
+                        const label procID = gblIdx_.whichProcID (gblIdx);
+                        needed[procID].insert(gblIdx);
                     }
                 }
             }
@@ -170,7 +164,7 @@ void Foam::zoneDistribute::setUpCommforZone
                 // Put data into send buffer
                 UOPstream toDomain(domain, pBufs);
 
-                toDomain << needed_[domain];
+                toDomain << needed[domain];
             }
         }
 
@@ -193,13 +187,6 @@ void Foam::zoneDistribute::setUpCommforZone
 }
 
 
-void Foam::zoneDistribute::updateMesh(const mapPolyMesh& mpm)
-{
-    if (mesh_.topoChanging())
-    {
-        coupledBoundaryPoints_ = coupledFacesPatch()().meshPoints();
-    }
-}
 
 
 // ************************************************************************* //
diff --git a/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.H b/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.H
index dc432c7745..4edaf8b62b 100644
--- a/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.H
+++ b/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.H
@@ -64,12 +64,10 @@ namespace Foam
 
 class zoneDistribute
 :
-    public MeshObject<fvMesh, UpdateableMeshObject, zoneDistribute>
+    public MeshObject<fvMesh, TopologicalMeshObject, zoneDistribute>
 {
     // Private Data
 
-        //- cell-point-cell stencil elements are in global addressing
-        zoneCPCStencil stencil_;
 
         //- labels of the points on coupled patches
         labelList coupledBoundaryPoints_;
@@ -80,6 +78,10 @@ class zoneDistribute
         //- Return patch of all coupled faces.
         autoPtr<indirectPrimitivePatch> coupledFacesPatch() const;
 
+        zoneCPCStencil& stencil_;
+
+        const globalIndex& gblIdx_;
+
 
         //- Gives patchNumber and patchFaceNumber for a given
         //- Geometric volume field
@@ -118,6 +120,11 @@ public:
         static zoneDistribute& New(const fvMesh&);
 
 
+    //- Destructor
+
+        virtual ~zoneDistribute() = default;
+
+
     // Member Functions
 
         //- Update stencil with boolList the size has to match mesh nCells
@@ -127,15 +134,15 @@ public:
         void updateStencil(const boolList& zone);
 
         //- Stencil reference
-        const labelListList& getStencil()
+        const labelListList& getStencil() noexcept
         {
             return stencil_;
         }
 
         //- Addressing reference
-        const globalIndex& globalNumbering() const
+        const globalIndex& globalNumbering() const noexcept
         {
-            return stencil_.globalNumbering();
+            return gblIdx_;
         }
 
         //- Gives patchNumber and patchFaceNumber for a given
@@ -164,13 +171,7 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& phi
         );
 
-        virtual void updateMesh(const mapPolyMesh& mpm);
 
-        virtual bool movePoints()
-        {
-            // do nothing
-            return false;
-        }
 };
 
 
diff --git a/src/finiteVolume/fvMesh/zoneDistribute/zoneDistributeI.H b/src/finiteVolume/fvMesh/zoneDistribute/zoneDistributeI.H
index 2c733ec6bb..4dd0309075 100644
--- a/src/finiteVolume/fvMesh/zoneDistribute/zoneDistributeI.H
+++ b/src/finiteVolume/fvMesh/zoneDistribute/zoneDistributeI.H
@@ -43,11 +43,8 @@ Type Foam::zoneDistribute::getLocalValue
     {
         return phi[localIdx];
     }
-    else
-    {
-        return faceValue(phi,localIdx);
-    }
 
+    return faceValue(phi,localIdx);
 }
 
 
@@ -88,9 +85,9 @@ Type Foam::zoneDistribute::getValue
     const label gblIdx
 ) const
 {
-    if (globalNumbering().isLocal(gblIdx))
+    if (gblIdx_.isLocal(gblIdx))
     {
-        const label idx = globalNumbering().toLocal(gblIdx);
+        const label idx = gblIdx_.toLocal(gblIdx);
         return getLocalValue(phi,idx);
     }
     else // from other proc
@@ -115,7 +112,6 @@ Foam::Map<Foam::Field<Type>> Foam::zoneDistribute::getFields
             << "do not match. Did the mesh change?"
             << exit(FatalError);
 
-        return Map<Field<Type>>();
     }
 
 
@@ -160,7 +156,6 @@ Foam::Map<Type> Foam::zoneDistribute::getDatafromOtherProc
             << "do not match. Did the mesh change?"
             << exit(FatalError);
 
-        return Map<Type>();
     }
 
 
@@ -177,7 +172,7 @@ Foam::Map<Type> Foam::zoneDistribute::getDatafromOtherProc
                 sendValues[domaini].insert
                 (
                     sendIdx,
-                    getLocalValue(phi,globalNumbering().toLocal(sendIdx))
+                    getLocalValue(phi,gblIdx_.toLocal(sendIdx))
                 );
             }
         }
diff --git a/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCPCStencil.C b/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCPCStencil.C
index 1c91f5478a..c32751b733 100644
--- a/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCPCStencil.C
+++ b/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCPCStencil.C
@@ -123,6 +123,7 @@ void Foam::zoneCPCStencil::calcPointBoundaryData
 
 Foam::zoneCPCStencil::zoneCPCStencil(const fvMesh& mesh)
 :
+    MeshObject<fvMesh, Foam::TopologicalMeshObject, zoneCPCStencil>(mesh),
     zoneCellStencils(mesh),
     nonEmptyBoundaryPoints_(nonEmptyFacesPatch()().meshPoints()),
     uptodate_(mesh.nCells(), false)
@@ -131,6 +132,19 @@ Foam::zoneCPCStencil::zoneCPCStencil(const fvMesh& mesh)
     validBoundaryFaces(isValidBFace_);
 }
 
+Foam::zoneCPCStencil& Foam::zoneCPCStencil::New(const fvMesh& mesh)
+{
+    auto* ptr = mesh.thisDb().getObjectPtr<zoneCPCStencil>("zoneCPCStencil");
+
+    if (!ptr)
+    {
+        ptr = new zoneCPCStencil(mesh);
+        regIOobject::store(ptr);
+    }
+
+    return *ptr;
+}
+
 
 void Foam::zoneCPCStencil::calculateStencil
 (
@@ -224,19 +238,4 @@ void Foam::zoneCPCStencil::calculateStencil
 }
 
 
-void Foam::zoneCPCStencil::updateMesh(const mapPolyMesh& mpm)
-{
-    if (mesh_.topoChanging())
-    {
-        // resize map and globalIndex
-        zoneCellStencils::updateMesh(mpm);
-
-        nonEmptyBoundaryPoints_ = nonEmptyFacesPatch()().meshPoints();
-        uptodate_.resize(mesh_.nCells());
-        uptodate_ = false;
-        validBoundaryFaces(isValidBFace_);
-    }
-}
-
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCPCStencil.H b/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCPCStencil.H
index 318fa38d97..23fcb624c4 100644
--- a/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCPCStencil.H
+++ b/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCPCStencil.H
@@ -58,6 +58,12 @@ namespace Foam
 
 class zoneCPCStencil
 :
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        zoneCPCStencil
+    >,
     public zoneCellStencils
 {
     // Private Data
@@ -97,6 +103,11 @@ class zoneCPCStencil
             labelListList& globalCellCells
         );
 
+        //- No copy construct
+        zoneCPCStencil(const zoneCPCStencil&) = delete;
+
+        //- No copy assignment
+        void operator=(const zoneCPCStencil&) = delete;
 
 public:
 
@@ -109,13 +120,10 @@ public:
         //- Construct from all cells and boundary faces
         explicit zoneCPCStencil(const fvMesh&);
 
+    // Selectors
 
-    // Member Functions
+        static zoneCPCStencil& New(const fvMesh&);
 
-        //- Calculates per cell the neighbour data
-        //-  (= cell or boundary in global numbering).
-        //-  First element is always cell itself!
-        virtual void updateMesh(const mapPolyMesh& mpm);
 };
 
 
diff --git a/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCellStencils.C b/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCellStencils.C
index a397e036ab..9a20a57cc7 100644
--- a/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCellStencils.C
+++ b/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCellStencils.C
@@ -43,7 +43,7 @@ namespace Foam
 Foam::autoPtr<Foam::indirectPrimitivePatch>
 Foam::zoneCellStencils::nonEmptyFacesPatch() const
 {
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const polyBoundaryMesh& patches = meshRef_.boundaryMesh();
 
     label nNonEmpty = 0;
 
@@ -74,10 +74,10 @@ Foam::zoneCellStencils::nonEmptyFacesPatch() const
     (
         IndirectList<face>
         (
-            mesh_.faces(),
+            meshRef_.faces(),
             nonEmptyFaces
         ),
-        mesh_.points()
+        meshRef_.points()
     );
 }
 
@@ -85,7 +85,7 @@ Foam::zoneCellStencils::nonEmptyFacesPatch() const
 Foam::autoPtr<Foam::indirectPrimitivePatch>
 Foam::zoneCellStencils::allCoupledFacesPatch() const
 {
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const polyBoundaryMesh& patches = meshRef_.boundaryMesh();
 
     label nCoupled = 0;
 
@@ -116,19 +116,19 @@ Foam::zoneCellStencils::allCoupledFacesPatch() const
     (
         IndirectList<face>
         (
-            mesh_.faces(),
+            meshRef_.faces(),
             coupledFaces
         ),
-        mesh_.points()
+        meshRef_.points()
     );
 }
 
 
 void Foam::zoneCellStencils::validBoundaryFaces(boolList& isValidBFace) const
 {
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const polyBoundaryMesh& patches = meshRef_.boundaryMesh();
 
-    isValidBFace.setSize(mesh_.nBoundaryFaces());
+    isValidBFace.setSize(meshRef_.nBoundaryFaces());
 
     isValidBFace = true;
 
@@ -136,7 +136,7 @@ void Foam::zoneCellStencils::validBoundaryFaces(boolList& isValidBFace) const
     {
         if (pp.coupled() || isA<emptyPolyPatch>(pp))
         {
-            label bFacei = pp.start()-mesh_.nInternalFaces();
+            label bFacei = pp.start() - meshRef_.nInternalFaces();
             forAll(pp, i)
             {
                 isValidBFace[bFacei++] = false;
@@ -190,22 +190,20 @@ void Foam::zoneCellStencils::insertFaceCells
     labelHashSet& globals
 ) const
 {
-    const labelList& own = mesh_.faceOwner();
-    const labelList& nei = mesh_.faceNeighbour();
+    const labelList& own = meshRef_.faceOwner();
+    const labelList& nei = meshRef_.faceNeighbour();
 
-    forAll(faceLabels, i)
+    for (const label facei : faceLabels)
     {
-        label facei = faceLabels[i];
-
-        label globalOwn = globalNumbering().toGlobal(own[facei]);
+        const label globalOwn = globalNumbering().toGlobal(own[facei]);
         if (globalOwn != exclude0 && globalOwn != exclude1)
         {
             globals.insert(globalOwn);
         }
 
-        if (mesh_.isInternalFace(facei))
+        if (meshRef_.isInternalFace(facei))
         {
-            label globalNei = globalNumbering().toGlobal(nei[facei]);
+            const label globalNei = globalNumbering().toGlobal(nei[facei]);
             if (globalNei != exclude0 && globalNei != exclude1)
             {
                 globals.insert(globalNei);
@@ -213,13 +211,14 @@ void Foam::zoneCellStencils::insertFaceCells
         }
         else
         {
-            const label bFacei = facei-mesh_.nInternalFaces();
+            const label bFacei = facei - meshRef_.nInternalFaces();
 
             if (isValidBFace[bFacei])
             {
                 label globalI = globalNumbering().toGlobal
                 (
-                    mesh_.nCells() + bFacei
+                    meshRef_.nCells()
+                  + bFacei
                 );
 
                 if (globalI != exclude0 && globalI != exclude1)
@@ -258,25 +257,12 @@ Foam::labelList Foam::zoneCellStencils::calcFaceCells
 
 Foam::zoneCellStencils::zoneCellStencils(const fvMesh& mesh)
 :
-    MeshObject<fvMesh, Foam::UpdateableMeshObject, zoneCellStencils>(mesh),
     labelListList(mesh.nCells()),
+    meshRef_(mesh),
     needComm_(),
-    globalNumbering_(mesh_.nCells() + mesh_.nBoundaryFaces())
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::zoneCellStencils::updateMesh(const mapPolyMesh& mpm)
+    globalNumbering_(meshRef_.nCells() + meshRef_.nBoundaryFaces())
 {
-    if (mesh_.topoChanging())
-    {
-        globalNumbering_ =
-            globalIndex(mesh_.nCells() + mesh_.nBoundaryFaces());
 
-        static_cast<labelListList&>(*this) = labelListList(mesh_.nCells());
-        needComm_.clear();
-    }
 }
 
 
diff --git a/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCellStencils.H b/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCellStencils.H
index 16a1cde11b..fd937329d8 100644
--- a/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCellStencils.H
+++ b/src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCellStencils.H
@@ -58,15 +58,19 @@ namespace Foam
 
 class zoneCellStencils
 :
-    public MeshObject<fvMesh, UpdateableMeshObject, zoneCellStencils>,
     public labelListList
 {
 protected:
 
-    // Protected Members
+    // Protected Data
 
+        //- const reference to fvMesh
+        const fvMesh& meshRef_;
+
+        //- cells requiring processor communciation
         labelHashSet needComm_;
 
+        //- Global numbering for cells and boundary faces
         globalIndex globalNumbering_;
 
 
@@ -114,6 +118,12 @@ protected:
         ) = 0;
 
 
+        //- No copy construct
+        zoneCellStencils(const zoneCellStencils&) = delete;
+
+        //- No copy assignment
+        void operator=(const zoneCellStencils&) = delete;
+
 public:
 
         // Declare name of the class and its debug switch
@@ -125,12 +135,15 @@ public:
         //- Construct from all cells and boundary faces
         explicit zoneCellStencils(const fvMesh&);
 
-
-    // Member Functions
-
         //- Calculates per cell the neighbour data
         //  (= cell or boundary in global numbering).
         //  First element is always cell itself!
+        //- Destructor
+        virtual ~zoneCellStencils() = default;
+
+
+    // Member Functions
+
         void updateStencil
         (
             const boolList& zone
@@ -139,23 +152,20 @@ public:
             calculateStencil(zone,*this);
         }
 
-        const labelHashSet& needsComm()
+        const labelHashSet& needsComm() noexcept
         {
             return needComm_;
         }
 
-        //- Global numbering for cells and boundary faces
-        const globalIndex& globalNumbering() const
+        const polyMesh& mesh() const noexcept
         {
-             return globalNumbering_;
+            return meshRef_;
         }
 
-        virtual void updateMesh(const mapPolyMesh& mpm);
-
-        virtual bool movePoints()
+        //- Global numbering for cells and boundary faces
+        const globalIndex& globalNumbering() const noexcept
         {
-            // Do nothing
-            return false;
+             return globalNumbering_;
         }
 };
 
diff --git a/src/functionObjects/field/setFlow/setFlow.C b/src/functionObjects/field/setFlow/setFlow.C
index ede4083de6..fbd45ac896 100644
--- a/src/functionObjects/field/setFlow/setFlow.C
+++ b/src/functionObjects/field/setFlow/setFlow.C
@@ -291,48 +291,42 @@ bool Foam::functionObjects::setFlow::execute()
             Uc.replace(vector::Z, sin(2*pi*x)*sqr(sin(pi*z)));
 
             U = U & R_;
+            U.correctBoundaryConditions();
 
-            // Calculating incompressible flux based on stream function
+            // Calculating phi
             // Note: R_ rotation not implemented in phi calculation
-            const scalarField xp(mesh_.points().component(0) - origin_[0]);
-            const scalarField yp(mesh_.points().component(1) - origin_[1]);
-            const scalarField zp(mesh_.points().component(2) - origin_[2]);
-            const scalarField psi((1.0/pi)*sqr(sin(pi*xp))*sqr(sin(pi*zp)));
+            const vectorField Cf(mesh_.Cf().primitiveField() - origin_);
+            const scalarField Xf(Cf.component(vector::X));
+            const scalarField Yf(Cf.component(vector::Y));
+            const scalarField Zf(Cf.component(vector::Z));
+            vectorField Uf(Xf.size());
+            Uf.replace(vector::X, -sin(2*pi*Zf)*sqr(sin(pi*Xf)));
+            Uf.replace(vector::Y, scalar(0));
+            Uf.replace(vector::Z, sin(2*pi*Xf)*sqr(sin(pi*Zf)));
 
             scalarField& phic = phi.primitiveFieldRef();
-            forAll(phic, fi)
-            {
-                phic[fi] = 0;
-                const face& f = mesh_.faces()[fi];
-                const label nPoints = f.size();
-
-                forAll(f, fpi)
-                {
-                    const label p1 = f[fpi];
-                    const label p2 = f[(fpi + 1) % nPoints];
-                    phic[fi] += 0.5*(psi[p1] + psi[p2])*(yp[p2] - yp[p1]);
-                }
-            }
+            const vectorField& Sfc = mesh_.Sf().primitiveField();
+            phic = Uf & Sfc;
 
             surfaceScalarField::Boundary& phibf = phi.boundaryFieldRef();
+            const surfaceVectorField::Boundary& Sfbf =
+                mesh_.Sf().boundaryField();
+            const surfaceVectorField::Boundary& Cfbf =
+                mesh_.Cf().boundaryField();
+
             forAll(phibf, patchi)
             {
                 scalarField& phif = phibf[patchi];
-                const label start = mesh_.boundaryMesh()[patchi].start();
-
-                forAll(phif, fi)
-                {
-                    phif[fi] = 0;
-                    const face& f = mesh_.faces()[start + fi];
-                    const label nPoints = f.size();
-
-                    forAll(f, fpi)
-                    {
-                        const label p1 = f[fpi];
-                        const label p2 = f[(fpi + 1) % nPoints];
-                        phif[fi] += 0.5*(psi[p1] + psi[p2])*(yp[p2] - yp[p1]);
-                    }
-                }
+                const vectorField& Sff = Sfbf[patchi];
+                const vectorField& Cff = Cfbf[patchi];
+                const scalarField xf(Cff.component(vector::X));
+                const scalarField yf(Cff.component(vector::Y));
+                const scalarField zf(Cff.component(vector::Z));
+                vectorField Ufb(xf.size());
+                Ufb.replace(vector::X, -sin(2*pi*zf)*sqr(sin(pi*xf)));
+                Ufb.replace(vector::Y, scalar(0));
+                Ufb.replace(vector::Z, sin(2*pi*xf)*sqr(sin(pi*zf)));
+                phif = Ufb & Sff;
             }
 
             break;
diff --git a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C
index e64d80dbee..13d80e8250 100644
--- a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C
+++ b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C
@@ -38,6 +38,9 @@ License
 #include "cellSet.H"
 #include "meshTools.H"
 #include "OBJstream.H"
+#include "syncTools.H"
+
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -88,6 +91,7 @@ Foam::isoAdvection::isoAdvection
         dimensionedScalar(dimVol/dimTime, Zero)
     ),
     advectionTime_(0),
+    timeIndex_(-1),
 
     // Tolerances and solution controls
     nAlphaBounds_(dict_.getOrDefault<label>("nAlphaBounds", 3)),
@@ -126,28 +130,69 @@ Foam::isoAdvection::isoAdvection
         mesh_.cells();
 
         // Get boundary mesh and resize the list for parallel comms
-        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-        surfaceCellFacesOnProcPatches_.resize(patches.size());
-
-        // Append all processor patch labels to the list
-        forAll(patches, patchi)
-        {
-            if
-            (
-                isA<processorPolyPatch>(patches[patchi])
-             && patches[patchi].size() > 0
-            )
-            {
-                procPatchLabels_.append(patchi);
-            }
-        }
+        setProcessorPatches();
     }
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::isoAdvection::setProcessorPatches()
+{
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    surfaceCellFacesOnProcPatches_.clear();
+    surfaceCellFacesOnProcPatches_.resize(patches.size());
+
+    // Append all processor patch labels to the list
+    procPatchLabels_.clear();
+    forAll(patches, patchi)
+    {
+        if
+        (
+            isA<processorPolyPatch>(patches[patchi])
+         && !patches[patchi].empty()
+        )
+        {
+            procPatchLabels_.append(patchi);
+        }
+    }
+}
+
+
+void Foam::isoAdvection::extendMarkedCells
+(
+    bitSet& markedCell
+) const
+{
+    // Mark faces using any marked cell
+    bitSet markedFace(mesh_.nFaces());
+
+    for (const label celli : markedCell)
+    {
+        markedFace.set(mesh_.cells()[celli]);  // set multiple faces
+    }
+
+    syncTools::syncFaceList(mesh_, markedFace, orEqOp<unsigned int>());
+
+    // Update cells using any markedFace
+    for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
+    {
+        if (markedFace.test(facei))
+        {
+            markedCell.set(mesh_.faceOwner()[facei]);
+            markedCell.set(mesh_.faceNeighbour()[facei]);
+        }
+    }
+    for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); ++facei)
+    {
+        if (markedFace.test(facei))
+        {
+            markedCell.set(mesh_.faceOwner()[facei]);
+        }
+    }
+}
+
+
 void Foam::isoAdvection::timeIntegratedFlux()
 {
     // Get time step
@@ -196,8 +241,7 @@ void Foam::isoAdvection::timeIntegratedFlux()
 
             // Cell is cut
             const point x0 = surf_->centre()[celli];
-            vector n0 = -surf_->normal()[celli];
-            n0 /= (mag(n0));
+            const vector n0(normalised(-surf_->normal()[celli]));
 
             // Get the speed of the isoface by interpolating velocity and
             // dotting it with isoface unit normal
@@ -212,10 +256,8 @@ void Foam::isoAdvection::timeIntegratedFlux()
             // Note: looping over all cell faces - in reduced-D, some of
             //       these faces will be on empty patches
             const cell& celliFaces = cellFaces[celli];
-            forAll(celliFaces, fi)
+            for (const label facei : celliFaces)
             {
-                const label facei = celliFaces[fi];
-
                 if (mesh_.isInternalFace(facei))
                 {
                     bool isDownwindFace = false;
@@ -279,7 +321,7 @@ void Foam::isoAdvection::timeIntegratedFlux()
         const label patchi = boundaryMesh.patchID()[facei - nInternalFaces];
         const label start = boundaryMesh[patchi].start();
 
-        if (phib[patchi].size())
+        if (!phib[patchi].empty())
         {
             const label patchFacei = facei - start;
             const scalar phiP = phib[patchi][patchFacei];
@@ -332,10 +374,9 @@ void Foam::isoAdvection::setDownwindFaces
     downwindFaces.clear();
 
     // Check all faces of the cell
-    forAll(c, fi)
+    for (const label facei: c)
     {
         // Get face and corresponding flux
-        const label facei = c[fi];
         const scalar phi = faceValue(phi_, facei);
 
         if (own[facei] == celli)
@@ -369,9 +410,8 @@ Foam::scalar Foam::isoAdvection::netFlux
     // Get mesh data
     const labelList& own = mesh_.faceOwner();
 
-    forAll(c, fi)
+    for (const label facei : c)
     {
-        const label facei = c[fi];
         const scalar dVff = faceValue(dVf, facei);
 
         if (own[facei] == celli)
@@ -403,10 +443,8 @@ Foam::DynamicList<Foam::label>  Foam::isoAdvection::syncProcPatches
         PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
 
         // Send
-        forAll(procPatchLabels_, i)
+        for (const label patchi : procPatchLabels_)
         {
-            const label patchi = procPatchLabels_[i];
-
             const processorPolyPatch& procPatch =
                 refCast<const processorPolyPatch>(patches[patchi]);
 
@@ -429,10 +467,8 @@ Foam::DynamicList<Foam::label>  Foam::isoAdvection::syncProcPatches
 
 
         // Receive and combine
-        forAll(procPatchLabels_, patchLabeli)
+        for (const label patchi : procPatchLabels_)
         {
-            const label patchi = procPatchLabels_[patchLabeli];
-
             const processorPolyPatch& procPatch =
                 refCast<const processorPolyPatch>(patches[patchi]);
 
@@ -466,7 +502,7 @@ Foam::DynamicList<Foam::label>  Foam::isoAdvection::syncProcPatches
             {
                 const label facei = faceIDs[i];
                 localFlux[facei] = - nbrdVfs[i];
-                if (debug && mag(localFlux[facei] + nbrdVfs[i]) > 10*SMALL)
+                if (debug && mag(localFlux[facei] + nbrdVfs[i]) > ROOTVSMALL)
                 {
                     Pout<< "localFlux[facei] = " << localFlux[facei]
                         << " and nbrdVfs[i] = " << nbrdVfs[i]
@@ -505,7 +541,7 @@ void Foam::isoAdvection::checkIfOnProcPatch(const label facei)
         const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
         const label patchi = pbm.patchID()[facei - mesh_.nInternalFaces()];
 
-        if (isA<processorPolyPatch>(pbm[patchi]) && pbm[patchi].size())
+        if (isA<processorPolyPatch>(pbm[patchi]) && !pbm[patchi].empty())
         {
             const label patchFacei = pbm[patchi].whichFace(facei);
             surfaceCellFacesOnProcPatches_[patchi].append(patchFacei);
@@ -514,8 +550,6 @@ void Foam::isoAdvection::checkIfOnProcPatch(const label facei)
 }
 
 
-
-
 void Foam::isoAdvection::applyBruteForceBounding()
 {
     bool alpha1Changed = false;
@@ -568,6 +602,7 @@ void Foam::isoAdvection::writeSurfaceCells() const
     }
 }
 
+
 void Foam::isoAdvection::writeIsoFaces
 (
     const DynamicList<List<point>>& faces
@@ -603,10 +638,8 @@ void Foam::isoAdvection::writeIsoFaces
                 const DynamicList<List<point>>& procFacePts =
                     allProcFaces[proci];
 
-                forAll(procFacePts, i)
+                for (const List<point>& facePts : procFacePts)
                 {
-                    const List<point>& facePts = procFacePts[i];
-
                     if (facePts.size() != f.size())
                     {
                         f = face(identity(facePts.size()));
@@ -625,10 +658,8 @@ void Foam::isoAdvection::writeIsoFaces
             << os.name() << nl << endl;
 
         face f;
-        forAll(faces, i)
+        for (const List<point>& facePts : faces)
         {
-            const List<point>& facePts = faces[i];
-
             if (facePts.size() != f.size())
             {
                 f = face(identity(facePts.size()));
diff --git a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.H b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.H
index 41743a30a2..20a415a8f3 100644
--- a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.H
+++ b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.H
@@ -115,6 +115,8 @@ class isoAdvection
         //- Time spent performing interface advection
         scalar advectionTime_;
 
+        //- store timeIndex
+        label timeIndex_;
 
         // Switches and tolerances. Tolerances need to go into toleranceSwitches
 
@@ -179,6 +181,11 @@ class isoAdvection
 
         // Advection functions
 
+            //- Extend markedCell with cell-face-cell.
+            void extendMarkedCells(bitSet& markedCell) const;
+
+            void setProcessorPatches();
+
             //- For each face calculate volumetric face transport during dt
             void timeIntegratedFlux();
 
@@ -202,7 +209,7 @@ class isoAdvection
             template < class SpType, class SuType >
             void boundFlux
             (
-                const scalarField& alpha1,
+                const bitSet& nextToInterface,
                 surfaceScalarField& dVfcorrectionValues,
                 DynamicLabelList& correctedFaces,
                 const SpType& Sp,
@@ -274,6 +281,8 @@ class isoAdvection
             //  list of surface cell faces on processor patches
             void checkIfOnProcPatch(const label facei);
 
+            //- Apply the bounding based on user inputs
+            void applyBruteForceBounding();
 
 public:
 
@@ -303,25 +312,22 @@ public:
         template < class SpType, class SuType >
         void advect(const SpType& Sp, const SuType& Su);
 
-        //- Apply the bounding based on user inputs
-        void applyBruteForceBounding();
-
         // Access functions
 
             //- Return reconstructionSchemes
-            reconstructionSchemes& surf()
+            reconstructionSchemes& surf() noexcept
             {
                 return surf_();
             }
 
             //- Return alpha field
-            const volScalarField& alpha() const
+            const volScalarField& alpha() const noexcept
             {
                 return alpha1_;
             }
 
             //- Return the controls dictionary
-            const dictionary& dict() const
+            const dictionary& dict() const noexcept
             {
                 return dict_;
             }
@@ -365,13 +371,13 @@ public:
             }
 
             //- reference to alphaPhi
-            const surfaceScalarField& alphaPhi() const
+            const surfaceScalarField& alphaPhi() const noexcept
             {
                 return alphaPhi_;
             }
 
             //- time spend in the advection step
-            scalar advectionTime() const
+            scalar advectionTime() const noexcept
             {
                 return advectionTime_;
             }
diff --git a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C
index 0ca1f1b382..f45426186e 100644
--- a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C
+++ b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C
@@ -120,7 +120,7 @@ void Foam::isoAdvection::limitFluxes
 {
     DebugInFunction << endl;
 
-    const scalar aTol = 1.0e-12;          // Note: tolerances
+    const scalar aTol = 100*SMALL;          // Note: tolerances
     scalar maxAlphaMinus1 = gMax(alpha1_) - 1;      // max(alphaNew - 1);
     scalar minAlpha = gMin(alpha1_);           // min(alphaNew);
     const label nOvershoots = 20;         // sum(pos0(alphaNew - 1 - aTol));
@@ -134,6 +134,10 @@ void Foam::isoAdvection::limitFluxes
 
     surfaceScalarField dVfcorrectionValues("dVfcorrectionValues", dVf_*0.0);
 
+    bitSet needBounding(mesh_.nCells(),false);
+    needBounding.set(surfCells_);
+
+    extendMarkedCells(needBounding);
 
     // Loop number of bounding steps
     for (label n = 0; n < nAlphaBounds_; n++)
@@ -143,8 +147,8 @@ void Foam::isoAdvection::limitFluxes
             DebugInfo << "boundAlpha... " << endl;
 
             DynamicList<label> correctedFaces(3*nOvershoots);
-            dVfcorrectionValues = dimensionedScalar(dimVolume, Zero);
-            boundFlux(alpha1In_, dVfcorrectionValues, correctedFaces,Sp,Su);
+            dVfcorrectionValues = dimensionedScalar("0",dimVolume,0.0);
+            boundFlux(needBounding, dVfcorrectionValues, correctedFaces,Sp,Su);
 
             correctedFaces.append
             (
@@ -152,25 +156,24 @@ void Foam::isoAdvection::limitFluxes
             );
 
             labelHashSet alreadyUpdated;
-            forAll(correctedFaces, fi)
+            for (const label facei : correctedFaces)
             {
-                label facei = correctedFaces[fi];
                 if (alreadyUpdated.insert(facei))
                 {
                     checkIfOnProcPatch(facei);
                     const label own = owner[facei];
 
-                    alpha1_[own] +=
-                        -faceValue(dVfcorrectionValues, facei)/mesh_.V()[own];
+                    alpha1_[own] -=
+                        faceValue(dVfcorrectionValues, facei)/mesh_.V()[own];
                     if (mesh_.isInternalFace(facei))
                     {
                         const label nei = neighbour[facei];
-                        alpha1_[nei] -=
-                            -faceValue(dVfcorrectionValues, facei)/mesh_.V()[nei];
+                        alpha1_[nei] +=
+                            faceValue(dVfcorrectionValues, facei)/mesh_.V()[nei];
                     }
 
                     // Change to treat boundaries consistently
-                    scalar corrVf =
+                    const scalar corrVf =
                         faceValue(dVf_, facei)
                       + faceValue(dVfcorrectionValues, facei);
 
@@ -212,7 +215,7 @@ void Foam::isoAdvection::limitFluxes
 template<class SpType, class SuType>
 void Foam::isoAdvection::boundFlux
 (
-    const scalarField& alpha1,
+    const bitSet& nextToInterface,
     surfaceScalarField& dVfcorrectionValues,
     DynamicList<label>& correctedFaces,
     const SpType& Sp,
@@ -223,7 +226,7 @@ void Foam::isoAdvection::boundFlux
     scalar rDeltaT = 1/mesh_.time().deltaTValue();
 
     correctedFaces.clear();
-    scalar aTol = 10*SMALL; // Note: tolerances
+    const scalar aTol = 100*SMALL; // Note: tolerances
 
     const scalarField& meshV = mesh_.cellVolumes();
     const scalar dt = mesh_.time().deltaTValue();
@@ -235,13 +238,15 @@ void Foam::isoAdvection::boundFlux
     const volScalarField& alphaOld = alpha1_.oldTime();
 
     // Loop through alpha cell centred field
-    forAll(alpha1, celli)
+    for (label celli: nextToInterface)
     {
-        if (alpha1[celli] < -aTol || alpha1[celli] > 1 + aTol)
+        if (alpha1_[celli] < -aTol || alpha1_[celli] > 1 + aTol)
         {
             const scalar Vi = meshV[celli];
-            scalar alphaOvershoot = pos0(alpha1[celli]-1)*(alpha1[celli]-1)
-                + neg0(alpha1[celli])*alpha1[celli];
+            scalar alphaOvershoot =
+                pos0(alpha1_[celli] - 1)*(alpha1_[celli] - 1)
+              + neg0(alpha1_[celli])*alpha1_[celli];
+
             scalar fluidToPassOn = alphaOvershoot*Vi;
             label nFacesToPassFluidThrough = 1;
 
@@ -251,7 +256,7 @@ void Foam::isoAdvection::boundFlux
             // not filled and to which dVf < phi*dt
             for (label iter=0; iter<10; iter++)
             {
-                if(mag(alphaOvershoot) < aTol || nFacesToPassFluidThrough == 0)
+                if (mag(alphaOvershoot) < aTol || nFacesToPassFluidThrough == 0)
                 {
                     break;
                 }
@@ -271,9 +276,8 @@ void Foam::isoAdvection::boundFlux
                 scalar dVftot = 0;
                 nFacesToPassFluidThrough = 0;
 
-                forAll(downwindFaces, fi)
+                for (const label facei : downwindFaces)
                 {
-                    const label facei = downwindFaces[fi];
                     const scalar phif = faceValue(phi_, facei);
 
                     const scalar dVff =
@@ -343,7 +347,7 @@ void Foam::isoAdvection::boundFlux
 
                 scalar alpha1New =
                 (
-                    alphaOld[celli]*rDeltaT  + Su[celli]
+                    alphaOld[celli]*rDeltaT + Su[celli]
                   - netFlux(dVf_, celli)/Vi*rDeltaT
                   - netFlux(dVfcorrectionValues, celli)/Vi*rDeltaT
                 )
@@ -351,7 +355,7 @@ void Foam::isoAdvection::boundFlux
                 (rDeltaT - Sp[celli]);
 
                 alphaOvershoot =
-                    pos0(alpha1New-1)*(alpha1New-1)
+                    pos0(alpha1New - 1)*(alpha1New - 1)
                   + neg0(alpha1New)*alpha1New;
 
                 fluidToPassOn = alphaOvershoot*Vi;
@@ -372,13 +376,25 @@ void Foam::isoAdvection::advect(const SpType& Sp, const SuType& Su)
 {
     DebugInFunction << endl;
 
+    if (mesh_.topoChanging())
+    {
+        setProcessorPatches();
+    }
+
     scalar advectionStartTime = mesh_.time().elapsedCpuTime();
 
-    scalar rDeltaT = 1/mesh_.time().deltaTValue();
+    const scalar rDeltaT = 1/mesh_.time().deltaTValue();
 
     // reconstruct the interface
     surf_->reconstruct();
 
+    if (timeIndex_ < mesh_.time().timeIndex())
+    {
+        timeIndex_= mesh_.time().timeIndex();
+        surf_->normal().oldTime() = surf_->normal();
+        surf_->centre().oldTime() = surf_->centre();
+    }
+
     // Initialising dVf with upwind values
     // i.e. phi[facei]*alpha1[upwindCell[facei]]*dt
     dVf_ = upwind<scalar>(mesh_, phi_).flux(alpha1_)*mesh_.time().deltaT();
diff --git a/src/transportModels/geometricVoF/cellCuts/cutCell/cutCell.C b/src/transportModels/geometricVoF/cellCuts/cutCell/cutCell.C
index 05dbc8677f..080b801f22 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutCell/cutCell.C
+++ b/src/transportModels/geometricVoF/cellCuts/cutCell/cutCell.C
@@ -37,8 +37,17 @@ int Foam::cutCell::debug = 0;
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::cutCell::cutCell(const fvMesh&)
-{}
+Foam::cutCell::cutCell
+(
+    const fvMesh& mesh
+)
+{
+    // required as otherwise setAlphaFields might not work in parallel
+    mesh.C();
+    mesh.V();
+    mesh.Cf();
+    mesh.magSf();
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -47,7 +56,8 @@ void Foam::cutCell::calcCellData
 (
     const DynamicList<point>& cutFaceCentres,
     const DynamicList<vector>& cutFaceAreas,
-    vector& subCellCentre, scalar& subCellVolume
+    vector& subCellCentre,
+    scalar& subCellVolume
 )
 {
     // Clear the fields for accumulation
@@ -159,11 +169,27 @@ void Foam::cutCell::calcIsoFacePointsFromEdges
     DynamicList<point>& facePoints
 )
 {
+    if (mag(faceArea) < VSMALL)
+    {
+        facePoints.clear();
+        return;
+    }
     const vector zhat = normalised(faceArea);
     vector xhat = faceEdges[0][0] - faceCentre;
     xhat = (xhat - (xhat & zhat)*zhat);
     xhat.normalise();
+    if (mag(xhat) == 0)
+    {
+        facePoints.clear();
+        return;
+    }
     vector yhat = normalised(zhat ^ xhat);
+    if (mag(yhat) == 0)
+    {
+        facePoints.clear();
+        return;
+    }
+    yhat.normalise();
 
     // Calculating all intersection points
     DynamicList<point> unsortedFacePoints(3 * faceEdges.size());
diff --git a/src/transportModels/geometricVoF/cellCuts/cutCell/cutCell.H b/src/transportModels/geometricVoF/cellCuts/cutCell/cutCell.H
index 7ae86d81a5..29f183c1fe 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutCell/cutCell.H
+++ b/src/transportModels/geometricVoF/cellCuts/cutCell/cutCell.H
@@ -99,7 +99,7 @@ public:
     // Constructors
 
         //- Construct from fvMesh
-        explicit cutCell(const fvMesh& unused);
+        explicit cutCell(const fvMesh& mesh);
 };
 
 
diff --git a/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellIso.C b/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellIso.C
index feb54df002..0ea299829c 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellIso.C
+++ b/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellIso.C
@@ -115,7 +115,7 @@ Foam::label Foam::cutCellIso::calcSubCell
         // In the rare but occuring cases where a cell is only touched at a
         // point or a line the isoFaceArea_ will have zero length and here the
         // cell should be treated as either completely empty or full.
-        if (mag(faceArea_) < 10*SMALL)
+        if (mag(faceArea_) < ROOTVSMALL)
         {
             if (nFaceBelowInterface == 0)
             {
@@ -170,21 +170,9 @@ Foam::label Foam::cutCellIso::calcSubCell
 }
 
 
-const Foam::point& Foam::cutCellIso::subCellCentre() const
-{
-    return subCellCentre_;
-}
-
-
-Foam::scalar Foam::cutCellIso::subCellVolume() const
-{
-    return subCellVolume_;
-}
-
-
 const Foam::DynamicList<Foam::point>& Foam::cutCellIso::facePoints()
 {
-    if (facePoints_.size() == 0)
+    if (facePoints_.empty())
     {
         // get face points in sorted order
         calcIsoFacePointsFromEdges
@@ -200,36 +188,6 @@ const Foam::DynamicList<Foam::point>& Foam::cutCellIso::facePoints()
 }
 
 
-const Foam::point& Foam::cutCellIso::faceCentre() const
-{
-    return faceCentre_;
-}
-
-
-const Foam::vector& Foam::cutCellIso::faceArea() const
-{
-    return faceArea_;
-}
-
-
-Foam::label Foam::cutCellIso::cellStatus() const
-{
-    return cellStatus_;
-}
-
-
-Foam::scalar Foam::cutCellIso::VolumeOfFluid() const
-{
-    return VOF_;
-}
-
-
-Foam::scalar Foam::cutCellIso::cutValue() const
-{
-    return cutValue_;
-}
-
-
 void Foam::cutCellIso::clearStorage()
 {
     cellI_ = -1;
diff --git a/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellIso.H b/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellIso.H
index 3e905ff70c..bf982952ed 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellIso.H
+++ b/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellIso.H
@@ -148,28 +148,49 @@ class cutCellIso
         label calcSubCell(const label cellI, const scalar cutValue);
 
         //- Returns subCellCentre
-        const point& subCellCentre() const;
+        const point& subCellCentre() const noexcept
+        {
+            return subCellCentre_;
+        }
 
         //- Returns subCellVolume
-        scalar subCellVolume() const;
+        scalar subCellVolume() const noexcept
+        {
+            return subCellVolume_;
+        }
 
-        //- Returns the points of the cutting isoface
+        //- Returns the points of the cutting PLICface
         const DynamicList<point>& facePoints();
 
-        //- Returns the centre of the cutting isoface
-        const point& faceCentre() const;
+        //- Returns the centre of the cutting PLICface
+        const point& faceCentre() const noexcept
+        {
+            return faceCentre_;
+        }
 
-        //- Returns the area normal vector of the cutting isoface
-        const vector& faceArea() const;
+        //- Returns the area normal vector of the cutting PLICface
+        const vector& faceArea() const noexcept
+        {
+            return faceArea_;
+        }
 
         //- Returns cellStatus
-        label cellStatus() const;
+        label cellStatus() const noexcept
+        {
+            return cellStatus_;
+        }
 
         //- Returns volume of fluid value
-        scalar VolumeOfFluid() const;
+        scalar VolumeOfFluid() const noexcept
+        {
+            return VOF_;
+        }
 
         //- Returns cutValue
-        scalar cutValue() const;
+        scalar cutValue() const noexcept
+        {
+            return cutValue_;
+        }
 
         //- Resets internal values
         void clearStorage();
diff --git a/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellPLIC.C b/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellPLIC.C
index ca5fcd9f8f..7a014ddde9 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellPLIC.C
+++ b/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellPLIC.C
@@ -117,7 +117,7 @@ Foam::label Foam::cutCellPLIC::calcSubCell
         // In the rare but occuring cases where a cell is only touched at a
         // point or a line the isoFaceArea_ will have zero length and here the
         // cell should be treated as either completely empty or full.
-        if (mag(faceArea_) < 10*SMALL)
+        if (mag(faceArea_) < ROOTVSMALL)
         {
             if (nFaceBelowInterface == 0)
             {
@@ -172,21 +172,9 @@ Foam::label Foam::cutCellPLIC::calcSubCell
 }
 
 
-const Foam::point& Foam::cutCellPLIC::subCellCentre() const
-{
-    return subCellCentre_;
-}
-
-
-Foam::scalar Foam::cutCellPLIC::subCellVolume() const
-{
-    return subCellVolume_;
-}
-
-
 const Foam::DynamicList<Foam::point>& Foam::cutCellPLIC::facePoints()
 {
-    if (facePoints_.size() == 0)
+    if (facePoints_.empty())
     {
         // get face points in sorted order
         calcIsoFacePointsFromEdges
@@ -202,36 +190,6 @@ const Foam::DynamicList<Foam::point>& Foam::cutCellPLIC::facePoints()
 }
 
 
-const Foam::point& Foam::cutCellPLIC::faceCentre() const
-{
-    return faceCentre_;
-}
-
-
-const Foam::vector& Foam::cutCellPLIC::faceArea() const
-{
-    return faceArea_;
-}
-
-
-Foam::scalar Foam::cutCellPLIC::VolumeOfFluid() const
-{
-    return VOF_;
-}
-
-
-Foam::label Foam::cutCellPLIC::cellStatus() const
-{
-    return cellStatus_;
-}
-
-
-Foam::scalar Foam::cutCellPLIC::cutValue() const
-{
-    return cutValue_;
-}
-
-
 void Foam::cutCellPLIC::clearStorage()
 {
     cellI_ = -1;
diff --git a/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellPLIC.H b/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellPLIC.H
index c5aadab2f9..4f00032390 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellPLIC.H
+++ b/src/transportModels/geometricVoF/cellCuts/cutCell/cutCellPLIC.H
@@ -147,28 +147,49 @@ class cutCellPLIC
         );
 
         //- Returns subCellCentre
-        const point& subCellCentre() const;
+        const point& subCellCentre() const noexcept
+        {
+            return subCellCentre_;
+        }
 
         //- Returns subCellVolume
-        scalar subCellVolume() const;
+        scalar subCellVolume() const noexcept
+        {
+            return subCellVolume_;
+        }
 
         //- Returns the points of the cutting PLICface
         const DynamicList<point>& facePoints();
 
         //- Returns the centre of the cutting PLICface
-        const point& faceCentre() const;
+        const point& faceCentre() const noexcept
+        {
+            return faceCentre_;
+        }
 
         //- Returns the area normal vector of the cutting PLICface
-        const vector& faceArea() const;
+        const vector& faceArea() const noexcept
+        {
+            return faceArea_;
+        }
 
         //- Returns cellStatus
-        label cellStatus() const;
+        label cellStatus() const noexcept
+        {
+            return cellStatus_;
+        }
 
         //- Returns volume of fluid value
-        scalar VolumeOfFluid() const;
+        scalar VolumeOfFluid() const noexcept
+        {
+            return VOF_;
+        }
 
         //- Returns cutValue
-        scalar cutValue() const;
+        scalar cutValue() const noexcept
+        {
+            return cutValue_;
+        }
 
         //- Resets internal values
         void clearStorage();
diff --git a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFace.C b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFace.C
index 7b86efff17..8f3fce2988 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFace.C
+++ b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFace.C
@@ -83,8 +83,8 @@ void Foam::cutFace::calcSubFace
 
         if
         (
-            (pointStatus[idx] < 0 && pointStatus[nextIdx] > 0) ||
-            (pointStatus[idx] > 0 && pointStatus[nextIdx] < 0)
+            (pointStatus[idx] < 0 && pointStatus[nextIdx] > 0)
+         || (pointStatus[idx] > 0 && pointStatus[nextIdx] < 0)
         ) // cut on edge cut Value is zero
         {
             label nextP = f.nextLabel(idx);
@@ -161,8 +161,8 @@ void Foam::cutFace::calcSubFace
 
         if
         (
-            (pointStatus[idx] < 0 && pointStatus[nextIdx] > 0) ||
-            (pointStatus[idx] > 0 && pointStatus[nextIdx] < 0)
+            (pointStatus[idx] < 0 && pointStatus[nextIdx] > 0)
+         || (pointStatus[idx] > 0 && pointStatus[nextIdx] < 0)
         ) // cut on edge cut Value is zero
         {
             label nextP = f.nextLabel(idx);
@@ -232,8 +232,8 @@ void Foam::cutFace::calcSubFace
 
         if
         (
-            (pointStatus[idx] < 0 && pointStatus[nextIdx] > 0) ||
-            (pointStatus[idx] > 0 &&  pointStatus[nextIdx] < 0)
+            (pointStatus[idx] < 0 && pointStatus[nextIdx] > 0)
+         || (pointStatus[idx] > 0 &&  pointStatus[nextIdx] < 0)
         )
         {
             label nextP = f.nextLabel(idx);
@@ -328,7 +328,10 @@ void Foam::cutFace::calcSubFaceCentreAndArea
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::cutFace::cutFace(const fvMesh& mesh)
+Foam::cutFace::cutFace
+(
+    const fvMesh& mesh
+)
 :
     mesh_(mesh)
 {}
diff --git a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceAdvect.C b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceAdvect.C
index c5bd8b589a..0091977e79 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceAdvect.C
+++ b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceAdvect.C
@@ -246,7 +246,7 @@ Foam::scalar Foam::cutFaceAdvect::timeIntegratedFaceFlux
 
         // Check if pTimes changes direction more than twice when looping face
         label nShifts = 0;
-        forAll(pTimes_, pi) // i have no clue what this is
+        forAll(pTimes_, pi)
         {
             const label oldEdgeSign =
                 sign(pTimes_[(pi + 1) % nPoints] - pTimes_[pi]);
@@ -356,7 +356,7 @@ Foam::scalar Foam::cutFaceAdvect::timeIntegratedFaceFlux
 
         // Check if pTimes changes direction more than twice when looping face
         label nShifts = 0;
-        forAll(pTimes_, pi) // i have no clue what this is
+        forAll(pTimes_, pi)
         {
             const label oldEdgeSign =
                 sign(pTimes_[(pi + 1) % nPoints] - pTimes_[pi]);
@@ -984,30 +984,6 @@ void Foam::cutFaceAdvect::cutPoints
 }
 
 
-const Foam::point& Foam::cutFaceAdvect::subFaceCentre() const
-{
-    return subFaceCentre_;
-}
-
-
-const Foam::vector& Foam::cutFaceAdvect::subFaceArea() const
-{
-    return subFaceArea_;
-}
-
-
-const Foam::DynamicList<Foam::point>& Foam::cutFaceAdvect::subFacePoints() const
-{
-    return subFacePoints_;
-}
-
-
-const Foam::DynamicList<Foam::point>& Foam::cutFaceAdvect::surfacePoints() const
-{
-    return surfacePoints_;
-}
-
-
 void Foam::cutFaceAdvect::clearStorage()
 {
     subFaceCentre_ = Zero;
diff --git a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceAdvect.H b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceAdvect.H
index 51f3538193..7332de506c 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceAdvect.H
+++ b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceAdvect.H
@@ -232,16 +232,29 @@ class cutFaceAdvect
         );
 
         //- Returns centre of cutted face
-        const point& subFaceCentre() const;
+        const point& subFaceCentre() const noexcept
+        {
+            return subFaceCentre_;
+        }
 
         //- Returns area vector of cutted face
-        const vector& subFaceArea() const;
+        const vector& subFaceArea() const noexcept
+        {
+            return subFaceArea_;
+        }
 
         //- Returns the cut edge of the cutted face
-        const DynamicList<point>& subFacePoints() const;
+        const DynamicList<point>& subFacePoints() const noexcept
+        {
+            return subFacePoints_;
+        }
+
+        //- Returns point of the face in sorted of cutted face
+        const DynamicList<point>& surfacePoints() const noexcept
+        {
+            return surfacePoints_;
+        }
 
-        //- Returns point of the cutted face in sorted order
-        const DynamicList<point>& surfacePoints() const;
 
         //- Resets internal variables
         void clearStorage();
diff --git a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceIso.C b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceIso.C
index 778109da11..3db29bc82a 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceIso.C
+++ b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceIso.C
@@ -111,30 +111,6 @@ Foam::label Foam::cutFaceIso::calcSubFace
 }
 
 
-const Foam::point& Foam::cutFaceIso::subFaceCentre() const
-{
-    return subFaceCentre_;
-}
-
-
-const Foam::vector& Foam::cutFaceIso::subFaceArea() const
-{
-    return subFaceArea_;
-}
-
-
-const Foam::DynamicList<Foam::point>& Foam::cutFaceIso::subFacePoints() const
-{
-    return subFacePoints_;
-}
-
-
-const Foam::DynamicList<Foam::point>& Foam::cutFaceIso::surfacePoints() const
-{
-    return surfacePoints_;
-}
-
-
 void Foam::cutFaceIso::clearStorage()
 {
     subFaceCentre_ = Zero;
diff --git a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceIso.H b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceIso.H
index 8051dc80d8..9ce3fcfa5b 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceIso.H
+++ b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceIso.H
@@ -124,16 +124,28 @@ public:
         );
 
         //- Returns centre of cutted face
-        const point& subFaceCentre() const;
+        const point& subFaceCentre() const noexcept
+        {
+            return subFaceCentre_;
+        }
 
         //- Returns area vector of cutted face
-        const vector& subFaceArea() const;
+        const vector& subFaceArea() const noexcept
+        {
+            return subFaceArea_;
+        }
 
         //- Returns the cut edge of the cutted face
-        const DynamicList<point>& subFacePoints() const;
+        const DynamicList<point>& subFacePoints() const noexcept
+        {
+            return subFacePoints_;
+        }
 
         //- Returns point of the face in sorted of cutted face
-        const DynamicList<point>& surfacePoints() const;
+        const DynamicList<point>& surfacePoints() const noexcept
+        {
+            return surfacePoints_;
+        }
 
         //- Resets internal variables
         void clearStorage();
diff --git a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFacePLIC.C b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFacePLIC.C
index 446b026e61..02dcb3ebca 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFacePLIC.C
+++ b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFacePLIC.C
@@ -115,30 +115,6 @@ Foam::label Foam::cutFacePLIC::calcSubFace
 }
 
 
-const Foam::point& Foam::cutFacePLIC::subFaceCentre() const
-{
-    return subFaceCentre_;
-}
-
-
-const Foam::vector& Foam::cutFacePLIC::subFaceArea() const
-{
-    return subFaceArea_;
-}
-
-
-const Foam::DynamicList<Foam::point>& Foam::cutFacePLIC::subFacePoints() const
-{
-    return subFacePoints_;
-}
-
-
-const Foam::DynamicList<Foam::point>& Foam::cutFacePLIC::surfacePoints() const
-{
-    return surfacePoints_;
-}
-
-
 void Foam::cutFacePLIC::clearStorage()
 {
     subFaceCentre_ = Zero;
diff --git a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFacePLIC.H b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFacePLIC.H
index 75b76a44d8..183553cb06 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFacePLIC.H
+++ b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFacePLIC.H
@@ -123,16 +123,28 @@ public:
         );
 
         //- Returns centre of cutted face
-        const point& subFaceCentre() const;
+        const point& subFaceCentre() const noexcept
+        {
+            return subFaceCentre_;
+        }
 
         //- Returns area vector of cutted face
-        const vector& subFaceArea() const;
+        const vector& subFaceArea() const noexcept
+        {
+            return subFaceArea_;
+        }
 
         //- Returns the cut edge of the cutted face
-        const DynamicList<point>& subFacePoints() const;
+        const DynamicList<point>& subFacePoints() const noexcept
+        {
+            return subFacePoints_;
+        }
 
         //- Returns point of the face in sorted of cutted face
-        const DynamicList<point>& surfacePoints() const;
+        const DynamicList<point>& surfacePoints() const noexcept
+        {
+            return surfacePoints_;
+        }
 
         //- Resets internal variables
         void clearStorage();
diff --git a/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C b/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C
index 30777d8bd0..fedd64a909 100644
--- a/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C
+++ b/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C
@@ -104,7 +104,7 @@ void Foam::reconstructedDistanceFunction::markCellsNearSurf
     // do coupled face first
     Map<bool> syncMap;
 
-    for (int level=0;level<=neiRingLevel;level++)
+    for (label level=0;level<=neiRingLevel;level++)
     {
         // parallel
         if (level > 0)
@@ -280,9 +280,8 @@ const Foam::volScalarField&  Foam::reconstructedDistanceFunction::constructRDF
                 scalar avgWeight = 0;
                 const point p = mesh_.C()[celli];
 
-                forAll(stencil[celli], i)
+                for (const label gblIdx : stencil[celli])
                 {
-                    const label gblIdx = stencil[celli][i];
                     vector n = -distribute.getValue(normal, mapNormal, gblIdx);
                     if (mag(n) != 0)
                     {
diff --git a/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.H b/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.H
index 278087e6d1..2939cea3e9 100644
--- a/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.H
+++ b/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.H
@@ -108,12 +108,12 @@ public:
             surfaceVectorField::Boundary& nHatb
         );
 
-        const volScalarField& cellDistLevel() const
+        const volScalarField& cellDistLevel() const noexcept
         {
             return cellDistLevel_;
         }
 
-        const boolList& nextToInterface() const
+        const boolList& nextToInterface() const noexcept
         {
             return nextToInterface_;
         }
diff --git a/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/gradAlpha/gradAlpha.C b/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/gradAlpha/gradAlpha.C
index b65b323946..ea414981c8 100644
--- a/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/gradAlpha/gradAlpha.C
+++ b/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/gradAlpha/gradAlpha.C
@@ -48,21 +48,23 @@ void Foam::reconstruction::gradAlpha::gradSurf(const volScalarField& phi)
 {
     leastSquareGrad<scalar> lsGrad("polyDegree1",mesh_.geometricD());
 
-    exchangeFields_.setUpCommforZone(interfaceCell_,true);
+    zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
+
+    exchangeFields.setUpCommforZone(interfaceCell_,true);
 
     Map<vector> mapCC
     (
-        exchangeFields_.getDatafromOtherProc(interfaceCell_, mesh_.C())
+        exchangeFields.getDatafromOtherProc(interfaceCell_, mesh_.C())
     );
     Map<scalar> mapPhi
     (
-        exchangeFields_.getDatafromOtherProc(interfaceCell_, phi)
+        exchangeFields.getDatafromOtherProc(interfaceCell_, phi)
     );
 
     DynamicField<vector> cellCentre(100);
     DynamicField<scalar> phiValues(100);
 
-    const labelListList& stencil = exchangeFields_.getStencil();
+    const labelListList& stencil = exchangeFields.getStencil();
 
     forAll(interfaceLabels_, i)
     {
@@ -75,11 +77,11 @@ void Foam::reconstruction::gradAlpha::gradSurf(const volScalarField& phi)
         {
             cellCentre.append
             (
-                exchangeFields_.getValue(mesh_.C(), mapCC, gblIdx)
+                exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
             );
             phiValues.append
             (
-                exchangeFields_.getValue(phi, mapPhi, gblIdx)
+                exchangeFields.getValue(phi, mapPhi, gblIdx)
             );
         }
 
@@ -111,7 +113,6 @@ Foam::reconstruction::gradAlpha::gradAlpha
     interfaceNormal_(fvc::grad(alpha1)),
     isoFaceTol_(modelDict().getOrDefault<scalar>("isoFaceTol", 1e-8)),
     surfCellTol_(modelDict().getOrDefault<scalar>("surfCellTol", 1e-8)),
-    exchangeFields_(zoneDistribute::New(mesh_)),
     sIterPLIC_(mesh_,surfCellTol_)
 {
     reconstruct();
diff --git a/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/gradAlpha/gradAlpha.H b/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/gradAlpha/gradAlpha.H
index 4933ddd005..f895ff625c 100644
--- a/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/gradAlpha/gradAlpha.H
+++ b/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/gradAlpha/gradAlpha.H
@@ -83,13 +83,9 @@ class gradAlpha
         //  Those with surfCellTol_ < alpha1 < 1 - surfCellTol_
         scalar surfCellTol_;
 
-        //- Provides stencil and map
-        zoneDistribute& exchangeFields_;
-
         //- SurfaceIterator finds the plane centre for specified VOF value
         surfaceIteratorPLIC sIterPLIC_;
 
-
     // Private Member Functions
 
         //- Compute gradient at the surfaces
diff --git a/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C b/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C
index 03ab5d4285..fd766484c2 100644
--- a/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C
+++ b/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C
@@ -48,27 +48,28 @@ namespace reconstruction
 void Foam::reconstruction::plicRDF::interpolateNormal()
 {
     scalar dt = mesh_.time().deltaTValue();
+    zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
 
     leastSquareGrad<scalar> lsGrad("polyDegree1",mesh_.geometricD());
 
-    exchangeFields_.setUpCommforZone(interfaceCell_,false);
+    exchangeFields.setUpCommforZone(interfaceCell_,false);
 
     Map<vector> mapCentre
     (
-        exchangeFields_.getDatafromOtherProc(interfaceCell_, centre_)
+        exchangeFields.getDatafromOtherProc(interfaceCell_, centre_)
     );
     Map<vector> mapNormal
     (
-        exchangeFields_.getDatafromOtherProc(interfaceCell_, normal_)
+        exchangeFields.getDatafromOtherProc(interfaceCell_, normal_)
     );
 
     Map<vector> mapCC
     (
-        exchangeFields_.getDatafromOtherProc(interfaceCell_, mesh_.C())
+        exchangeFields.getDatafromOtherProc(interfaceCell_, mesh_.C())
     );
     Map<scalar> mapAlpha
     (
-        exchangeFields_.getDatafromOtherProc(interfaceCell_, alpha1_)
+        exchangeFields.getDatafromOtherProc(interfaceCell_, alpha1_)
     );
 
     DynamicField<vector > cellCentre(100);
@@ -76,7 +77,7 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
 
     DynamicList<vector> foundNormals(30);
 
-    const labelListList& stencil = exchangeFields_.getStencil();
+    const labelListList& stencil = exchangeFields.getStencil();
 
     forAll(interfaceLabels_, i)
     {
@@ -84,17 +85,16 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
         vector estimatedNormal{Zero};
         scalar weight{0};
         foundNormals.clear();
-        forAll(stencil[celli], i)
+        for (const label gblIdx : stencil[celli])
         {
-            const label gblIdx = stencil[celli][i];
             vector n =
-                exchangeFields_.getValue(normal_, mapNormal, gblIdx);
+                exchangeFields.getValue(normal_, mapNormal, gblIdx);
             point p = mesh_.C()[celli]-U_[celli]*dt;
             if (mag(n) != 0)
             {
                 n /= mag(n);
                 vector centre =
-                    exchangeFields_.getValue(centre_, mapCentre, gblIdx);
+                    exchangeFields.getValue(centre_, mapCentre, gblIdx);
                 vector distanceToIntSeg = (tensor::I- n*n) & (p - centre);
                 estimatedNormal += n /max(mag(distanceToIntSeg), SMALL);
                 weight += 1/max(mag(distanceToIntSeg), SMALL);
@@ -143,11 +143,11 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
                 const label gblIdx = stencil[celli][i];
                 cellCentre.append
                 (
-                    exchangeFields_.getValue(mesh_.C(), mapCC, gblIdx)
+                    exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
                 );
                 alphaValues.append
                 (
-                    exchangeFields_.getValue(alpha1_, mapAlpha, gblIdx)
+                    exchangeFields.getValue(alpha1_, mapAlpha, gblIdx)
                 );
             }
             cellCentre -= mesh_.C()[celli];
@@ -160,22 +160,23 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
 void Foam::reconstruction::plicRDF::gradSurf(const volScalarField& phi)
 {
     leastSquareGrad<scalar> lsGrad("polyDegree1", mesh_.geometricD());
+    zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
 
-    exchangeFields_.setUpCommforZone(interfaceCell_, false);
+    exchangeFields.setUpCommforZone(interfaceCell_, false);
 
     Map<vector> mapCC
     (
-        exchangeFields_.getDatafromOtherProc(interfaceCell_, mesh_.C())
+        exchangeFields.getDatafromOtherProc(interfaceCell_, mesh_.C())
     );
     Map<scalar> mapPhi
     (
-        exchangeFields_.getDatafromOtherProc(interfaceCell_, phi)
+        exchangeFields.getDatafromOtherProc(interfaceCell_, phi)
     );
 
     DynamicField<vector> cellCentre(100);
     DynamicField<scalar> phiValues(100);
 
-    const labelListList& stencil = exchangeFields_.getStencil();
+    const labelListList& stencil = exchangeFields.getStencil();
 
     forAll(interfaceLabels_, i)
     {
@@ -188,11 +189,11 @@ void Foam::reconstruction::plicRDF::gradSurf(const volScalarField& phi)
         {
             cellCentre.append
             (
-                exchangeFields_.getValue(mesh_.C(), mapCC, gblIdx)
+                exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
             );
             phiValues.append
             (
-                exchangeFields_.getValue(phi, mapPhi, gblIdx)
+                exchangeFields.getValue(phi, mapPhi, gblIdx)
             );
         }
 
@@ -204,6 +205,8 @@ void Foam::reconstruction::plicRDF::gradSurf(const volScalarField& phi)
 
 void Foam::reconstruction::plicRDF::setInitNormals(bool interpolate)
 {
+    zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
+
     interfaceLabels_.clear();
 
     forAll(alpha1_, celli)
@@ -218,7 +221,7 @@ void Foam::reconstruction::plicRDF::setInitNormals(bool interpolate)
 
     RDF_.markCellsNearSurf(interfaceCell_, 1);
     const boolList& nextToInterface_ = RDF_.nextToInterface();
-    exchangeFields_.updateStencil(nextToInterface_);
+    exchangeFields.updateStencil(nextToInterface_);
 
     if (interpolate)
     {
@@ -236,15 +239,15 @@ void Foam::reconstruction::plicRDF::calcResidual
     List<normalRes>& normalResidual
 )
 {
-    exchangeFields_.setUpCommforZone(interfaceCell_,false);
+    zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
+    exchangeFields.setUpCommforZone(interfaceCell_,false);
 
     Map<vector> mapNormal
     (
-        exchangeFields_.getDatafromOtherProc(interfaceCell_, normal_)
+        exchangeFields.getDatafromOtherProc(interfaceCell_, normal_)
     );
 
-    const labelListList& stencil = exchangeFields_.getStencil();
-
+    const labelListList& stencil = exchangeFields.getStencil();
 
     forAll(interfaceLabels_, i)
     {
@@ -265,7 +268,7 @@ void Foam::reconstruction::plicRDF::calcResidual
         {
             const label gblIdx = stencil[celli][j];
             vector normal =
-                exchangeFields_.getValue(normal_, mapNormal, gblIdx);
+                exchangeFields.getValue(normal_, mapNormal, gblIdx);
 
             if (mag(normal) != 0 && j != 0)
             {
@@ -328,7 +331,6 @@ Foam::reconstruction::plicRDF::plicRDF
     iteration_(modelDict().getOrDefault<label>("iterations", 5)),
     interpolateNormal_(modelDict().getOrDefault("interpolateNormal", true)),
     RDF_(mesh_),
-    exchangeFields_(zoneDistribute::New(mesh_)),
     sIterPLIC_(mesh_,surfCellTol_)
 {
     setInitNormals(false);
@@ -375,6 +377,7 @@ Foam::reconstruction::plicRDF::plicRDF
 
 void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
 {
+    zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
     const bool uptodate = alreadyReconstructed(forceUpdate);
 
     if (uptodate && !forceUpdate)
@@ -403,7 +406,7 @@ void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
 
     bitSet tooCoarse(mesh_.nCells(),false);
 
-    for (int iter=0; iter<iteration_; ++iter)
+    for (label iter=0; iter<iteration_; ++iter)
     {
         forAll(interfaceLabels_, i)
         {
@@ -452,7 +455,7 @@ void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
                 nextToInterface_,
                 centre_,
                 normal_,
-                exchangeFields_,
+                exchangeFields,
                 false
             );
             RDF_.updateContactAngle(alpha1_, U_, nHatb);
@@ -492,7 +495,6 @@ void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
                 resCounter++;
 
             }
-
         }
 
         reduce(avgRes,sumOp<scalar>());
diff --git a/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.H b/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.H
index cbdb19d8cf..02a1bdf139 100644
--- a/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.H
+++ b/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.H
@@ -118,9 +118,6 @@ class plicRDF
         //- Calculates the RDF function
         reconstructedDistanceFunction RDF_;
 
-        //- Provides stencil and map
-        zoneDistribute& exchangeFields_;
-
         //- surfaceIterator finds the plane centre for specified VOF value
         surfaceIteratorPLIC sIterPLIC_;
 
diff --git a/src/transportModels/geometricVoF/reconstructionSchemes/reconstructionSchemes.C b/src/transportModels/geometricVoF/reconstructionSchemes/reconstructionSchemes.C
index 78f81981a0..5ff4c82482 100644
--- a/src/transportModels/geometricVoF/reconstructionSchemes/reconstructionSchemes.C
+++ b/src/transportModels/geometricVoF/reconstructionSchemes/reconstructionSchemes.C
@@ -106,11 +106,13 @@ Foam::reconstructionSchemes::reconstructionSchemes
     (
         IOobject
         (
-            "recon::normal",
+            IOobject::groupName("interfaceNormal", alpha1.group()),
             alpha1_.mesh().time().timeName(),
             alpha1_.mesh(),
             IOobject::NO_READ,
-            IOobject::AUTO_WRITE
+            dict.getOrDefault("writeFields",false)
+          ? IOobject::AUTO_WRITE
+          : IOobject::NO_WRITE
         ),
         alpha1_.mesh(),
         dimensionedVector(dimArea, Zero)
@@ -119,11 +121,13 @@ Foam::reconstructionSchemes::reconstructionSchemes
     (
         IOobject
         (
-            "recon::centre",
+            IOobject::groupName("interfaceCentre", alpha1.group()),
             alpha1_.mesh().time().timeName(),
             alpha1_.mesh(),
             IOobject::NO_READ,
-            IOobject::AUTO_WRITE
+            dict.getOrDefault("writeFields",false)
+          ? IOobject::AUTO_WRITE
+          : IOobject::NO_WRITE
         ),
         alpha1_.mesh(),
         dimensionedVector(dimLength, Zero)
diff --git a/src/transportModels/geometricVoF/reconstructionSchemes/reconstructionSchemes.H b/src/transportModels/geometricVoF/reconstructionSchemes/reconstructionSchemes.H
index b2b5d900b3..5ea5a43486 100644
--- a/src/transportModels/geometricVoF/reconstructionSchemes/reconstructionSchemes.H
+++ b/src/transportModels/geometricVoF/reconstructionSchemes/reconstructionSchemes.H
@@ -209,13 +209,13 @@ public:
     // Member Functions
 
         //- Access to the model dictionary
-        dictionary& modelDict()
+        dictionary& modelDict() noexcept
         {
             return reconstructionSchemesCoeffs_;
         }
 
         //- Const access to the model dictionary
-        const dictionary& modelDict() const
+        const dictionary& modelDict() const noexcept
         {
             return reconstructionSchemesCoeffs_;
         }
@@ -223,26 +223,38 @@ public:
         //- Reconstruct the interface
         virtual void reconstruct(bool forceUpdate = true) = 0;
 
+        //- const-Reference to interface normals
+        const volVectorField& normal() const noexcept
+        {
+            return normal_;
+        }
+
+        //- const-Reference to interface centres
+        const volVectorField& centre() const noexcept
+        {
+            return centre_;
+        }
+
         //- Reference to interface normals
-        const volVectorField& normal() const
+        volVectorField& normal() noexcept
         {
             return normal_;
         }
 
         //- Reference to interface centres
-        const volVectorField& centre() const
+        volVectorField& centre() noexcept
         {
             return centre_;
         }
 
         //- Does the cell contain interface
-        const boolList& interfaceCell() const
+        const boolList& interfaceCell() const noexcept
         {
             return interfaceCell_;
         }
 
         //- List of cells with an interface
-        const DynamicField<label>& interfaceLabels() const
+        const DynamicField<label>& interfaceLabels() const noexcept
         {
             return interfaceLabels_;
         }
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/system/fvSolution b/tutorials/multiphase/interIsoFoam/damBreak/system/fvSolution
index 580eb9aee1..b7f08b9a11 100644
--- a/tutorials/multiphase/interIsoFoam/damBreak/system/fvSolution
+++ b/tutorials/multiphase/interIsoFoam/damBreak/system/fvSolution
@@ -25,6 +25,7 @@ solvers
         snapTol         1e-12;
         clip            true;
         reconstructionScheme isoAlpha;
+        writeFields     true;
 
         nAlphaSubCycles 1;
         cAlpha          1; // Note: cAlpha is not used by isoAdvector but must