diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
index f44176b09e..c66b642475 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -132,58 +132,76 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::initialisePatch()
     patchNormal_ /= mag(patchNormal_) + ROOTVSMALL;
 
 
-    // Decompose the patch faces into triangles to enable point search
-
     const polyPatch& patch = this->patch().patch();
     const pointField& points = patch.points();
 
     // Triangulate the patch faces and create addressing
-    DynamicList<label> triToFace(2*patch.size());
-    DynamicList<scalar> triMagSf(2*patch.size());
-    DynamicList<face> triFace(2*patch.size());
-    DynamicList<face> tris(5);
-
-    // Set zero value at the start of the tri area list
-    triMagSf.append(0.0);
-
-    forAll(patch, faceI)
     {
-        const face& f = patch[faceI];
-
-        tris.clear();
-        f.triangles(points, tris);
-
-        forAll(tris, i)
+        label nTris = 0;
+        for (const face& f : patch)
         {
-            triToFace.append(faceI);
-            triFace.append(tris[i]);
-            triMagSf.append(tris[i].mag(points));
+            nTris += f.nTriangles();
+        }
+
+        DynamicList<labelledTri> dynTriFace(nTris);
+        DynamicList<face> tris(8);  // work array
+
+        forAll(patch, facei)
+        {
+            const face& f = patch[facei];
+
+            tris.clear();
+            f.triangles(points, tris);
+
+            for (const auto& t : tris)
+            {
+                dynTriFace.emplace_back(t[0], t[1], t[2], facei);
+            }
+        }
+
+        // Transfer to persistent storage
+        triFace_.transfer(dynTriFace);
+    }
+
+
+    const label myProci = UPstream::myProcNo();
+    const label numProc = UPstream::nProcs();
+
+    // Calculate the cumulative triangle weights
+    {
+        triCumulativeMagSf_.resize_nocopy(triFace_.size()+1);
+
+        auto iter = triCumulativeMagSf_.begin();
+
+        // Set zero value at the start of the tri area/weight list
+        scalar patchArea = 0;
+        *iter++ = patchArea;
+
+        // Calculate cumulative and total area
+        for (const auto& t : triFace_)
+        {
+            patchArea += t.mag(points);
+            *iter++ = patchArea;
+        }
+
+        sumTriMagSf_.resize_nocopy(numProc+1);
+        sumTriMagSf_[0] = 0;
+
+        {
+            scalarList::subList slice(sumTriMagSf_, numProc, 1);
+            slice[myProci] = patchArea;
+            Pstream::allGatherList(slice);
+        }
+
+        // Convert to cumulative sum of areas per proc
+        for (label i = 1; i < sumTriMagSf_.size(); ++i)
+        {
+            sumTriMagSf_[i] += sumTriMagSf_[i-1];
         }
     }
 
-    sumTriMagSf_ = Zero;
-    sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
-
-    Pstream::listCombineReduce(sumTriMagSf_, maxEqOp<scalar>());
-
-    for (label i = 1; i < triMagSf.size(); ++i)
-    {
-        triMagSf[i] += triMagSf[i-1];
-    }
-
-    // Transfer to persistent storage
-    triFace_.transfer(triFace);
-    triToFace_.transfer(triToFace);
-    triCumulativeMagSf_.transfer(triMagSf);
-
-    // Convert sumTriMagSf_ into cumulative sum of areas per proc
-    for (label i = 1; i < sumTriMagSf_.size(); ++i)
-    {
-        sumTriMagSf_[i] += sumTriMagSf_[i-1];
-    }
-
     // Global patch area (over all processors)
-    patchArea_ = sumTriMagSf_.last();
+    patchArea_ = sumTriMagSf_.back();
 
     // Local patch bounds (this processor)
     patchBounds_ = boundBox(patch.localPoints(), false);
@@ -239,33 +257,32 @@ Foam::pointIndexHit Foam::turbulentDFSEMInletFvPatchVectorField::setNewPosition
     const bool global
 )
 {
-    // Initialise to miss
-    pointIndexHit pos(false, vector::max, -1);
-
     const polyPatch& patch = this->patch().patch();
     const pointField& points = patch.points();
 
+    label triI = -1;
+
     if (global)
     {
         const scalar areaFraction =
             rndGen_.globalPosition<scalar>(0, patchArea_);
 
         // Determine which processor to use
-        label procI = 0;
+        label proci = 0;
         forAllReverse(sumTriMagSf_, i)
         {
             if (areaFraction >= sumTriMagSf_[i])
             {
-                procI = i;
+                proci = i;
                 break;
             }
         }
 
-        if (Pstream::myProcNo() == procI)
+        if (UPstream::myProcNo() == proci)
         {
             // Find corresponding decomposed face triangle
-            label triI = 0;
-            const scalar offset = sumTriMagSf_[procI];
+            triI = 0;
+            const scalar offset = sumTriMagSf_[proci];
             forAllReverse(triCumulativeMagSf_, i)
             {
                 if (areaFraction > triCumulativeMagSf_[i] + offset)
@@ -274,20 +291,13 @@ Foam::pointIndexHit Foam::turbulentDFSEMInletFvPatchVectorField::setNewPosition
                     break;
                 }
             }
-
-            // Find random point in triangle
-            const face& tf = triFace_[triI];
-            const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
-
-            pos.hitPoint(tri.randomPoint(rndGen_));
-            pos.setIndex(triToFace_[triI]);
         }
     }
     else
     {
         // Find corresponding decomposed face triangle on local processor
-        label triI = 0;
-        const scalar maxAreaLimit = triCumulativeMagSf_.last();
+        triI = 0;
+        const scalar maxAreaLimit = triCumulativeMagSf_.back();
         const scalar areaFraction = rndGen_.position<scalar>(0, maxAreaLimit);
 
         forAllReverse(triCumulativeMagSf_, i)
@@ -298,16 +308,22 @@ Foam::pointIndexHit Foam::turbulentDFSEMInletFvPatchVectorField::setNewPosition
                 break;
             }
         }
-
-        // Find random point in triangle
-        const face& tf = triFace_[triI];
-        const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
-
-        pos.hitPoint(tri.randomPoint(rndGen_));
-        pos.setIndex(triToFace_[triI]);
     }
 
-    return pos;
+
+    if (triI >= 0)
+    {
+        return pointIndexHit
+        (
+            true,
+            // Find random point in triangle
+            triFace_[triI].tri(points).randomPoint(rndGen_),
+            triFace_[triI].index()
+        );
+    }
+
+    // No hit
+    return pointIndexHit(false, vector::max, -1);
 }
 
 
@@ -573,7 +589,6 @@ turbulentDFSEMInletFvPatchVectorField
 
     patchArea_(-1),
     triFace_(),
-    triToFace_(),
     triCumulativeMagSf_(),
     sumTriMagSf_(Pstream::nProcs() + 1, Zero),
     patchNormal_(Zero),
@@ -615,7 +630,6 @@ turbulentDFSEMInletFvPatchVectorField
 
     patchArea_(ptf.patchArea_),
     triFace_(ptf.triFace_),
-    triToFace_(ptf.triToFace_),
     triCumulativeMagSf_(ptf.triCumulativeMagSf_),
     sumTriMagSf_(ptf.sumTriMagSf_),
     patchNormal_(ptf.patchNormal_),
@@ -656,7 +670,6 @@ turbulentDFSEMInletFvPatchVectorField
 
     patchArea_(-1),
     triFace_(),
-    triToFace_(),
     triCumulativeMagSf_(),
     sumTriMagSf_(Pstream::nProcs() + 1, Zero),
     patchNormal_(Zero),
@@ -684,10 +697,11 @@ turbulentDFSEMInletFvPatchVectorField
 Foam::turbulentDFSEMInletFvPatchVectorField::
 turbulentDFSEMInletFvPatchVectorField
 (
-    const turbulentDFSEMInletFvPatchVectorField& ptf
+    const turbulentDFSEMInletFvPatchVectorField& ptf,
+    const DimensionedField<vector, volMesh>& iF
 )
 :
-    fixedValueFvPatchField<vector>(ptf),
+    fixedValueFvPatchField<vector>(ptf, iF),
     U_(ptf.U_.clone(patch().patch())),
     R_(ptf.R_.clone(patch().patch())),
     L_(ptf.L_.clone(patch().patch())),
@@ -702,7 +716,6 @@ turbulentDFSEMInletFvPatchVectorField
 
     patchArea_(ptf.patchArea_),
     triFace_(ptf.triFace_),
-    triToFace_(ptf.triToFace_),
     triCumulativeMagSf_(ptf.triCumulativeMagSf_),
     sumTriMagSf_(ptf.sumTriMagSf_),
     patchNormal_(ptf.patchNormal_),
@@ -723,40 +736,10 @@ turbulentDFSEMInletFvPatchVectorField
 Foam::turbulentDFSEMInletFvPatchVectorField::
 turbulentDFSEMInletFvPatchVectorField
 (
-    const turbulentDFSEMInletFvPatchVectorField& ptf,
-    const DimensionedField<vector, volMesh>& iF
+    const turbulentDFSEMInletFvPatchVectorField& ptf
 )
 :
-    fixedValueFvPatchField<vector>(ptf, iF),
-    U_(ptf.U_.clone(patch().patch())),
-    R_(ptf.R_.clone(patch().patch())),
-    L_(ptf.L_.clone(patch().patch())),
-    delta_(ptf.delta_),
-    d_(ptf.d_),
-    kappa_(ptf.kappa_),
-    Uref_(ptf.Uref_),
-    Lref_(ptf.Lref_),
-    scale_(ptf.scale_),
-    m_(ptf.m_),
-    nCellPerEddy_(ptf.nCellPerEddy_),
-
-    patchArea_(ptf.patchArea_),
-    triFace_(ptf.triFace_),
-    triToFace_(ptf.triToFace_),
-    triCumulativeMagSf_(ptf.triCumulativeMagSf_),
-    sumTriMagSf_(ptf.sumTriMagSf_),
-    patchNormal_(ptf.patchNormal_),
-    patchBounds_(ptf.patchBounds_),
-
-    eddies_(ptf.eddies_),
-    v0_(ptf.v0_),
-    rndGen_(ptf.rndGen_),
-    sigmax_(ptf.sigmax_),
-    maxSigmaX_(ptf.maxSigmaX_),
-    nEddy_(ptf.nEddy_),
-    curTimeIndex_(-1),
-    singleProc_(ptf.singleProc_),
-    writeEddies_(ptf.writeEddies_)
+    turbulentDFSEMInletFvPatchVectorField(ptf, ptf.internalField())
 {}
 
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H
index 862684eb1e..6755c1c5b6 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -146,6 +146,7 @@ SourceFiles
 #include "eddy.H"
 #include "pointIndexHit.H"
 #include "PatchFunction1.H"
+#include "labelledTri.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -204,13 +205,11 @@ class turbulentDFSEMInletFvPatchVectorField
             //- Patch area - total across all processors
             scalar patchArea_;
 
-            //- Decomposed patch faces as a list of triangles
-            faceList triFace_;
+            //- The polyPatch faces as triangles, the index of each corresponds
+            //- to the undecomposed patch face index.
+            List<labelledTri> triFace_;
 
-            //- Addressing from per triangle to patch face
-            labelList triToFace_;
-
-            //- Cumulative triangle area per triangle face
+            //- Cumulative triangle area per triangle face (processor-local)
             scalarList triCumulativeMagSf_;
 
             //- Cumulative area fractions per processor
diff --git a/src/lagrangian/basic/Cloud/CloudIO.C b/src/lagrangian/basic/Cloud/CloudIO.C
index d804f755c7..3cb8569df6 100644
--- a/src/lagrangian/basic/Cloud/CloudIO.C
+++ b/src/lagrangian/basic/Cloud/CloudIO.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017, 2020 OpenFOAM Foundation
-    Copyright (C) 2017-2023 OpenCFD Ltd.
+    Copyright (C) 2017-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -101,11 +101,9 @@ void Foam::Cloud<ParticleType>::writeCloudUniformProperties() const
         )
     );
 
-    labelList np(Pstream::nProcs(), Zero);
-    np[Pstream::myProcNo()] = ParticleType::particleCount_;
-
-    // FIXME: replace with  Pstream::allGatherList(np);
-    Pstream::listCombineReduce(np, maxEqOp<label>());
+    labelList np(UPstream::nProcs(), Foam::zero{});
+    np[UPstream::myProcNo()] = ParticleType::particleCount_;
+    Pstream::allGatherList(np);
 
     uniformPropsDict.add
     (
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C
index 36b5a7d874..80b34fb203 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2017 OpenFOAM Foundation
+    Copyright (C) 2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,13 +44,12 @@ Foam::patchInjectionBase::patchInjectionBase
 :
     patchName_(patchName),
     patchId_(mesh.boundaryMesh().findPatchID(patchName_)),
-    patchArea_(0.0),
+    patchArea_(0),
     patchNormal_(),
     cellOwners_(),
     triFace_(),
-    triToFace_(),
     triCumulativeMagSf_(),
-    sumTriMagSf_(Pstream::nProcs() + 1, Zero)
+    sumTriMagSf_()
 {
     if (patchId_ < 0)
     {
@@ -71,7 +71,6 @@ Foam::patchInjectionBase::patchInjectionBase(const patchInjectionBase& pib)
     patchNormal_(pib.patchNormal_),
     cellOwners_(pib.cellOwners_),
     triFace_(pib.triFace_),
-    triToFace_(pib.triToFace_),
     triCumulativeMagSf_(pib.triCumulativeMagSf_),
     sumTriMagSf_(pib.sumTriMagSf_)
 {}
@@ -88,48 +87,68 @@ void Foam::patchInjectionBase::updateMesh(const polyMesh& mesh)
     cellOwners_ = patch.faceCells();
 
     // Triangulate the patch faces and create addressing
-    DynamicList<label> triToFace(2*patch.size());
-    DynamicList<scalar> triMagSf(2*patch.size());
-    DynamicList<face> triFace(2*patch.size());
-    DynamicList<face> tris(5);
-
-    // Set zero value at the start of the tri area list
-    triMagSf.append(0.0);
-
-    forAll(patch, facei)
     {
-        const face& f = patch[facei];
-
-        tris.clear();
-        f.triangles(points, tris);
-
-        forAll(tris, i)
+        label nTris = 0;
+        for (const face& f : patch)
         {
-            triToFace.append(facei);
-            triFace.append(tris[i]);
-            triMagSf.append(tris[i].mag(points));
+            nTris += f.nTriangles();
         }
+
+        DynamicList<labelledTri> dynTriFace(nTris);
+        DynamicList<face> tris(8);  // work array
+
+        forAll(patch, facei)
+        {
+            const face& f = patch[facei];
+
+            tris.clear();
+            f.triangles(points, tris);
+
+            for (const auto& t : tris)
+            {
+                dynTriFace.emplace_back(t[0], t[1], t[2], facei);
+            }
+        }
+
+        // Transfer to persistent storage
+        triFace_.transfer(dynTriFace);
     }
 
-    sumTriMagSf_ = Zero;
-    sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
 
-    Pstream::listCombineReduce(sumTriMagSf_, maxEqOp<scalar>());
+    const label myProci = UPstream::myProcNo();
+    const label numProc = UPstream::nProcs();
 
-    for (label i = 1; i < triMagSf.size(); i++)
+    // Calculate the cumulative triangle weights
     {
-        triMagSf[i] += triMagSf[i-1];
-    }
+        triCumulativeMagSf_.resize_nocopy(triFace_.size()+1);
 
-    // Transfer to persistent storage
-    triFace_.transfer(triFace);
-    triToFace_.transfer(triToFace);
-    triCumulativeMagSf_.transfer(triMagSf);
+        auto iter = triCumulativeMagSf_.begin();
 
-    // Convert sumTriMagSf_ into cumulative sum of areas per proc
-    for (label i = 1; i < sumTriMagSf_.size(); i++)
-    {
-        sumTriMagSf_[i] += sumTriMagSf_[i-1];
+        // Set zero value at the start of the tri area/weight list
+        scalar patchArea = 0;
+        *iter++ = patchArea;
+
+        // Calculate cumulative and total area
+        for (const auto& t : triFace_)
+        {
+            patchArea += t.mag(points);
+            *iter++ = patchArea;
+        }
+
+        sumTriMagSf_.resize_nocopy(numProc+1);
+        sumTriMagSf_[0] = 0;
+
+        {
+            scalarList::subList slice(sumTriMagSf_, numProc, 1);
+            slice[myProci] = patchArea;
+            Pstream::allGatherList(slice);
+        }
+
+        // Convert to cumulative
+        for (label i = 1; i < sumTriMagSf_.size(); ++i)
+        {
+            sumTriMagSf_[i] += sumTriMagSf_[i-1];
+        }
     }
 
     const scalarField magSf(mag(patch.faceAreas()));
@@ -152,12 +171,12 @@ Foam::label Foam::patchInjectionBase::setPositionAndCell
 {
     label facei = -1;
 
-    if (cellOwners_.size() > 0)
+    if (!cellOwners_.empty())
     {
         // Determine which processor to inject from
         const label proci = whichProc(fraction01);
 
-        if (Pstream::myProcNo() == proci)
+        if (UPstream::myProcNo() == proci)
         {
             const scalar areaFraction = fraction01*patchArea_;
 
@@ -174,15 +193,13 @@ Foam::label Foam::patchInjectionBase::setPositionAndCell
             }
 
             // Set cellOwner
-            facei = triToFace_[trii];
+            facei = triFace_[trii].index();
             cellOwner = cellOwners_[facei];
 
             // Find random point in triangle
             const polyPatch& patch = mesh.boundaryMesh()[patchId_];
             const pointField& points = patch.points();
-            const face& tf = triFace_[trii];
-            const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
-            const point pf(tri.randomPoint(rnd));
+            const point pf = triFace_[trii].tri(points).randomPoint(rnd);
 
             // Position perturbed away from face (into domain)
             const scalar a = rnd.position(scalar(0.1), scalar(0.5));
@@ -239,17 +256,9 @@ Foam::label Foam::patchInjectionBase::setPositionAndCell
                 tetPti = cellTetIs[teti].tetPt();
             }
         }
-        else
-        {
-            cellOwner = -1;
-            tetFacei = -1;
-            tetPti = -1;
-
-            // Dummy position
-            position = pTraits<vector>::max;
-        }
     }
-    else
+
+    if (facei == -1)
     {
         cellOwner = -1;
         tetFacei = -1;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.H
index 6c2336670a..c7bb7f2ec4 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,21 +40,22 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef patchInjectionBase_H
-#define patchInjectionBase_H
+#ifndef Foam_patchInjectionBase_H
+#define Foam_patchInjectionBase_H
 
 #include "word.H"
 #include "labelList.H"
 #include "scalarList.H"
 #include "vectorList.H"
 #include "faceList.H"
+#include "labelledTri.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// Forward declarations
+// Forward Declarations
 class polyMesh;
 class fvMesh;
 class Random;
@@ -66,7 +68,7 @@ class patchInjectionBase
 {
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Patch name
         const word patchName_;
@@ -83,13 +85,11 @@ protected:
         //- List of cell labels corresponding to injector positions
         labelList cellOwners_;
 
-        //- Decomposed patch faces as a list of triangles
-        faceList triFace_;
+        //- The polyPatch faces as triangles, the index of each corresponds
+        //- to the undecomposed patch face index.
+        List<labelledTri> triFace_;
 
-        //- Addressing from per triangle to patch face
-        labelList triToFace_;
-
-        //- Cumulative triangle area per triangle face
+        //- Cumulative triangle area per triangle face (processor-local)
         scalarList triCumulativeMagSf_;
 
         //- Cumulative area fractions per processor
diff --git a/src/meshTools/triangulatedPatch/triangulatedPatch.C b/src/meshTools/triangulatedPatch/triangulatedPatch.C
index f80ec25595..ab2763b9af 100644
--- a/src/meshTools/triangulatedPatch/triangulatedPatch.C
+++ b/src/meshTools/triangulatedPatch/triangulatedPatch.C
@@ -171,31 +171,32 @@ void Foam::triangulatedPatch::update()
 
     // Set zero value at the start of the tri area/weight list
     scalar patchArea = 0;
-    *iter = patchArea;
-    ++iter;
+    *iter++ = patchArea;
 
     // Calculate cumulative and total area (processor-local at this point)
     for (const auto& t : triFace_)
     {
         patchArea += t.mag(points);
-
-        *iter = patchArea;
-        ++iter;
+        *iter++ = patchArea;
     }
 
-    // FIXME: use allGatherList of subslice
-    scalarList procSumWght(numProc+1, Foam::zero{});
-    procSumWght[myProci+1] = patchArea;
-    Pstream::listCombineReduce(procSumWght, maxEqOp<scalar>());
+    scalarList procSumArea(numProc+1);
+    procSumArea[0] = 0;
+
+    {
+        scalarList::subList slice(procSumArea, numProc, 1);
+        slice[myProci] = patchArea;
+        Pstream::allGatherList(slice);
+    }
 
     // Convert to cumulative
-    for (label i = 1; i < procSumWght.size(); ++i)
+    for (label i = 1; i < procSumArea.size(); ++i)
     {
-        procSumWght[i] += procSumWght[i-1];
+        procSumArea[i] += procSumArea[i-1];
     }
 
-    const scalar offset = procSumWght[myProci];
-    const scalar totalArea = procSumWght.back();
+    const scalar offset = procSumArea[myProci];
+    const scalar totalArea = procSumArea.back();
 
     // Apply processor offset and normalise - for a global 0-1 interval
     for (scalar& w : triWght_)
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.C
index 1fc071bbd7..bc33b041e5 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/topOVariablesBase/topOVariablesBase.C
@@ -669,19 +669,13 @@ void Foam::topOVariablesBase::writeFluidSolidInterface
     //    invertOneToMany(surfFaces_.size(), changedFaceToCutFace);
 
     // Transform origin cut faces to a global numbering
-    labelList cuttingFacesPerProc(Pstream::nProcs(), Zero);
-    cuttingFacesPerProc[Pstream::myProcNo()] = changedFaces.size();
-    Pstream::listCombineReduce(cuttingFacesPerProc, plusEqOp<label>());
 
-    labelList passedFaces(Pstream::nProcs(), Zero);
-    for (label i = 1; i < Pstream::nProcs(); ++i)
-    {
-        passedFaces[i] = passedFaces[i - 1] + cuttingFacesPerProc[i - 1];
-    }
+    const label changedFacesOffset =
+        globalIndex::calcOffset(changedFaces.size());
 
-    forAll(changedFacesInfo, facei)
+    for (auto& info : changedFacesInfo)
     {
-        changedFacesInfo[facei].data() += passedFaces[Pstream::myProcNo()];
+        info.data() += changedFacesOffset;
     }
 }