diff --git a/src/fileFormats/ensight/mesh/ensightMesh.C b/src/fileFormats/ensight/mesh/ensightMesh.C
index b889b524b5..8681c14abb 100644
--- a/src/fileFormats/ensight/mesh/ensightMesh.C
+++ b/src/fileFormats/ensight/mesh/ensightMesh.C
@@ -314,7 +314,7 @@ void Foam::ensightMesh::correct()
     }
 
 
-    // Face exclusion for empty/processor types
+    // Face exclusion for empty types and neighbour side of coupled
     // so they are ignored for face zones
 
     if (fzoneIds.size())
@@ -324,15 +324,16 @@ void Foam::ensightMesh::correct()
 
         for (const polyPatch& p : mesh_.boundaryMesh())
         {
-            const auto* procPatch = isA<processorPolyPatch>(p);
+            const auto* cpp = isA<coupledPolyPatch>(p);
 
-            if (isA<emptyPolyPatch>(p))
+            if
+            (
+                isA<emptyPolyPatch>(p)
+             || (cpp && !cpp->owner())
+            )
             {
-                excludeFace.set(p.range());
-            }
-            else if (procPatch && !procPatch->owner())
-            {
-                // Exclude neighbour-side, retain owner-side only
+                // Ignore empty patch
+                // Ignore neighbour side of coupled
                 excludeFace.set(p.range());
             }
         }
diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
index 5ddf3663c2..22d691e199 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
@@ -168,9 +168,9 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::setFaceZoneFaces()
     }
     #endif
 
-    faceId_.resize(numFaces);
-    facePatchId_.resize(numFaces);
-    faceFlip_.resize(numFaces);
+    faceId_.resize_nocopy(numFaces);
+    facePatchId_.resize_nocopy(numFaces);
+    faceFlip_.resize_nocopy(numFaces);
 
     numFaces = 0;
 
@@ -192,21 +192,20 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::setFaceZoneFaces()
             {
                 facePatchId = mesh_.boundaryMesh().whichPatch(meshFacei);
                 const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
+
+                if (isA<emptyPolyPatch>(pp))
+                {
+                    continue;  // Ignore empty patch
+                }
+
                 const auto* cpp = isA<coupledPolyPatch>(pp);
 
-                if (cpp)
+                if (cpp && !cpp->owner())
                 {
-                    faceId = (cpp->owner() ? pp.whichFace(meshFacei) : -1);
-                }
-                else if (!isA<emptyPolyPatch>(pp))
-                {
-                    faceId = pp.whichFace(meshFacei);
-                }
-                else
-                {
-                    faceId = -1;
-                    facePatchId = -1;
+                    continue;  // Ignore neighbour side
                 }
+
+                faceId = pp.whichFace(meshFacei);
             }
 
             if (faceId >= 0)
diff --git a/src/functionObjects/field/fluxSummary/fluxSummary.C b/src/functionObjects/field/fluxSummary/fluxSummary.C
index d87ab130e6..2969be6906 100644
--- a/src/functionObjects/field/fluxSummary/fluxSummary.C
+++ b/src/functionObjects/field/fluxSummary/fluxSummary.C
@@ -201,56 +201,62 @@ void Foam::functionObjects::fluxSummary::initialiseFaceZone
     names.append(faceZoneName);
     directions.append(Zero); // dummy value
 
-    DynamicList<label> faceIDs(fZone.size());
-    DynamicList<label> facePatchIDs(fZone.size());
-    DynamicList<bool>  flips(fZone.size());
+    labelList faceIds(fZone.size());
+    labelList facePatchIds(fZone.size());
+    boolList  faceFlips(fZone.size());
+
+    // Total number of faces selected
+    label numFaces = 0;
 
     forAll(fZone, i)
     {
-        label facei = fZone[i];
+        const label meshFacei = fZone[i];
         const bool isFlip = fZone.flipMap()[i];
 
-        label faceID = -1;
-        label facePatchID = -1;
-        if (mesh_.isInternalFace(facei))
+        // Internal faces
+        label faceId = meshFacei;
+        label facePatchId = -1;
+
+        // Boundary faces
+        if (!mesh_.isInternalFace(meshFacei))
         {
-            faceID = facei;
-            facePatchID = -1;
-        }
-        else
-        {
-            facePatchID = mesh_.boundaryMesh().whichPatch(facei);
-            const polyPatch& pp = mesh_.boundaryMesh()[facePatchID];
+            facePatchId = mesh_.boundaryMesh().whichPatch(meshFacei);
+            const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
+
+            if (isA<emptyPolyPatch>(pp))
+            {
+                continue;  // Ignore empty patch
+            }
+
             const auto* cpp = isA<coupledPolyPatch>(pp);
 
-            if (cpp)
+            if (cpp && !cpp->owner())
             {
-                faceID = (cpp->owner() ? pp.whichFace(facei) : -1);
-            }
-            else if (!isA<emptyPolyPatch>(pp))
-            {
-                faceID = pp.whichFace(facei);
-            }
-            else
-            {
-                faceID = -1;
-                facePatchID = -1;
+                continue;  // Ignore neighbour side
             }
+
+            faceId = pp.whichFace(meshFacei);
         }
 
-        if (faceID >= 0)
+        if (faceId >= 0)
         {
             // Orientation set by faceZone flip map
-            flips.append(isFlip);
-            faceIDs.append(faceID);
-            facePatchIDs.append(facePatchID);
+            faceIds[numFaces] = faceId;
+            facePatchIds[numFaces] = facePatchId;
+            faceFlips[numFaces] = isFlip;
+
+            ++numFaces;
         }
     }
 
-    // could reduce some copying here
-    faceID.append(faceIDs);
-    facePatchID.append(facePatchIDs);
-    faceFlip.append(flips);
+    // Shrink to size used
+    faceIds.resize(numFaces);
+    facePatchIds.resize(numFaces);
+    faceFlips.resize(numFaces);
+
+    faceID.append(std::move(faceIds));
+    facePatchID.append(std::move(facePatchIds));
+    faceFlip.append(std::move(faceFlips));
 }
 
 
@@ -281,78 +287,84 @@ void Foam::functionObjects::fluxSummary::initialiseFaceZoneAndDirection
     names.append(faceZoneName);
     directions.append(refDir);
 
-    DynamicList<label> faceIDs(fZone.size());
-    DynamicList<label> facePatchIDs(fZone.size());
-    DynamicList<bool>  flips(fZone.size());
+    labelList faceIds(fZone.size());
+    labelList facePatchIds(fZone.size());
+    boolList  faceFlips(fZone.size());
 
     const surfaceVectorField& Sf = mesh_.Sf();
     const surfaceScalarField& magSf = mesh_.magSf();
 
     vector n(Zero);
 
+    // Total number of faces selected
+    label numFaces = 0;
+
     forAll(fZone, i)
     {
-        label facei = fZone[i];
+        const label meshFacei = fZone[i];
 
-        label faceID = -1;
-        label facePatchID = -1;
-        if (mesh_.isInternalFace(facei))
+        // Internal faces
+        label faceId = meshFacei;
+        label facePatchId = -1;
+
+        // Boundary faces
+        if (!mesh_.isInternalFace(meshFacei))
         {
-            faceID = facei;
-            facePatchID = -1;
-        }
-        else
-        {
-            facePatchID = mesh_.boundaryMesh().whichPatch(facei);
-            const polyPatch& pp = mesh_.boundaryMesh()[facePatchID];
+            facePatchId = mesh_.boundaryMesh().whichPatch(meshFacei);
+            const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
+
+            if (isA<emptyPolyPatch>(pp))
+            {
+                continue;  // Ignore empty patch
+            }
+
             const auto* cpp = isA<coupledPolyPatch>(pp);
 
-            if (cpp)
+            if (cpp && !cpp->owner())
             {
-                faceID = (cpp->owner() ? pp.whichFace(facei) : -1);
-            }
-            else if (!isA<emptyPolyPatch>(pp))
-            {
-                faceID = pp.whichFace(facei);
-            }
-            else
-            {
-                faceID = -1;
-                facePatchID = -1;
+                continue;  // Ignore neighbour side
             }
+
+            faceId = pp.whichFace(meshFacei);
         }
 
-        if (faceID >= 0)
+        if (faceId >= 0)
         {
-            // orientation set by comparison with reference direction
-            if (facePatchID != -1)
+            // Orientation set by comparison with reference direction
+            if (facePatchId != -1)
             {
-                n = Sf.boundaryField()[facePatchID][faceID]
-                   /(magSf.boundaryField()[facePatchID][faceID] + ROOTVSMALL);
+                n = Sf.boundaryField()[facePatchId][faceId]
+                   /(magSf.boundaryField()[facePatchId][faceId] + ROOTVSMALL);
             }
             else
             {
-                n = Sf[faceID]/(magSf[faceID] + ROOTVSMALL);
+                n = Sf[faceId]/(magSf[faceId] + ROOTVSMALL);
             }
 
             if ((n & refDir) > tolerance_)
             {
-                flips.append(false);
+                faceFlips[numFaces] = false;
             }
             else
             {
-                flips.append(true);
+                faceFlips[numFaces] = true;
             }
 
-            faceIDs.append(faceID);
-            facePatchIDs.append(facePatchID);
+            faceIds[numFaces] = faceId;
+            facePatchIds[numFaces] = facePatchId;
+
+            ++numFaces;
         }
     }
 
-    // could reduce copying here
-    faceID.append(faceIDs);
-    facePatchID.append(facePatchIDs);
-    faceFlip.append(flips);
+    // Shrink to size used
+    faceIds.resize(numFaces);
+    facePatchIds.resize(numFaces);
+    faceFlips.resize(numFaces);
+
+    faceID.append(std::move(faceIds));
+    facePatchID.append(std::move(facePatchIds));
+    faceFlip.append(std::move(faceFlips));
 }
 
 
@@ -498,7 +510,7 @@ void Foam::functionObjects::fluxSummary::initialiseCellZoneAndDirection
         << "Starting walk to split patch into faceZones"
         << endl;
 
-    globalIndex globalFaces(patch.size());
+    const globalIndex globalFaces(patch.size());
 
     label oldFaceID = 0;
     label regioni = 0;
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/faceZoneReferenceTemperature/faceZoneReferenceTemperature.C b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/faceZoneReferenceTemperature/faceZoneReferenceTemperature.C
index 4cd8854a94..66f6a390a6 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/faceZoneReferenceTemperature/faceZoneReferenceTemperature.C
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/faceZoneReferenceTemperature/faceZoneReferenceTemperature.C
@@ -70,7 +70,7 @@ setFaceZoneFaces(const dictionary& dict)
 
     label numFaces = fZone.size();
 
-    if (!returnReduce(bool(numFaces), orOp<bool>()))
+    if (!returnReduceOr(numFaces))
     {
         FatalIOErrorInFunction(dict)
             << "referenceFaceZone: " << faceZoneName
@@ -78,49 +78,55 @@ setFaceZoneFaces(const dictionary& dict)
             << exit(FatalIOError);
     }
 
-    faceId_.resize(numFaces);
-    facePatchId_.resize(numFaces);
+    faceId_.resize_nocopy(numFaces);
+    facePatchId_.resize_nocopy(numFaces);
 
     numFaces = 0;
 
-    forAll(fZone, i)
+    // TDB: handle multiple zones
     {
-        const label meshFacei = fZone[i];
-
-        // Internal faces
-        label faceId = meshFacei;
-        label facePatchId = -1;
-
-        // Boundary faces
-        if (!mesh.isInternalFace(meshFacei))
+        forAll(fZone, i)
         {
-            facePatchId = mesh.boundaryMesh().whichPatch(meshFacei);
-            const polyPatch& pp = mesh.boundaryMesh()[facePatchId];
-            const auto* cpp = isA<coupledPolyPatch>(pp);
+            const label meshFacei = fZone[i];
 
-            if (cpp)
-            {
-                faceId = (cpp->owner() ? pp.whichFace(meshFacei) : -1);
-            }
-            else if (!isA<emptyPolyPatch>(pp))
+            // Internal faces
+            label faceId = meshFacei;
+            label facePatchId = -1;
+
+            // Boundary faces
+            if (!mesh.isInternalFace(meshFacei))
             {
+                facePatchId = mesh.boundaryMesh().whichPatch(meshFacei);
+                const polyPatch& pp = mesh.boundaryMesh()[facePatchId];
+
+                if (isA<emptyPolyPatch>(pp))
+                {
+                    continue;  // Ignore empty patch
+                }
+
+                const auto* cpp = isA<coupledPolyPatch>(pp);
+
+                if (cpp && !cpp->owner())
+                {
+                    continue;  // Ignore neighbour side
+                }
+
                 faceId = pp.whichFace(meshFacei);
             }
-            else
+
+            if (faceId >= 0)
             {
-                faceId = -1;
-                facePatchId = -1;
+                faceId_[numFaces] = faceId;
+                facePatchId_[numFaces] = facePatchId;
+
+                ++numFaces;
             }
         }
-
-        if (faceId >= 0)
-        {
-            faceId_[numFaces] = faceId;
-            facePatchId_[numFaces] = facePatchId;
-
-            ++numFaces;
-        }
     }
+
+    // Shrink to size used
+    faceId_.resize(numFaces);
+    facePatchId_.resize(numFaces);
 }
 
 
diff --git a/src/fvOptions/sources/derived/directionalPressureGradientExplicitSource/directionalPressureGradientExplicitSource.C b/src/fvOptions/sources/derived/directionalPressureGradientExplicitSource/directionalPressureGradientExplicitSource.C
index 3a0b6fe59c..7671078e5e 100644
--- a/src/fvOptions/sources/derived/directionalPressureGradientExplicitSource/directionalPressureGradientExplicitSource.C
+++ b/src/fvOptions/sources/derived/directionalPressureGradientExplicitSource/directionalPressureGradientExplicitSource.C
@@ -76,51 +76,58 @@ void Foam::fv::directionalPressureGradientExplicitSource::initialise()
 {
     const faceZone& fZone = mesh_.faceZones()[zoneID_];
 
-    faceId_.setSize(fZone.size());
-    facePatchId_.setSize(fZone.size());
+    // Total number of faces selected
+    label numFaces = fZone.size();
 
-    label count = 0;
-    forAll(fZone, i)
+    faceId_.resize_nocopy(numFaces);
+    facePatchId_.resize_nocopy(numFaces);
+
+    numFaces = 0;
+
+    // TDB: handle multiple zones
     {
-        const label faceI = fZone[i];
-
-        label faceId = -1;
-        label facePatchId = -1;
-        if (mesh_.isInternalFace(faceI))
+        forAll(fZone, i)
         {
-            faceId = faceI;
-            facePatchId = -1;
-        }
-        else
-        {
-            facePatchId = mesh_.boundaryMesh().whichPatch(faceI);
-            const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
-            const auto* cpp = isA<coupledPolyPatch>(pp);
+            const label meshFacei = fZone[i];
 
-            if (cpp)
-            {
-                faceId = (cpp->owner() ? pp.whichFace(faceI) : -1);
-            }
-            else if (!isA<emptyPolyPatch>(pp))
-            {
-                faceId = pp.whichFace(faceI);
-            }
-            else
-            {
-                faceId = -1;
-                facePatchId = -1;
-            }
-        }
+            // Internal faces
+            label faceId = meshFacei;
+            label facePatchId = -1;
 
-        if (faceId >= 0)
-        {
-            facePatchId_[count] = facePatchId;
-            faceId_[count] = faceId;
-            count++;
+            // Boundary faces
+            if (!mesh_.isInternalFace(meshFacei))
+            {
+                facePatchId = mesh_.boundaryMesh().whichPatch(meshFacei);
+                const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
+
+                if (isA<emptyPolyPatch>(pp))
+                {
+                    continue;  // Ignore empty patch
+                }
+
+                const auto* cpp = isA<coupledPolyPatch>(pp);
+
+                if (cpp && !cpp->owner())
+                {
+                    continue;  // Ignore neighbour side
+                }
+
+                faceId = pp.whichFace(meshFacei);
+            }
+
+            if (faceId >= 0)
+            {
+                faceId_[numFaces] = faceId;
+                facePatchId_[numFaces] = facePatchId;
+
+                ++numFaces;
+            }
         }
     }
-    faceId_.setSize(count);
-    facePatchId_.setSize(count);
+
+    // Shrink to size used
+    faceId_.resize(numFaces);
+    facePatchId_.resize(numFaces);
 }
 
 
diff --git a/src/fvOptions/sources/derived/heatExchangerSource/heatExchangerModels/heatExchangerModel/heatExchangerModel.C b/src/fvOptions/sources/derived/heatExchangerSource/heatExchangerModels/heatExchangerModel/heatExchangerModel.C
index 9034a02b91..a9f3670093 100644
--- a/src/fvOptions/sources/derived/heatExchangerSource/heatExchangerModels/heatExchangerModel/heatExchangerModel.C
+++ b/src/fvOptions/sources/derived/heatExchangerSource/heatExchangerModels/heatExchangerModel/heatExchangerModel.C
@@ -78,60 +78,62 @@ void Foam::heatExchangerModel::initialise()
 
     const faceZone& fZone = mesh_.faceZones()[zoneID];
 
-    faceId_.setSize(fZone.size());
-    facePatchId_.setSize(fZone.size());
-    faceSign_.setSize(fZone.size());
+    // Total number of faces selected
+    label numFaces = fZone.size();
 
-    label count = 0;
-    forAll(fZone, i)
+    faceId_.resize_nocopy(numFaces);
+    facePatchId_.resize_nocopy(numFaces);
+    faceSign_.resize_nocopy(numFaces);
+
+    numFaces = 0;
+
+    // TDB: handle multiple zones
     {
-        const label facei = fZone[i];
-        label faceId = -1;
-        label facePatchId = -1;
-        if (mesh_.isInternalFace(facei))
+        forAll(fZone, i)
         {
-            faceId = facei;
-            facePatchId = -1;
-        }
-        else
-        {
-            facePatchId = mesh_.boundaryMesh().whichPatch(facei);
-            const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
-            const auto* cpp = isA<coupledPolyPatch>(pp);
+            const label meshFacei = fZone[i];
+            const label flipSign = (fZone.flipMap()[i] ? -1 : 1);
 
-            if (cpp)
-            {
-                faceId = (cpp->owner() ? pp.whichFace(facei) : -1);
-            }
-            else if (!isA<emptyPolyPatch>(pp))
-            {
-                faceId = pp.whichFace(facei);
-            }
-            else
-            {
-                faceId = -1;
-                facePatchId = -1;
-            }
-        }
+            // Internal faces
+            label faceId = meshFacei;
+            label facePatchId = -1;
 
-        if (faceId >= 0)
-        {
-            if (fZone.flipMap()[i])
+            // Boundary faces
+            if (!mesh_.isInternalFace(meshFacei))
             {
-                faceSign_[count] = -1;
+                facePatchId = mesh_.boundaryMesh().whichPatch(meshFacei);
+                const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
+
+                if (isA<emptyPolyPatch>(pp))
+                {
+                    continue;  // Ignore empty patch
+                }
+
+                const auto* cpp = isA<coupledPolyPatch>(pp);
+
+                if (cpp && !cpp->owner())
+                {
+                    continue;  // Ignore neighbour side
+                }
+
+                faceId = pp.whichFace(meshFacei);
             }
-            else
+
+            if (faceId >= 0)
             {
-                faceSign_[count] = 1;
+                faceId_[numFaces] = faceId;
+                facePatchId_[numFaces] = facePatchId;
+                faceSign_[numFaces] = flipSign;
+
+                ++numFaces;
             }
-            faceId_[count] = faceId;
-            facePatchId_[count] = facePatchId;
-            count++;
         }
     }
-    faceId_.setSize(count);
-    facePatchId_.setSize(count);
-    faceSign_.setSize(count);
+
+    // Shrink to size used
+    faceId_.resize(numFaces);
+    facePatchId_.resize(numFaces);
+    faceSign_.resize(numFaces);
 }
 
 
diff --git a/src/sampling/sampledSurface/sampledFaceZone/sampledFaceZone.C b/src/sampling/sampledSurface/sampledFaceZone/sampledFaceZone.C
index 9f5c6b167c..6ae6955171 100644
--- a/src/sampling/sampledSurface/sampledFaceZone/sampledFaceZone.C
+++ b/src/sampling/sampledSurface/sampledFaceZone/sampledFaceZone.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -151,19 +151,18 @@ bool Foam::sampledFaceZone::update()
     // Could also check numFaces
 
     // The mesh face or local patch face and the patch id
-    faceId_.resize(numFaces);
-    facePatchId_.resize(numFaces);
+    faceId_.resize_nocopy(numFaces);
+    facePatchId_.resize_nocopy(numFaces);
 
     IndirectList<face> selectedFaces(mesh().faces(), labelList());
     labelList& meshFaceIds = selectedFaces.addressing();
-    meshFaceIds.resize(numFaces);
+    meshFaceIds.resize_nocopy(numFaces);
 
     numFaces = 0;
 
-    forAll(zoneIDs(), idx)
+    for (const label zoneId : zoneIDs())
     {
-        const label zonei = zoneIDs()[idx];
-        const faceZone& fZone = mesh().faceZones()[zonei];
+        const faceZone& fZone = mesh().faceZones()[zoneId];
 
         for (const label meshFacei : fZone)
         {
@@ -179,26 +178,17 @@ bool Foam::sampledFaceZone::update()
 
                 if (isA<emptyPolyPatch>(pp))
                 {
-                    // Do not sample an empty patch
-                    continue;
-                }
-
-                const auto* procPatch = isA<processorPolyPatch>(pp);
-                if (procPatch && !procPatch->owner())
-                {
-                    // Do not sample neighbour-side, retain owner-side only
-                    continue;
+                    continue;  // Ignore empty patch
                 }
 
                 const auto* cpp = isA<coupledPolyPatch>(pp);
-                if (cpp)
+
+                if (cpp && !cpp->owner())
                 {
-                    faceId = (cpp->owner() ? pp.whichFace(meshFacei) : -1);
-                }
-                else
-                {
-                    faceId = pp.whichFace(meshFacei);
+                    continue;  // Ignore neighbour side
                 }
+
+                faceId = pp.whichFace(meshFacei);
             }
 
             if (faceId >= 0)
@@ -206,6 +196,7 @@ bool Foam::sampledFaceZone::update()
                 faceId_[numFaces] = faceId;
                 facePatchId_[numFaces] = facePatchId;
                 meshFaceIds[numFaces] = meshFacei;
+
                 ++numFaces;
             }
         }