diff --git a/src/dynamicFaMesh/interfaceTrackingFvMesh/freeSurfacePointDisplacement.C b/src/dynamicFaMesh/interfaceTrackingFvMesh/freeSurfacePointDisplacement.C
index c11c4f91bb17360133cefa103f300574eba4fcaf..05ac7f392f0cfe7aafab9a1a3f2e2baabf9b7928 100644
--- a/src/dynamicFaMesh/interfaceTrackingFvMesh/freeSurfacePointDisplacement.C
+++ b/src/dynamicFaMesh/interfaceTrackingFvMesh/freeSurfacePointDisplacement.C
@@ -315,7 +315,7 @@ Foam::interfaceTrackingFvMesh::pointDisplacement()
     }
 
     // Calculate displacement of processor patch points
-    forAll (aMesh().boundary(), patchI)
+    forAll(aMesh().boundary(), patchI)
     {
         if
         (
diff --git a/src/functionObjects/field/externalCoupled/externalCoupledTemperatureMixed/externalCoupledTemperatureMixedFvPatchScalarField.C b/src/functionObjects/field/externalCoupled/externalCoupledTemperatureMixed/externalCoupledTemperatureMixedFvPatchScalarField.C
index 2d5209d7b1a6bba1d62e438b42ff4ce1813cea49..2500429a6e3199b31546b13ba485ebd51fd5bcea 100644
--- a/src/functionObjects/field/externalCoupled/externalCoupledTemperatureMixed/externalCoupledTemperatureMixedFvPatchScalarField.C
+++ b/src/functionObjects/field/externalCoupled/externalCoupledTemperatureMixed/externalCoupledTemperatureMixedFvPatchScalarField.C
@@ -286,10 +286,10 @@ void Foam::externalCoupledTemperatureMixedFvPatchScalarField::writeData
     // Heat transfer coefficient [W/m2/K]
     // This htc needs to be always larger or equal to zero
     //const scalarField htc(qDot/max(Twall - Tfluid, 1e-3));
-    scalarField htc(qDot.size(), 0);
-    forAll (qDot, i)
+    scalarField htc(qDot.size(), Zero);
+    forAll(qDot, i)
     {
-        scalar deltaT = mag(Twall[i] - Tfluid[i]);
+        const scalar deltaT = mag(Twall[i] - Tfluid[i]);
         if (deltaT > 1e-3)
         {
             htc[i] = sign(qDot[i])*qDot[i]/deltaT;
diff --git a/src/meshTools/sets/faceZoneSources/planeToFaceZone/planeToFaceZone.C b/src/meshTools/sets/faceZoneSources/planeToFaceZone/planeToFaceZone.C
deleted file mode 100644
index e316ba0cc396fd4aee486797e4c8f1bc4c699861..0000000000000000000000000000000000000000
--- a/src/meshTools/sets/faceZoneSources/planeToFaceZone/planeToFaceZone.C
+++ /dev/null
@@ -1,413 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenFOAM Foundation
-    Copyright (C) 2020 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 "planeToFaceZone.H"
-#include "polyMesh.H"
-#include "faceZoneSet.H"
-#include "uindirectPrimitivePatch.H"
-#include "PatchTools.H"
-#include "syncTools.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(planeToFaceZone, 0);
-    addToRunTimeSelectionTable(topoSetSource, planeToFaceZone, word);
-    addToRunTimeSelectionTable(topoSetSource, planeToFaceZone, istream);
-}
-
-
-Foam::topoSetSource::addToUsageTable Foam::planeToFaceZone::usage_
-(
-    planeToFaceZone::typeName,
-    "\n    Usage: planeToFaceZone (px py pz) (nx ny nz) include\n\n"
-    "    Select faces for which the adjacent cell centres lie on opposite "
-    " of a plane\n\n"
-);
-
-
-const Foam::Enum
-<
-    Foam::planeToFaceZone::faceZoneAction
->
-Foam::planeToFaceZone::faceZoneActionNames_
-({
-    { faceZoneAction::ALL, "all" },
-    { faceZoneAction::CLOSEST, "closest" },
-});
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::planeToFaceZone::combine(faceZoneSet& fzSet, const bool add) const
-{
-    // Mark all cells with centres above the plane
-    boolList cellIsAbovePlane(mesh_.nCells());
-    forAll(mesh_.cells(), celli)
-    {
-        cellIsAbovePlane[celli] =
-            ((mesh_.cellCentres()[celli] - point_) & normal_) > 0;
-    }
-
-    // Mark all faces that sit between cells above and below the plane
-    boolList faceIsOnPlane(mesh_.nFaces());
-    forAll(mesh_.faceNeighbour(), facei)
-    {
-        faceIsOnPlane[facei] =
-            cellIsAbovePlane[mesh_.faceOwner()[facei]]
-         != cellIsAbovePlane[mesh_.faceNeighbour()[facei]];
-    }
-    forAll(mesh_.boundaryMesh(), patchi)
-    {
-        const polyPatch& patch = mesh_.boundaryMesh()[patchi];
-        forAll(patch, patchFacei)
-        {
-            const label facei = patch.start() + patchFacei;
-            faceIsOnPlane[facei] =
-                patch.coupled() && cellIsAbovePlane[mesh_.faceOwner()[facei]];
-        }
-    }
-    syncTools::syncFaceList(mesh_, faceIsOnPlane, notEqOp<bool>());
-
-    // Convert marked faces to a list of indices
-    labelList newSetFaces(findIndices(faceIsOnPlane, true));
-
-    // If constructing a single contiguous set, remove all faces except those
-    // connected to the contiguous region closest to the specified point
-    if (option_ == faceZoneAction::CLOSEST)
-    {
-        // Step 1: Get locally contiguous regions for the new face set and the
-        // total number of regions across all processors.
-        labelList newSetFaceRegions(newSetFaces.size(), -1);
-        label nRegions = -1;
-        {
-            // Create a patch of the set faces
-            const uindirectPrimitivePatch newSetPatch
-            (
-                UIndirectList<face>(mesh_.faces(), newSetFaces),
-                mesh_.points()
-            );
-
-            // Get the region ID-s and store the total number of regions on
-            // each processor
-            labelList procNRegions(Pstream::nProcs(), -1);
-            procNRegions[Pstream::myProcNo()] =
-                PatchTools::markZones
-                (
-                    newSetPatch,
-                    boolList(newSetPatch.nEdges(), false),
-                    newSetFaceRegions
-                );
-            Pstream::gatherList(procNRegions);
-            Pstream::scatterList(procNRegions);
-
-            // Cumulative sum the number of regions on each processor to get an
-            // offset which makes the local region ID-s globally unique
-            labelList procRegionOffset(Pstream::nProcs(), 0);
-            for (label proci = 1; proci < Pstream::nProcs(); ++ proci)
-            {
-                procRegionOffset[proci] +=
-                    procRegionOffset[proci - 1]
-                  + procNRegions[proci - 1];
-            }
-
-            // Apply the offset
-            forAll(newSetFaces, newSetFacei)
-            {
-                newSetFaceRegions[newSetFacei] +=
-                    procRegionOffset[Pstream::myProcNo()];
-            }
-
-            // Store the total number of regions across all processors
-            nRegions = procRegionOffset.last() + procNRegions.last();
-        }
-
-        // Step 2: Create a region map which combines regions which are
-        // connected across coupled interfaces
-        labelList regionMap(identity(nRegions));
-        {
-            // Put region labels on connected boundary edges and synchronise to
-            // create a list of all regions connected to a given edge
-            labelListList meshEdgeRegions(mesh_.nEdges(), labelList());
-            forAll(newSetFaces, newSetFacei)
-            {
-                const label facei = newSetFaces[newSetFacei];
-                const label regioni = newSetFaceRegions[newSetFacei];
-                forAll(mesh_.faceEdges()[facei], faceEdgei)
-                {
-                    const label edgei = mesh_.faceEdges()[facei][faceEdgei];
-                    meshEdgeRegions[edgei] = labelList(1, regioni);
-                }
-            }
-            syncTools::syncEdgeList
-            (
-                mesh_,
-                meshEdgeRegions,
-                globalMeshData::ListPlusEqOp<labelList>(),
-                labelList()
-            );
-
-            // Combine edge regions to create a list of what regions a given
-            // region is connected to
-            List<labelHashSet> regionRegions(nRegions);
-            forAll(newSetFaces, newSetFacei)
-            {
-                const label facei = newSetFaces[newSetFacei];
-                const label regioni = newSetFaceRegions[newSetFacei];
-                forAll(mesh_.faceEdges()[facei], faceEdgei)
-                {
-                    const label edgei = mesh_.faceEdges()[facei][faceEdgei];
-                    forAll(meshEdgeRegions[edgei], edgeRegioni)
-                    {
-                        if (meshEdgeRegions[edgei][edgeRegioni] != regioni)
-                        {
-                            regionRegions[regioni].insert
-                            (
-                                meshEdgeRegions[edgei][edgeRegioni]
-                            );
-                        }
-                    }
-                }
-            }
-            Pstream::listCombineGather(regionRegions, plusEqOp<labelHashSet>());
-            Pstream::listCombineScatter(regionRegions);
-
-            // Collapse the region connections into a map between each region
-            // and the lowest numbered region that it connects to
-            forAll(regionRegions, regioni)
-            {
-                forAllConstIter(labelHashSet, regionRegions[regioni], iter)
-                {
-                    regionMap[iter.key()] =
-                        min(regionMap[iter.key()], regionMap[regioni]);
-                }
-            }
-        }
-
-        // Step 3: Combine connected regions
-        labelList regionNFaces;
-        {
-            // Remove duplicates from the region map
-            label regioni0 = 0;
-            forAll(regionMap, regioni)
-            {
-                if (regionMap[regioni] > regioni0)
-                {
-                    ++ regioni0;
-                    regionMap[regioni] = regioni0;
-                }
-            }
-
-            // Recompute the number of regions
-            nRegions = regioni0 + 1;
-
-            // Renumber the face region ID-s
-            newSetFaceRegions =
-                IndirectList<label>(regionMap, newSetFaceRegions);
-
-            // Report the final number and size of the regions
-            regionNFaces = labelList(nRegions, 0);
-            forAll(newSetFaces, newSetFacei)
-            {
-                regionNFaces[newSetFaceRegions[newSetFacei]] ++;
-            }
-            Pstream::listCombineGather(regionNFaces, plusEqOp<label>());
-            Pstream::listCombineScatter(regionNFaces);
-            Info<< "    Found " << nRegions << " contiguous regions with "
-                << regionNFaces << " faces" << endl;
-        }
-
-        // Step 4: Choose the closest region to output
-        label selectedRegioni = -1;
-        {
-            // Compute the region centres
-            scalarField regionMagAreas(nRegions, 0);
-            pointField regionCentres(nRegions, Zero);
-            forAll(newSetFaces, newSetFacei)
-            {
-                const label facei = newSetFaces[newSetFacei];
-                const label regioni = newSetFaceRegions[newSetFacei];
-
-                const vector& a = mesh_.faceAreas()[facei];
-                const point& c = mesh_.faceCentres()[facei];
-
-                regionMagAreas[regioni] += mag(a);
-                regionCentres[regioni] += mag(a)*c;
-            }
-            Pstream::listCombineGather(regionMagAreas, plusEqOp<scalar>());
-            Pstream::listCombineGather(regionCentres, plusEqOp<point>());
-            Pstream::listCombineScatter(regionMagAreas);
-            Pstream::listCombineScatter(regionCentres);
-            regionCentres /= regionMagAreas;
-
-            // Find the region centroid closest to the reference point
-            selectedRegioni =
-                returnReduce
-                (
-                    findMin(mag(regionCentres - point_)()),
-                    minOp<label>()
-                );
-
-            // Report the selection
-            Info<< "    Selecting region " << selectedRegioni << " with "
-                << regionNFaces[selectedRegioni]
-                << " faces as the closest to point " << point_ << endl;
-        }
-
-        // Step 5: Remove any faces from the set list not in the selected region
-        {
-            // Remove faces from the list by shuffling up and resizing
-            label newSetFacei0 = 0;
-            forAll(newSetFaces, newSetFacei)
-            {
-                newSetFaces[newSetFacei0] = newSetFaces[newSetFacei];
-
-                if (newSetFaceRegions[newSetFacei] == selectedRegioni)
-                {
-                    ++ newSetFacei0;
-                }
-            }
-            newSetFaces.resize(newSetFacei0);
-        }
-    }
-
-    // Modify the face zone set
-    DynamicList<label> newAddressing;
-    DynamicList<bool> newFlipMap;
-    if (add)
-    {
-        // Start from copy
-        newAddressing = DynamicList<label>(fzSet.addressing());
-        newFlipMap = DynamicList<bool>(fzSet.flipMap());
-
-        // Add anything from the new set that is not already in the zone set
-        forAll(newSetFaces, newSetFacei)
-        {
-            const label facei = newSetFaces[newSetFacei];
-
-            if (!fzSet.found(facei))
-            {
-                newAddressing.append(facei);
-                newFlipMap.append(cellIsAbovePlane[mesh_.faceOwner()[facei]]);
-            }
-        }
-    }
-    else
-    {
-        // Start from empty
-        newAddressing = DynamicList<label>(fzSet.addressing().size());
-        newFlipMap = DynamicList<bool>(fzSet.flipMap().size());
-
-        // Add everything from the zone set that is not also in the new set
-        labelHashSet newSet(newSetFaces);
-        forAll(fzSet.addressing(), i)
-        {
-            const label facei = fzSet.addressing()[i];
-
-            if (!newSet.found(facei))
-            {
-                newAddressing.append(facei);
-                newFlipMap.append(cellIsAbovePlane[mesh_.faceOwner()[facei]]);
-            }
-        }
-    }
-    fzSet.addressing().transfer(newAddressing);
-    fzSet.flipMap().transfer(newFlipMap);
-    fzSet.updateSet();
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::planeToFaceZone::planeToFaceZone
-(
-    const polyMesh& mesh,
-    const dictionary& dict
-)
-:
-    topoSetSource(mesh),
-    point_(dict.get<vector>("point")),
-    normal_(dict.get<vector>("normal")),
-    option_
-    (
-        faceZoneActionNames_.getOrDefault("option", dict, faceZoneAction::ALL)
-    )
-{}
-
-
-Foam::planeToFaceZone::planeToFaceZone
-(
-    const polyMesh& mesh,
-    Istream& is
-)
-:
-    topoSetSource(mesh),
-    point_(checkIs(is)),
-    normal_(checkIs(is)),
-    option_(faceZoneActionNames_.read(checkIs(is)))
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::planeToFaceZone::applyToSet
-(
-    const topoSetSource::setAction action,
-    topoSet& set
-) const
-{
-    if (!isA<faceZoneSet>(set))
-    {
-        WarningInFunction
-            << "Operation only allowed on a faceZoneSet." << endl;
-    }
-    else
-    {
-        faceZoneSet& fzSet = refCast<faceZoneSet>(set);
-
-        if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
-        {
-            Info<< "    Adding faces which form a plane at " << point_
-                << " with normal " << normal_ << endl;
-
-            combine(fzSet, true);
-        }
-        else if (action == topoSetSource::DELETE)
-        {
-            Info<< "    Removing faces which form a plane at " << point_
-                << " with normal " << normal_ << endl;
-
-            combine(fzSet, false);
-        }
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DVolume/volBSplinesBase/volBSplinesBase.C b/src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DVolume/volBSplinesBase/volBSplinesBase.C
index 4f8a5f94259e5013dcf1af1220102bfb8e429377..462e9e077a26204d668c0cf9d235ff41e04985c3 100644
--- a/src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DVolume/volBSplinesBase/volBSplinesBase.C
+++ b/src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DVolume/volBSplinesBase/volBSplinesBase.C
@@ -260,7 +260,7 @@ void Foam::volBSplinesBase::boundControlPointMovement
         vectorField localControlPointsMovement(nb, Zero);
 
         // Set localControlPointsMovement
-        forAll (localControlPointsMovement, iCPM)
+        forAll(localControlPointsMovement, iCPM)
         {
             localControlPointsMovement[iCPM] =
                 controlPointsMovement[pastControlPoints + iCPM];
@@ -292,7 +292,7 @@ void Foam::volBSplinesBase::moveControlPoints
         vectorField localControlPointsMovement(nb, Zero);
 
         // Set localControlPointsMovement
-        forAll (localControlPointsMovement, iCPM)
+        forAll(localControlPointsMovement, iCPM)
         {
             localControlPointsMovement[iCPM] =
                 controlPointsMovement[pastControlPoints + iCPM];
diff --git a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
index ba17ff0a2fe389521c1b85f7cfb06a5c1321f906..2d673c7d2307b8e9a59d81cdf43141d3f6a38ca5 100644
--- a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
+++ b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
@@ -468,7 +468,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
             // Convective thermal diffusivity for single phase
             const scalarField alphatv(calcAlphat(*this));
 
-            forAll (*this, i)
+            forAll(*this, i)
             {
                 this->operator[](i) =
                 (
@@ -802,7 +802,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
             // Use to identify regimes per face
             labelField regimeTypes(A1.size(), -1);
 
-            forAll (*this, i)
+            forAll(*this, i)
             {
                 if (Tw[i] > Tsatw[i])
                 {
@@ -1004,7 +1004,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     scalarField nFilms(this->size(), 0);
                     scalarField nNonBoilings(this->size(), 0);
 
-                    forAll (*this, i)
+                    forAll(*this, i)
                     {
                         //faceRegimes[i] = regimeTypes[i];
                         switch (regimeTypes[i])
diff --git a/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C b/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C
index 7aba9c546de0ae1faf672ba56ea6a23f3270344c..30777d8bd0496908f0ae9408ca0f1a36b011fa71 100644
--- a/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C
+++ b/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C
@@ -280,7 +280,7 @@ const Foam::volScalarField&  Foam::reconstructedDistanceFunction::constructRDF
                 scalar avgWeight = 0;
                 const point p = mesh_.C()[celli];
 
-                forAll (stencil[celli],i)
+                forAll(stencil[celli], i)
                 {
                     const label gblIdx = stencil[celli][i];
                     vector n = -distribute.getValue(normal, mapNormal, gblIdx);
@@ -334,7 +334,7 @@ const Foam::volScalarField&  Foam::reconstructedDistanceFunction::constructRDF
                     scalar avgWeight = 0;
                     const point p = mesh_.C().boundaryField()[patchI][i];
 
-                    forAll (stencil[pCellI], j)
+                    forAll(stencil[pCellI], j)
                     {
                         const label gblIdx = stencil[pCellI][j];
                         vector n = -distribute.getValue(normal, mapNormal, gblIdx);