diff --git a/applications/test/checkDecomposePar/Test-checkDecomposePar.C b/applications/test/checkDecomposePar/Test-checkDecomposePar.C
index 8019fd0d45..903a22b0d4 100644
--- a/applications/test/checkDecomposePar/Test-checkDecomposePar.C
+++ b/applications/test/checkDecomposePar/Test-checkDecomposePar.C
@@ -123,15 +123,9 @@ int main(int argc, char *argv[])
 
         const label nDomains = max(cellToProc) + 1;
 
+        // Local mesh connectivity
         CompactListList<label> cellCells;
-        globalMeshData::calcCellCells
-        (
-            mesh,
-            identity(mesh.nCells()),
-            mesh.nCells(),
-            false,
-            cellCells
-        );
+        globalMeshData::calcCellCells(mesh, cellCells);
 
         decompositionInformation info
         (
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.C b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.C
index 36d04152d7..af5cfeca98 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.C
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -112,15 +112,9 @@ void Foam::domainDecompositionDryRun::execute
 
     decompositionMethod& method = model.decomposer();
 
+    // Local mesh connectivity
     CompactListList<label> cellCells;
-    globalMeshData::calcCellCells
-    (
-        mesh_,
-        identity(mesh_.nCells()),
-        mesh_.nCells(),
-        false,  // false = local only
-        cellCells
-    );
+    globalMeshData::calcCellCells(mesh_, cellCells);
 
     labelList cellToProc = method.decompose(mesh_, cellWeights);
 
diff --git a/src/OpenFOAM/meshes/bandCompression/bandCompression.C b/src/OpenFOAM/meshes/bandCompression/bandCompression.C
index ffcb574f48..b2f1a977ad 100644
--- a/src/OpenFOAM/meshes/bandCompression/bandCompression.C
+++ b/src/OpenFOAM/meshes/bandCompression/bandCompression.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2022-2023 OpenCFD Ltd.
+    Copyright (C) 2022-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +28,8 @@ License
 
 #include "bandCompression.H"
 #include "bitSet.H"
+#include "polyMesh.H"
+#include "globalMeshData.H"
 #include "CircularBuffer.H"
 #include "CompactListList.H"
 #include "DynamicList.H"
@@ -173,10 +175,10 @@ Foam::labelList Foam::meshTools::bandCompression
 )
 {
     // Protect against zero-sized offset list
-    const label nOldCells = max(0, (offsets.size()-1));
+    const label nOldCells = Foam::max(0, (offsets.size()-1));
 
     // Count number of neighbours
-    labelList numNbrs(nOldCells, Zero);
+    labelList numNbrs(nOldCells, Foam::zero{});
     for (label celli = 0; celli < nOldCells; ++celli)
     {
         const label beg = offsets[celli];
@@ -312,6 +314,16 @@ Foam::labelList Foam::meshTools::bandCompression
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+Foam::labelList Foam::meshTools::bandCompression(const polyMesh& mesh)
+{
+    // Local mesh connectivity
+    CompactListList<label> cellCells;
+    globalMeshData::calcCellCells(mesh, cellCells);
+
+    return cuthill_mckee_algorithm(cellCells);
+}
+
+
 Foam::labelList Foam::meshTools::bandCompression
 (
     const CompactListList<label>& cellCellAddressing
diff --git a/src/OpenFOAM/meshes/bandCompression/bandCompression.H b/src/OpenFOAM/meshes/bandCompression/bandCompression.H
index 805a0783bf..d90cca8e58 100644
--- a/src/OpenFOAM/meshes/bandCompression/bandCompression.H
+++ b/src/OpenFOAM/meshes/bandCompression/bandCompression.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2022-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,6 +50,7 @@ namespace Foam
 {
 
 // Forward Declarations
+class polyMesh;
 template<class T> class CompactListList;
 
 } // End namespace Foam
@@ -64,6 +65,12 @@ namespace meshTools
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+//- Renumber (mesh) addressing to reduce the band of the mesh connectivity,
+//- using the Cuthill-McKee algorithm.
+//
+//  \returns order in which the cells are to be visited (ordered to original)
+labelList bandCompression(const polyMesh& mesh);
+
 //- Renumber (mesh) addressing to reduce the band of the matrix,
 //- using the Cuthill-McKee algorithm.
 //
@@ -97,12 +104,14 @@ namespace Foam
 {
 
 //- Forward to meshTools::bandCompression
+//  \deprecated(2022-03) prefer meshTools::bandCompression()
 inline labelList bandCompression(const labelListList& cellCellAddressing)
 {
     return meshTools::bandCompression(cellCellAddressing);
 }
 
 //- Forward to meshTools::bandCompression
+//  \deprecated(2022-03) prefer meshTools::bandCompression()
 inline labelList bandCompression
 (
     const labelUList& cellCells,
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
index a1228f6545..a06582b94e 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -614,6 +614,7 @@ public:
 
         //- Determine (local or global) cellCells from mesh agglomeration.
         //  Agglomeration is local to the processor.
+        //
         //  - parallel = false
         //    Resulting connections are in local cell indices.
         //    Coupled across cyclics but not processor patches.
@@ -623,9 +624,14 @@ public:
         static void calcCellCells
         (
             const polyMesh& mesh,
-            const labelList& agglom,
+            //! The cell agglomeration, negative agglomeration are excluded
+            //! from the sub-mesh
+            const labelUList& agglom,
+            //! The number of unique, agglomerated cells
             const label nLocalCoarse,
-            const bool parallel,  //!< Use global cell ids in parallel
+            //! Use global cell ids in parallel
+            const bool parallel,
+            //! [out] the mesh connectivity as CSR
             CompactListList<label>& cellCells
         );
 
@@ -635,12 +641,77 @@ public:
         static void calcCellCells
         (
             const polyMesh& mesh,
-            const labelList& agglom,
+            //! The cell agglomeration, negative agglomeration are excluded
+            //! from the sub-mesh
+            const labelUList& agglom,
+            //! The number of unique, agglomerated cells
             const label nLocalCoarse,
-            const bool parallel,  //!< Use global cell ids in parallel
+            //! Use global cell ids in parallel
+            const bool parallel,
+            //! [out] the mesh connectivity as CSR
             CompactListList<label>& cellCells,
+            //! [out] the connectivity weights (face area)
             CompactListList<scalar>& cellCellWeights
         );
+
+        //- Determine (local or global) mesh connectivity
+        //
+        //  - parallel = false
+        //    Resulting connections are in local cell indices.
+        //    Coupled across cyclics but not processor patches.
+        //  - parallel = true
+        //    Resulting connections are in global cell indices.
+        //    Coupled across cyclics and processor patches.
+        static void calcCellCells
+        (
+            const polyMesh& mesh,
+            //! [out] the mesh connectivity as CSR
+            CompactListList<label>& cellCells,
+            //! Use global cell ids in parallel
+            const bool parallel = false
+        );
+
+        //- Determine (local or global) sub-mesh connectivity
+        //
+        //  - parallel = false
+        //    Resulting connections are in local cell indices.
+        //    Coupled across cyclics but not processor patches.
+        //  - parallel = true
+        //    Resulting connections are in global cell indices.
+        //    Coupled across cyclics and processor patches.
+        //
+        //  \return the cellMap
+        static labelList calcCellCells
+        (
+            const polyMesh& mesh,
+            //! The cell-subset
+            const bitSet& selectedCells,
+            //! [out] the mesh connectivity as CSR
+            CompactListList<label>& cellCells,
+            //! Use global cell ids in parallel
+            const bool parallel = false
+        );
+
+        //- Determine (local or global) sub-mesh connectivity
+        //
+        //  - parallel = false
+        //    Resulting connections are in local cell indices.
+        //    Coupled across cyclics but not processor patches.
+        //  - parallel = true
+        //    Resulting connections are in global cell indices.
+        //    Coupled across cyclics and processor patches.
+        //
+        //  \return the cellMap
+        static labelList calcCellCells
+        (
+            const polyMesh& mesh,
+            //! The cell-subset
+            const labelUList& selectedCells,
+            //! [out] the mesh connectivity as CSR
+            CompactListList<label>& cellCells,
+            //! Use global cell ids in parallel
+            const bool parallel = false
+        );
 };
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshDataTopology.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshDataTopology.C
index 7743b75714..3a26d28681 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshDataTopology.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshDataTopology.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2015-2023 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,33 +32,58 @@ License
 #include "processorPolyPatch.H"
 #include "syncTools.H"
 
-// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
-void Foam::globalMeshData::calcCellCells
+namespace Foam
+{
+
+// NOTE: the AgglomerationType is anything that behaves like a List with
+// an operator[] and provides coverage in the (0-nCells) range.
+// - identityOp() does this
+
+template<class AgglomerationType>
+static void calcCellCellsImpl
 (
     const polyMesh& mesh,
-    const labelList& agglom,
+    const AgglomerationType& agglom,
     const label nLocalCoarse,
     const bool parallel,
-    CompactListList<label>& cellCells
+    CompactListList<label>& cellCells,
+    CompactListList<scalar>* cellCellWeightsPtr = nullptr
 )
 {
     const labelList& faceOwner = mesh.faceOwner();
     const labelList& faceNeigh = mesh.faceNeighbour();
     const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
-    // FUTURE? treat empty agglomeration like an identity map
-
     // Global cell numbers (agglomerated numbering)
     const label myProci = UPstream::myProcNo(UPstream::worldComm);
     const globalIndex globalAgglom(nLocalCoarse, UPstream::worldComm, parallel);
 
-
     // The agglomerated owner per boundary faces (global numbering)
     // from the other side (of coupled patches)
 
-    labelList globalNeighbour(agglom, pbm.faceOwner());
-    globalAgglom.inplaceToGlobal(myProci, globalNeighbour);
+    labelList globalNeighbour;
+    {
+        const label myAgglomOffset = globalAgglom.localStart(myProci);
+
+        const labelList::subList bndFaceOwner = pbm.faceOwner();
+
+        const label nBoundaryFaces = bndFaceOwner.size();
+
+        globalNeighbour.resize(nBoundaryFaces);
+
+        for (label bfacei = 0; bfacei < nBoundaryFaces; ++bfacei)
+        {
+            label val = agglom[bndFaceOwner[bfacei]];
+            if (val >= 0)
+            {
+                // Only offset 'real' (non-negative) agglomerations
+                val += myAgglomOffset;
+            }
+            globalNeighbour[bfacei] = val;
+        }
+    }
     syncTools::swapBoundaryFaceList(mesh, globalNeighbour);
 
 
@@ -66,13 +91,16 @@ void Foam::globalMeshData::calcCellCells
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Number of faces per coarse cell
-    labelList nFacesPerCell(nLocalCoarse, Zero);
+    labelList nFacesPerCell(nLocalCoarse, Foam::zero{});
 
     for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
     {
         const label own = agglom[faceOwner[facei]];
         const label nei = agglom[faceNeigh[facei]];
 
+        // Negative agglomeration (exclude from subset)
+        if (own < 0 || nei < 0) continue;
+
         ++nFacesPerCell[own];
         ++nFacesPerCell[nei];
     }
@@ -81,12 +109,15 @@ void Foam::globalMeshData::calcCellCells
     {
         if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
         {
-            label bFacei = pp.start()-mesh.nInternalFaces();
+            const label bndOffset = mesh.nInternalFaces();
 
             for (const label facei : pp.range())
             {
                 const label own = agglom[faceOwner[facei]];
-                const label globalNei = globalNeighbour[bFacei];
+                const label globalNei = globalNeighbour[facei-bndOffset];
+
+                // Negative agglomeration (exclude from subset)
+                if (own < 0 || globalNei < 0) continue;
 
                 if
                 (
@@ -96,8 +127,6 @@ void Foam::globalMeshData::calcCellCells
                 {
                     ++nFacesPerCell[own];
                 }
-
-                ++bFacei;
             }
         }
     }
@@ -107,23 +136,45 @@ void Foam::globalMeshData::calcCellCells
     // ~~~~~~~~~~~~~~~~~~~~~~~
 
     cellCells.resize_nocopy(nFacesPerCell);
-
     nFacesPerCell = 0;  // Restart the count
 
-    labelList& m = cellCells.values();
+
+    // CSR indexing
     const labelList& offsets = cellCells.offsets();
 
+    // CSR connections
+    labelList& connect = cellCells.values();
+
+
+    // CSR connection weights
+    scalarList weights;
+    if (cellCellWeightsPtr)
+    {
+        cellCellWeightsPtr->clear();
+        weights.resize(cellCells.totalSize());
+    }
+
+
     // For internal faces is just offsetted owner and neighbour
     for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
     {
         const label own = agglom[faceOwner[facei]];
         const label nei = agglom[faceNeigh[facei]];
 
+        // Negative agglomeration (exclude from subset)
+        if (own < 0 || nei < 0) continue;
+
         const label ownIndex = offsets[own] + nFacesPerCell[own]++;
         const label neiIndex = offsets[nei] + nFacesPerCell[nei]++;
 
-        m[ownIndex] = globalAgglom.toGlobal(myProci, nei);
-        m[neiIndex] = globalAgglom.toGlobal(myProci, own);
+        connect[ownIndex] = globalAgglom.toGlobal(myProci, nei);
+        connect[neiIndex] = globalAgglom.toGlobal(myProci, own);
+
+        if (!weights.empty())
+        {
+            weights[ownIndex] = Foam::mag(mesh.faceAreas()[facei]);
+            weights[neiIndex] = weights[ownIndex];
+        }
     }
 
     // For boundary faces is offsetted coupled neighbour
@@ -131,12 +182,15 @@ void Foam::globalMeshData::calcCellCells
     {
         if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
         {
-            label bFacei = pp.start()-mesh.nInternalFaces();
+            const label bndOffset = mesh.nInternalFaces();
 
             for (const label facei : pp.range())
             {
                 const label own = agglom[faceOwner[facei]];
-                const label globalNei = globalNeighbour[bFacei];
+                const label globalNei = globalNeighbour[facei-bndOffset];
+
+                // Negative agglomeration (exclude from subset)
+                if (own < 0 || globalNei < 0) continue;
 
                 if
                 (
@@ -146,10 +200,13 @@ void Foam::globalMeshData::calcCellCells
                 {
                     const label ownIndex = offsets[own] + nFacesPerCell[own]++;
 
-                    m[ownIndex] = globalNei;
-                }
+                    connect[ownIndex] = globalNei;
 
-                ++bFacei;
+                    if (!weights.empty())
+                    {
+                        weights[ownIndex] = Foam::mag(mesh.faceAreas()[facei]);
+                    }
+                }
             }
         }
     }
@@ -164,13 +221,12 @@ void Foam::globalMeshData::calcCellCells
 
     if (!cellCells.empty())
     {
-        label newIndex = 0;
-        labelHashSet nbrCells;
-
-        labelList& m = cellCells.values();
+        // Need non-const access to CSR indexing
         labelList& offsets = cellCells.offsets();
 
         label startIndex = offsets[0];
+        label newIndex = 0;
+        labelHashSet nbrCells;
 
         const label nCellCells = cellCells.size();
 
@@ -189,9 +245,15 @@ void Foam::globalMeshData::calcCellCells
 
             for (label i = startIndex; i < endIndex; ++i)
             {
-                if (nbrCells.insert(m[i]))
+                if (nbrCells.insert(connect[i]))
                 {
-                    m[newIndex] = m[i];
+                    connect[newIndex] = connect[i];
+
+                    if (!weights.empty())
+                    {
+                        weights[newIndex] = weights[i];
+                    }
+
                     ++newIndex;
                 }
             }
@@ -199,9 +261,23 @@ void Foam::globalMeshData::calcCellCells
             offsets[celli+1] = newIndex;
         }
 
-        m.resize(newIndex);
+        connect.resize(newIndex);
+        if (!weights.empty())
+        {
+            weights.resize(newIndex);
+        }
     }
 
+
+    // CSR connection weights
+    // - addressing is identical to the connections
+    if (cellCellWeightsPtr)
+    {
+        cellCellWeightsPtr->offsets() = cellCells.offsets();
+        cellCellWeightsPtr->values() = std::move(weights);
+    }
+
+
     //forAll(cellCells, celli)
     //{
     //    Pout<< "Original: Coarse cell " << celli << endl;
@@ -218,183 +294,194 @@ void Foam::globalMeshData::calcCellCells
     //}
 }
 
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
 void Foam::globalMeshData::calcCellCells
 (
     const polyMesh& mesh,
-    const labelList& agglom,
+    const labelUList& agglom,
+    const label nLocalCoarse,
+    const bool parallel,
+    CompactListList<label>& cellCells
+)
+{
+    calcCellCellsImpl
+    (
+        mesh,
+        agglom,
+        nLocalCoarse,
+        parallel,
+        cellCells
+    );
+}
+
+
+void Foam::globalMeshData::calcCellCells
+(
+    const polyMesh& mesh,
+    const labelUList& agglom,
     const label nLocalCoarse,
     const bool parallel,
     CompactListList<label>& cellCells,
     CompactListList<scalar>& cellCellWeights
 )
 {
-    const labelList& faceOwner = mesh.faceOwner();
-    const labelList& faceNeigh = mesh.faceNeighbour();
-    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
-
-    // FUTURE? treat empty agglomeration like an identity map
-
-    // Global cell numbers (agglomerated numbering)
-    const label myProci = UPstream::myProcNo(UPstream::worldComm);
-    const globalIndex globalAgglom(nLocalCoarse, UPstream::worldComm, parallel);
+    calcCellCellsImpl
+    (
+        mesh,
+        agglom,
+        nLocalCoarse,
+        parallel,
+        cellCells,
+       &cellCellWeights
+    );
+}
 
 
-    // The agglomerated owner per boundary faces (global numbering)
-    // from the other side (of coupled patches)
+// Convenience forms
 
-    labelList globalNeighbour(agglom, pbm.faceOwner());
-    globalAgglom.inplaceToGlobal(myProci, globalNeighbour);
-    syncTools::swapBoundaryFaceList(mesh, globalNeighbour);
+void Foam::globalMeshData::calcCellCells
+(
+    const polyMesh& mesh,
+    CompactListList<label>& cellCells,
+    const bool parallel
+)
+{
+    calcCellCellsImpl
+    (
+        mesh,
+        Foam::identityOp{},
+        mesh.nCells(),
+        parallel,
+        cellCells
+    );
+}
 
 
-    // Count number of faces (internal + coupled)
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Foam::labelList Foam::globalMeshData::calcCellCells
+(
+    const polyMesh& mesh,
+    const bitSet& selectedCells,
+    CompactListList<label>& cellCells,
+    const bool parallel
+)
+{
+    const label nCells = mesh.nCells();
 
-    // Number of faces per coarse cell
-    labelList nFacesPerCell(nLocalCoarse, Zero);
+    labelList agglom(nCells, -1);
+    labelList cellMap;
 
-    for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
+    // First pass - sorted order without duplicates
+    label nCompact = 0;
+    for (const label celli : selectedCells)
     {
-        const label own = agglom[faceOwner[facei]];
-        const label nei = agglom[faceNeigh[facei]];
-
-        ++nFacesPerCell[own];
-        ++nFacesPerCell[nei];
+        // A bitSet has no negatives/duplicates, so just check the upper range
+        if (celli >= nCells)
+        {
+            break;
+        }
+        else
+        {
+            agglom[celli] = celli;
+            ++nCompact;
+        }
     }
 
-    for (const polyPatch& pp : pbm)
+    // Second pass - finalize mappings
+    if (nCompact)
     {
-        if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
+        cellMap.resize(nCompact);
+        nCompact = 0;
+
+        for (label& celli : agglom)
         {
-            label bFacei = pp.start()-mesh.nInternalFaces();
-
-            for (const label facei : pp.range())
+            if (celli >= 0)
             {
-                const label own = agglom[faceOwner[facei]];
-                const label globalNei = globalNeighbour[bFacei];
+                cellMap[nCompact] = celli;
+                celli = nCompact;
+                ++nCompact;
 
-                if
-                (
-                   !globalAgglom.isLocal(myProci, globalNei)
-                 || globalAgglom.toLocal(myProci, globalNei) != own
-                )
+                if (nCompact == cellMap.size())
                 {
-                    ++nFacesPerCell[own];
+                    break;  // Early termination
                 }
-
-                ++bFacei;
             }
         }
     }
 
+    globalMeshData::calcCellCells
+    (
+        mesh,
+        agglom,
+        nCompact,   // == cellMap.size()
+        parallel,
+        cellCells
+    );
 
-    // Fill in offset and data
-    // ~~~~~~~~~~~~~~~~~~~~~~~
+    return cellMap;
+}
 
-    cellCells.resize_nocopy(nFacesPerCell);
-    cellCellWeights.resize_nocopy(nFacesPerCell);
 
-    nFacesPerCell = 0;  // Restart the count
+Foam::labelList Foam::globalMeshData::calcCellCells
+(
+    const polyMesh& mesh,
+    const labelUList& selectedCells,
+    CompactListList<label>& cellCells,
+    const bool parallel
+)
+{
+    const label nCells = mesh.nCells();
 
-    labelList& m = cellCells.values();
-    scalarList& w = cellCellWeights.values();
-    const labelList& offsets = cellCells.offsets();
+    labelList agglom(nCells, -1);
+    labelList cellMap;
 
-    // For internal faces is just offsetted owner and neighbour
-    for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
+    // First pass - creates a sorted order without duplicates
+    label nCompact = 0;
+    for (const label celli : selectedCells)
     {
-        const label own = agglom[faceOwner[facei]];
-        const label nei = agglom[faceNeigh[facei]];
-
-        const label ownIndex = offsets[own] + nFacesPerCell[own]++;
-        const label neiIndex = offsets[nei] + nFacesPerCell[nei]++;
-
-        m[ownIndex] = globalAgglom.toGlobal(myProci, nei);
-        m[neiIndex] = globalAgglom.toGlobal(myProci, own);
-
-        w[ownIndex] = mag(mesh.faceAreas()[facei]);
-        w[neiIndex] = w[ownIndex];
+        // Check cell is in range, and squash out duplicates
+        if (celli >= 0 && celli < nCells && agglom[celli] < 0)
+        {
+            agglom[celli] = celli;
+            ++nCompact;
+        }
     }
 
-    // For boundary faces is offsetted coupled neighbour
-    for (const polyPatch& pp : pbm)
+    // Second pass - finalize mappings
+    if (nCompact)
     {
-        if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
+        cellMap.resize(nCompact);
+        nCompact = 0;
+
+        for (label& celli : agglom)
         {
-            label bFacei = pp.start()-mesh.nInternalFaces();
-
-            for (const label facei : pp.range())
+            if (celli >= 0)
             {
-                const label own = agglom[faceOwner[facei]];
-                const label globalNei = globalNeighbour[bFacei];
+                cellMap[nCompact] = celli;
+                celli = nCompact;
+                ++nCompact;
 
-                if
-                (
-                   !globalAgglom.isLocal(myProci, globalNei)
-                 || globalAgglom.toLocal(myProci, globalNei) != own
-                )
+                if (nCompact == cellMap.size())
                 {
-                    const label ownIndex = offsets[own] + nFacesPerCell[own]++;
-
-                    m[ownIndex] = globalNei;
-                    w[ownIndex] = mag(mesh.faceAreas()[facei]);
+                    break;  // Early termination
                 }
-
-                ++bFacei;
             }
         }
     }
 
+    globalMeshData::calcCellCells
+    (
+        mesh,
+        agglom,
+        nCompact,   // == cellMap.size()
+        parallel,
+        cellCells
+    );
 
-    // Check for duplicates connections between cells
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Done as postprocessing step since we now have cellCells.
-
-    // NB: Because of agglomeration, self-connections will occur
-    //     and must be filtered out.
-
-    if (!cellCells.empty())
-    {
-        label newIndex = 0;
-        labelHashSet nbrCells;
-
-        labelList& m = cellCells.values();
-        scalarList& w = cellCellWeights.values();
-        labelList& offsets = cellCells.offsets();
-
-        label startIndex = offsets[0];
-
-        const label nCellCells = cellCells.size();
-
-        for (label celli = 0; celli < nCellCells; ++celli)
-        {
-            const label self = globalAgglom.toGlobal(myProci, celli);
-
-            nbrCells.clear();
-            nbrCells.insert(self);
-
-            const label endIndex = offsets[celli+1];
-
-            for (label i = startIndex; i < endIndex; ++i)
-            {
-                if (nbrCells.insert(m[i]))
-                {
-                    m[newIndex] = m[i];
-                    w[newIndex] = w[i];
-                    ++newIndex;
-                }
-            }
-            startIndex = endIndex;
-            offsets[celli+1] = newIndex;
-        }
-
-        m.resize(newIndex);
-        w.resize(newIndex);
-
-        // Weights has identical offsets as cellCells
-        cellCellWeights.offsets() = cellCells.offsets();
-    }
+    return cellMap;
 }
 
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCells.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCells.C
index a1924d6388..e28eb7e53f 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCells.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCells.C
@@ -60,7 +60,7 @@ void Foam::primitiveMesh::calcCellCells() const
     {
         // 1. Count number of internal faces per cell
 
-        labelList ncc(nCells(), Zero);
+        labelList ncc(nCells(), Foam::zero{});
 
         const labelList& own = faceOwner();
         const labelList& nei = faceNeighbour();
@@ -73,17 +73,16 @@ void Foam::primitiveMesh::calcCellCells() const
 
         // Create the storage
         ccPtr_ = new labelListList(ncc.size());
-        labelListList& cellCellAddr = *ccPtr_;
-
+        auto& cellCellAddr = *ccPtr_;
 
 
         // 2. Size and fill cellFaceAddr
 
         forAll(cellCellAddr, celli)
         {
-            cellCellAddr[celli].setSize(ncc[celli]);
+            cellCellAddr[celli].resize(ncc[celli]);
+            ncc[celli] = 0;  // reset size counter
         }
-        ncc = 0;
 
         forAll(nei, facei)
         {
diff --git a/src/parallel/decompose/decompose/decompositionInformation.C b/src/parallel/decompose/decompose/decompositionInformation.C
index ba3bc31a2c..8b4a31b7b4 100644
--- a/src/parallel/decompose/decompose/decompositionInformation.C
+++ b/src/parallel/decompose/decompose/decompositionInformation.C
@@ -40,17 +40,16 @@ void Foam::decompositionInformation::populate
 {
     nDomains_ = nDomain;
 
-    distrib_.clear();
-    distrib_.resize(nDomain);
+    distrib_.resize_nocopy(nDomain);
 
     for (labelList& subdist : distrib_)
     {
-        subdist.clear();
-        subdist.resize(nDomain, Zero);
+        subdist.resize_nocopy(nDomain);
+        subdist = Foam::zero{};
     }
 
     // Protect against zero-sized offset list
-    const label numCells = max(0, (xadj.size()-1));
+    const label numCells = Foam::max(0, (xadj.size()-1));
 
     for (label celli = 0; celli < numCells; ++celli)
     {
@@ -75,9 +74,9 @@ void Foam::decompositionInformation::populate
 
     // Build summary
 
-    labelList cellsCount(nDomains_, Zero);
-    labelList neighCount(nDomains_, Zero);
-    labelList facesCount(nDomains_, Zero);
+    labelList cellsCount(nDomains_, Foam::zero{});
+    labelList neighCount(nDomains_, Foam::zero{});
+    labelList facesCount(nDomains_, Foam::zero{});
 
     forAll(distrib_, ownProc)
     {
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
index 355899fe40..ba67ac8902 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -413,7 +413,7 @@ Foam::labelList Foam::decompositionMethod::decompose
         mesh,
         fineToCoarse,
         coarsePoints.size(),
-        true,                       // use global cell labels
+        true,        // Global mesh connectivity
         coarseCellCells
     );
 
@@ -589,7 +589,7 @@ Foam::labelList Foam::decompositionMethod::decompose
     // If we average the region centre instead, cyclics could cause
     // the average domain centre to be outside of domain.
 
-    scalarField regionWeights(localRegion.nLocalRegions(), Zero);
+    scalarField regionWeights(localRegion.nLocalRegions(), Foam::zero{});
 
     pointField regionCentres(localRegion.nLocalRegions(), point::max);
 
diff --git a/src/parallel/decompose/decompositionMethods/metisLikeDecomp/metisLikeDecomp.C b/src/parallel/decompose/decompositionMethods/metisLikeDecomp/metisLikeDecomp.C
index 210dcde5f3..2a07ac2a04 100644
--- a/src/parallel/decompose/decompositionMethods/metisLikeDecomp/metisLikeDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/metisLikeDecomp/metisLikeDecomp.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2023 OpenCFD Ltd.
+    Copyright (C) 2017-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,7 +59,7 @@ Foam::label Foam::metisLikeDecomp::decomposeGeneral
     }
 
     // Protect against zero-sized offset list
-    const label numCells = max(0, (xadj.size()-1));
+    const label numCells = Foam::max(0, (xadj.size()-1));
 
     const globalIndex globalAdjncy(adjncy.size());
     const globalIndex globalCells(numCells);
@@ -224,17 +224,10 @@ Foam::labelList Foam::metisLikeDecomp::decompose
             << exit(FatalError);
     }
 
+    // Global mesh connectivity
     CompactListList<label> cellCells;
-    globalMeshData::calcCellCells
-    (
-        mesh,
-        identity(mesh.nCells()),
-        mesh.nCells(),
-        true,
-        cellCells
-    );
+    globalMeshData::calcCellCells(mesh, cellCells, true);
 
-    // Decompose using default weights
     labelList decomp;
     decomposeGeneral
     (
@@ -274,11 +267,10 @@ Foam::labelList Foam::metisLikeDecomp::decompose
         mesh,
         agglom,
         agglomPoints.size(),
-        true,
+        true,        // Global mesh connectivity
         cellCells
     );
 
-    // Decompose using default weights
     labelList decomp;
     decomposeGeneral
     (
@@ -313,7 +305,6 @@ Foam::labelList Foam::metisLikeDecomp::decompose
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
 
-    // Decompose using default weights
     labelList decomp;
     decomposeGeneral
     (
@@ -349,7 +340,6 @@ Foam::labelList Foam::metisLikeDecomp::decompose
 
     auto cellCells(CompactListList<label>::pack(globalCellCells));
 
-    // Decompose using default weights
     labelList decomp;
     decomposeGeneral
     (
diff --git a/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C b/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C
index d3d0f3469b..b91df7c6d6 100644
--- a/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017-2023 OpenCFD Ltd.
+    Copyright (C) 2017-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -383,9 +383,6 @@ void Foam::multiLevelDecomp::decompose
     {
         // Recurse
 
-        // Determine points per domain
-        labelListList domainToPoints(invertOneToMany(nCurrDomains, dist));
-
         // Extract processor+local index from point-point addressing
         if (debug && Pstream::master())
         {
@@ -510,7 +507,7 @@ void Foam::multiLevelDecomp::decompose
                 // Count the number inbetween blocks of nNext size
 
                 label nPoints = 0;
-                labelList nOutsideConnections(nCurrDomains, Zero);
+                labelList nOutsideConnections(nCurrDomains, Foam::zero{});
                 forAll(pointPoints, pointi)
                 {
                     if ((dist[pointi] / nNext) == blockI)
@@ -611,19 +608,19 @@ Foam::labelList Foam::multiLevelDecomp::decompose
     const scalarField& cWeights
 ) const
 {
+    labelList finalDecomp(cc.size(), Foam::zero{});
+    labelList cellMap(Foam::identity(cc.size()));
+
     CompactListList<label> cellCells;
     globalMeshData::calcCellCells
     (
         mesh,
-        identity(cc.size()),
-        cc.size(),
-        true,
+        cellMap,
+        cellMap.size(),
+        true,       // Global mesh connectivity
         cellCells
     );
 
-    labelList finalDecomp(cc.size(), Foam::zero{});
-    labelList cellMap(identity(cc.size()));
-
     decompose
     (
         cellCells.unpack(),
@@ -648,7 +645,7 @@ Foam::labelList Foam::multiLevelDecomp::decompose
 ) const
 {
     labelList finalDecomp(cc.size(), Foam::zero{});
-    labelList cellMap(identity(cc.size()));
+    labelList cellMap(Foam::identity(cc.size()));
 
     decompose
     (
@@ -674,7 +671,7 @@ Foam::labelList Foam::multiLevelDecomp::decompose
 ) const
 {
     labelList finalDecomp(cc.size(), Foam::zero{});
-    labelList cellMap(identity(cc.size()));
+    labelList cellMap(Foam::identity(cc.size()));
 
     decompose
     (
diff --git a/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.H b/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.H
index 468af9c45d..5baeca557d 100644
--- a/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.H
@@ -72,7 +72,7 @@ class multiLevelDecomp
 
 
         //- Given connectivity across processors work out connectivity
-        //  for a (consistent) subset
+        //- for a (consistent) subset
         void subsetGlobalCellCells
         (
             const label nDomains,
diff --git a/src/parallel/decompose/kahipDecomp/kahipDecomp.C b/src/parallel/decompose/kahipDecomp/kahipDecomp.C
index b889efa020..b2341cd69c 100644
--- a/src/parallel/decompose/kahipDecomp/kahipDecomp.C
+++ b/src/parallel/decompose/kahipDecomp/kahipDecomp.C
@@ -101,7 +101,7 @@ Foam::label Foam::kahipDecomp::decomposeSerial
     }
     #endif
 
-    int numCells = max(0, (xadj.size()-1));
+    int numCells = Foam::max(0, (xadj.size()-1));
 
     // Addressing
     ConstPrecisionAdaptor<int, label, List> adjncy_param(adjncy);
diff --git a/src/parallel/decompose/metisDecomp/metisDecomp.C b/src/parallel/decompose/metisDecomp/metisDecomp.C
index d62307cb48..4557c6effc 100644
--- a/src/parallel/decompose/metisDecomp/metisDecomp.C
+++ b/src/parallel/decompose/metisDecomp/metisDecomp.C
@@ -85,7 +85,7 @@ Foam::label Foam::metisDecomp::decomposeSerial
 
     const dictionary* coeffsDictPtr = decompDict_.findDict("metisCoeffs");
 
-    idx_t numCells = max(0, (xadj.size()-1));
+    idx_t numCells = Foam::max(0, (xadj.size()-1));
 
     // Decomposition options
     List<idx_t> options(METIS_NOPTIONS);
diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
index 1fb105cf4f..f5e1d9b0bf 100644
--- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
+++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2023 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -117,7 +117,7 @@ Foam::label Foam::ptscotchDecomp::decompose
     labelList& decomp
 ) const
 {
-    const SCOTCH_Num numCells = max(0, (xadj.size()-1));
+    const SCOTCH_Num numCells = Foam::max(0, (xadj.size()-1));
 
     // Addressing
     ConstPrecisionAdaptor<SCOTCH_Num, label, List> adjncy_param(adjncy);
@@ -559,17 +559,10 @@ Foam::labelList Foam::ptscotchDecomp::decompose
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
 
+    // Global mesh connectivity
     CompactListList<label> cellCells;
-    globalMeshData::calcCellCells
-    (
-        mesh,
-        identity(mesh.nCells()),
-        mesh.nCells(),
-        true,
-        cellCells
-    );
+    globalMeshData::calcCellCells(mesh, cellCells, true);
 
-    // Decompose using default weights
     labelList decomp;
     decompose
     (
@@ -606,13 +599,14 @@ Foam::labelList Foam::ptscotchDecomp::decompose
     // Make Metis CSR (Compressed Storage Format) storage
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
+
     CompactListList<label> cellCells;
     globalMeshData::calcCellCells
     (
         mesh,
         agglom,
         agglomPoints.size(),
-        true,
+        true,       // Global mesh connectivity
         cellCells
     );
 
@@ -654,7 +648,6 @@ Foam::labelList Foam::ptscotchDecomp::decompose
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
 
-    // Decompose using default weights
     labelList decomp;
     decompose
     (
@@ -692,7 +685,6 @@ Foam::labelList Foam::ptscotchDecomp::decompose
 
     auto cellCells(CompactListList<label>::pack(globalCellCells));
 
-    // Decompose using default weights
     labelList decomp;
     decompose
     (
diff --git a/src/parallel/decompose/scotchDecomp/scotchDecomp.C b/src/parallel/decompose/scotchDecomp/scotchDecomp.C
index c3637bd4dd..b2e4dc65d4 100644
--- a/src/parallel/decompose/scotchDecomp/scotchDecomp.C
+++ b/src/parallel/decompose/scotchDecomp/scotchDecomp.C
@@ -113,7 +113,7 @@ Foam::label Foam::scotchDecomp::decomposeSerial
     labelList& decomp
 ) const
 {
-    const SCOTCH_Num numCells = max(0, (xadj.size()-1));
+    const SCOTCH_Num numCells = Foam::max(0, (xadj.size()-1));
 
     // Addressing
     ConstPrecisionAdaptor<SCOTCH_Num, label, List> adjncy_param(adjncy);
diff --git a/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C b/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
index b2662f90f6..c5a867decd 100644
--- a/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
+++ b/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2020-2023 OpenCFD Ltd.
+    Copyright (C) 2020-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,7 +29,6 @@ License
 #include "CuthillMcKeeRenumber.H"
 #include "addToRunTimeSelectionTable.H"
 #include "bandCompression.H"
-#include "globalMeshData.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -67,17 +66,7 @@ Foam::labelList Foam::CuthillMcKeeRenumber::renumber
     const pointField& points
 ) const
 {
-    CompactListList<label> cellCells;
-    globalMeshData::calcCellCells
-    (
-        mesh,
-        identity(mesh.nCells()),
-        mesh.nCells(),
-        false,                      // local only
-        cellCells
-    );
-
-    labelList orderedToOld = meshTools::bandCompression(cellCells);
+    labelList orderedToOld = meshTools::bandCompression(mesh);
 
     if (reverse_)
     {
diff --git a/src/renumber/renumberMethods/renumberMethod/renumberMethod.C b/src/renumber/renumberMethods/renumberMethod/renumberMethod.C
index 2ea1cf0c7e..7be4654d92 100644
--- a/src/renumber/renumberMethods/renumberMethod/renumberMethod.C
+++ b/src/renumber/renumberMethods/renumberMethod/renumberMethod.C
@@ -73,15 +73,9 @@ Foam::labelList Foam::renumberMethod::renumber
     const pointField& points
 ) const
 {
+    // Local mesh connectivity
     CompactListList<label> cellCells;
-    globalMeshData::calcCellCells
-    (
-        mesh,
-        identity(mesh.nCells()),
-        mesh.nCells(),
-        false,                      // local only
-        cellCells
-    );
+    globalMeshData::calcCellCells(mesh, cellCells);
 
     return renumber(cellCells, points);
 }
@@ -122,7 +116,7 @@ Foam::labelList Foam::renumberMethod::renumber
         mesh,
         fineToCoarse,
         coarsePoints.size(),
-        false,                      // local only
+        false,  // local only (parallel = false)
         coarseCellCells
     );
 
diff --git a/src/renumber/renumberMethods/springRenumber/springRenumber.C b/src/renumber/renumberMethods/springRenumber/springRenumber.C
index f7b79dfc93..c331b23ad8 100644
--- a/src/renumber/renumberMethods/springRenumber/springRenumber.C
+++ b/src/renumber/renumberMethods/springRenumber/springRenumber.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -79,7 +79,7 @@ Foam::labelList Foam::springRenumber::renumberImpl
     // Sum force per cell. Also reused for the displacement
     scalarField sumForce(nOldCells);
 
-    labelList oldToNew(identity(nOldCells));
+    labelList oldToNew(Foam::identity(nOldCells));
 
     scalar maxCo = (maxCo_ * nOldCells);
 
@@ -163,15 +163,9 @@ Foam::labelList Foam::springRenumber::renumber
     const pointField&
 ) const
 {
+    // Local mesh connectivity
     CompactListList<label> cellCells;
-    globalMeshData::calcCellCells
-    (
-        mesh,
-        identity(mesh.nCells()),
-        mesh.nCells(),
-        false,                      // local only
-        cellCells
-    );
+    globalMeshData::calcCellCells(mesh, cellCells);
 
     return renumberImpl(cellCells);
 }