diff --git a/applications/test/checkDecomposePar/Test-checkDecomposePar.C b/applications/test/checkDecomposePar/Test-checkDecomposePar.C
index d616b013c7..8019fd0d45 100644
--- a/applications/test/checkDecomposePar/Test-checkDecomposePar.C
+++ b/applications/test/checkDecomposePar/Test-checkDecomposePar.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,6 +40,7 @@ Description
 #include "cpuTime.H"
 #include "IFstream.H"
 #include "regionProperties.H"
+#include "globalMeshData.H"
 #include "decompositionInformation.H"
 #include "decompositionModel.H"
 
@@ -123,7 +124,7 @@ int main(int argc, char *argv[])
         const label nDomains = max(cellToProc) + 1;
 
         CompactListList<label> cellCells;
-        decompositionMethod::calcCellCells
+        globalMeshData::calcCellCells
         (
             mesh,
             identity(mesh.nCells()),
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.C b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.C
index c80f791515..36d04152d7 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-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,6 +27,7 @@ License
 
 #include "domainDecompositionDryRun.H"
 #include "volFields.H"
+#include "globalMeshData.H"
 #include "decompositionModel.H"
 #include "decompositionInformation.H"
 
@@ -112,12 +113,12 @@ void Foam::domainDecompositionDryRun::execute
     decompositionMethod& method = model.decomposer();
 
     CompactListList<label> cellCells;
-    decompositionMethod::calcCellCells
+    globalMeshData::calcCellCells
     (
         mesh_,
         identity(mesh_.nCells()),
         mesh_.nCells(),
-        false,
+        false,  // false = local only
         cellCells
     );
 
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 54ffbe8a71..5aa90ccedc 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -614,6 +614,7 @@ $(polyBoundaryMesh)/polyBoundaryMeshEntries.C
 
 globalMeshData = $(polyMesh)/globalMeshData
 $(globalMeshData)/globalMeshData.C
+$(globalMeshData)/globalMeshDataTopology.C
 $(globalMeshData)/globalPoints.C
 
 $(polyMesh)/syncTools/syncTools.C
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
index 40af50afbe..a1228f6545 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
@@ -96,6 +96,7 @@ namespace Foam
 // Forward Declarations
 class polyMesh;
 class mapDistribute;
+template<class T> class CompactListList;
 template<class T> class EdgeMap;
 class globalIndex;
 class globalIndexAndTransform;
@@ -305,6 +306,9 @@ class globalMeshData
 
             void calcGlobalCoPointSlaves() const;
 
+public:
+
+    // Generated Methods
 
         //- No copy construct
         globalMeshData(const globalMeshData&) = delete;
@@ -313,8 +317,6 @@ class globalMeshData
         void operator=(const globalMeshData&) = delete;
 
 
-public:
-
     //- Runtime type information
     ClassName("globalMeshData");
 
@@ -606,6 +608,39 @@ public:
             //  full parallel analysis to determine shared points and
             //  boundaries.
             void updateMesh();
+
+
+    // Mesh Topology Calculation
+
+        //- 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.
+        //  - parallel = true
+        //    Resulting connections are in global cell indices.
+        //    Coupled across cyclics and processor patches.
+        static void calcCellCells
+        (
+            const polyMesh& mesh,
+            const labelList& agglom,
+            const label nLocalCoarse,
+            const bool parallel,  //!< Use global cell ids in parallel
+            CompactListList<label>& cellCells
+        );
+
+        //- Determine (local or global) cellCells and face weights
+        //- from mesh agglomeration.
+        //  Uses mag of faceArea as weights
+        static void calcCellCells
+        (
+            const polyMesh& mesh,
+            const labelList& agglom,
+            const label nLocalCoarse,
+            const bool parallel,  //!< Use global cell ids in parallel
+            CompactListList<label>& cellCells,
+            CompactListList<scalar>& cellCellWeights
+        );
 };
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshDataTopology.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshDataTopology.C
new file mode 100644
index 0000000000..7743b75714
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshDataTopology.C
@@ -0,0 +1,401 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2015-2023 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "globalMeshData.H"
+#include "globalIndex.H"
+#include "CompactListList.H"
+#include "processorPolyPatch.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+void Foam::globalMeshData::calcCellCells
+(
+    const polyMesh& mesh,
+    const labelList& agglom,
+    const label nLocalCoarse,
+    const bool parallel,
+    CompactListList<label>& cellCells
+)
+{
+    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);
+    syncTools::swapBoundaryFaceList(mesh, globalNeighbour);
+
+
+    // Count number of faces (internal + coupled)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Number of faces per coarse cell
+    labelList nFacesPerCell(nLocalCoarse, Zero);
+
+    for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
+    {
+        const label own = agglom[faceOwner[facei]];
+        const label nei = agglom[faceNeigh[facei]];
+
+        ++nFacesPerCell[own];
+        ++nFacesPerCell[nei];
+    }
+
+    for (const polyPatch& pp : pbm)
+    {
+        if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
+        {
+            label bFacei = pp.start()-mesh.nInternalFaces();
+
+            for (const label facei : pp.range())
+            {
+                const label own = agglom[faceOwner[facei]];
+                const label globalNei = globalNeighbour[bFacei];
+
+                if
+                (
+                   !globalAgglom.isLocal(myProci, globalNei)
+                 || globalAgglom.toLocal(myProci, globalNei) != own
+                )
+                {
+                    ++nFacesPerCell[own];
+                }
+
+                ++bFacei;
+            }
+        }
+    }
+
+
+    // Fill in offset and data
+    // ~~~~~~~~~~~~~~~~~~~~~~~
+
+    cellCells.resize_nocopy(nFacesPerCell);
+
+    nFacesPerCell = 0;  // Restart the count
+
+    labelList& m = cellCells.values();
+    const labelList& offsets = cellCells.offsets();
+
+    // 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]];
+
+        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);
+    }
+
+    // For boundary faces is offsetted coupled neighbour
+    for (const polyPatch& pp : pbm)
+    {
+        if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
+        {
+            label bFacei = pp.start()-mesh.nInternalFaces();
+
+            for (const label facei : pp.range())
+            {
+                const label own = agglom[faceOwner[facei]];
+                const label globalNei = globalNeighbour[bFacei];
+
+                if
+                (
+                   !globalAgglom.isLocal(myProci, globalNei)
+                 || globalAgglom.toLocal(myProci, globalNei) != own
+                )
+                {
+                    const label ownIndex = offsets[own] + nFacesPerCell[own]++;
+
+                    m[ownIndex] = globalNei;
+                }
+
+                ++bFacei;
+            }
+        }
+    }
+
+
+    // 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();
+        labelList& offsets = cellCells.offsets();
+
+        label startIndex = offsets[0];
+
+        const label nCellCells = cellCells.size();
+
+        for (label celli = 0; celli < nCellCells; ++celli)
+        {
+            // Could be done as combination of std::sort, std::unique_copy
+            // except need to block out 'self' and std::copy_if/remove_if
+            // are undefined for overlapping regions
+
+            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];
+                    ++newIndex;
+                }
+            }
+            startIndex = endIndex;
+            offsets[celli+1] = newIndex;
+        }
+
+        m.resize(newIndex);
+    }
+
+    //forAll(cellCells, celli)
+    //{
+    //    Pout<< "Original: Coarse cell " << celli << endl;
+    //    forAll(mesh.cellCells()[celli], i)
+    //    {
+    //        Pout<< "    nbr:" << mesh.cellCells()[celli][i] << endl;
+    //    }
+    //    Pout<< "Compacted: Coarse cell " << celli << endl;
+    //    const labelUList cCells = cellCells[celli];
+    //    forAll(cCells, i)
+    //    {
+    //        Pout<< "    nbr:" << cCells[i] << endl;
+    //    }
+    //}
+}
+
+
+void Foam::globalMeshData::calcCellCells
+(
+    const polyMesh& mesh,
+    const labelList& 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);
+
+
+    // The agglomerated owner per boundary faces (global numbering)
+    // from the other side (of coupled patches)
+
+    labelList globalNeighbour(agglom, pbm.faceOwner());
+    globalAgglom.inplaceToGlobal(myProci, globalNeighbour);
+    syncTools::swapBoundaryFaceList(mesh, globalNeighbour);
+
+
+    // Count number of faces (internal + coupled)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Number of faces per coarse cell
+    labelList nFacesPerCell(nLocalCoarse, Zero);
+
+    for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
+    {
+        const label own = agglom[faceOwner[facei]];
+        const label nei = agglom[faceNeigh[facei]];
+
+        ++nFacesPerCell[own];
+        ++nFacesPerCell[nei];
+    }
+
+    for (const polyPatch& pp : pbm)
+    {
+        if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
+        {
+            label bFacei = pp.start()-mesh.nInternalFaces();
+
+            for (const label facei : pp.range())
+            {
+                const label own = agglom[faceOwner[facei]];
+                const label globalNei = globalNeighbour[bFacei];
+
+                if
+                (
+                   !globalAgglom.isLocal(myProci, globalNei)
+                 || globalAgglom.toLocal(myProci, globalNei) != own
+                )
+                {
+                    ++nFacesPerCell[own];
+                }
+
+                ++bFacei;
+            }
+        }
+    }
+
+
+    // Fill in offset and data
+    // ~~~~~~~~~~~~~~~~~~~~~~~
+
+    cellCells.resize_nocopy(nFacesPerCell);
+    cellCellWeights.resize_nocopy(nFacesPerCell);
+
+    nFacesPerCell = 0;  // Restart the count
+
+    labelList& m = cellCells.values();
+    scalarList& w = cellCellWeights.values();
+    const labelList& offsets = cellCells.offsets();
+
+    // 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]];
+
+        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];
+    }
+
+    // For boundary faces is offsetted coupled neighbour
+    for (const polyPatch& pp : pbm)
+    {
+        if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
+        {
+            label bFacei = pp.start()-mesh.nInternalFaces();
+
+            for (const label facei : pp.range())
+            {
+                const label own = agglom[faceOwner[facei]];
+                const label globalNei = globalNeighbour[bFacei];
+
+                if
+                (
+                   !globalAgglom.isLocal(myProci, globalNei)
+                 || globalAgglom.toLocal(myProci, globalNei) != own
+                )
+                {
+                    const label ownIndex = offsets[own] + nFacesPerCell[own]++;
+
+                    m[ownIndex] = globalNei;
+                    w[ownIndex] = mag(mesh.faceAreas()[facei]);
+                }
+
+                ++bFacei;
+            }
+        }
+    }
+
+
+    // 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();
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/dummyThirdParty/MGridGen/Make/files b/src/dummyThirdParty/MGridGen/Make/files
index 684d67331a..05f7baef2c 100644
--- a/src/dummyThirdParty/MGridGen/Make/files
+++ b/src/dummyThirdParty/MGridGen/Make/files
@@ -1,3 +1,3 @@
-dummyMGridGen.C
+dummyMGridGen.cxx
 
 LIB = $(FOAM_LIBBIN)/dummy/libMGridGen
diff --git a/src/dummyThirdParty/MGridGen/dummyMGridGen.C b/src/dummyThirdParty/MGridGen/dummyMGridGen.cxx
similarity index 100%
rename from src/dummyThirdParty/MGridGen/dummyMGridGen.C
rename to src/dummyThirdParty/MGridGen/dummyMGridGen.cxx
diff --git a/src/dummyThirdParty/MGridGen/mgridgen.h b/src/dummyThirdParty/MGridGen/mgridgen.h
index 0c50137f85..990965add8 100644
--- a/src/dummyThirdParty/MGridGen/mgridgen.h
+++ b/src/dummyThirdParty/MGridGen/mgridgen.h
@@ -31,14 +31,14 @@ Description
     Only implements the absolute minimum we are using.
 
 SourceFiles
-    dummyMGridGen.C
+    dummyMGridGen.cxx
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef mgridgen_H
 #define mgridgen_H
 
-#include "scalar.H"
+#include "scalarFwd.H"
 
 #ifndef idxtype
 #define idxtype int
diff --git a/src/dummyThirdParty/kahipDecomp/Make/files b/src/dummyThirdParty/kahipDecomp/Make/files
index 1a02275676..fdf8f0abee 100644
--- a/src/dummyThirdParty/kahipDecomp/Make/files
+++ b/src/dummyThirdParty/kahipDecomp/Make/files
@@ -1,3 +1,3 @@
-dummyKahipDecomp.C
+dummyKahipDecomp.cxx
 
 LIB = $(FOAM_LIBBIN)/dummy/libkahipDecomp
diff --git a/src/dummyThirdParty/kahipDecomp/dummyKahipDecomp.C b/src/dummyThirdParty/kahipDecomp/dummyKahipDecomp.cxx
similarity index 95%
rename from src/dummyThirdParty/kahipDecomp/dummyKahipDecomp.C
rename to src/dummyThirdParty/kahipDecomp/dummyKahipDecomp.cxx
index fd7ebcf39c..d1d4aa88ad 100644
--- a/src/dummyThirdParty/kahipDecomp/dummyKahipDecomp.C
+++ b/src/dummyThirdParty/kahipDecomp/dummyKahipDecomp.cxx
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -70,6 +70,12 @@ Foam::label Foam::kahipDecomp::decomposeSerial
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::kahipDecomp::kahipDecomp(const label numDomains)
+:
+    metisLikeDecomp(numDomains)
+{}
+
+
 Foam::kahipDecomp::kahipDecomp
 (
     const dictionary& decompDict,
diff --git a/src/dummyThirdParty/metisDecomp/Make/files b/src/dummyThirdParty/metisDecomp/Make/files
index cdc9482e7d..a22272821a 100644
--- a/src/dummyThirdParty/metisDecomp/Make/files
+++ b/src/dummyThirdParty/metisDecomp/Make/files
@@ -1,3 +1,3 @@
-dummyMetisDecomp.C
+dummyMetisDecomp.cxx
 
 LIB = $(FOAM_LIBBIN)/dummy/libmetisDecomp
diff --git a/src/dummyThirdParty/metisDecomp/dummyMetisDecomp.C b/src/dummyThirdParty/metisDecomp/dummyMetisDecomp.cxx
similarity index 100%
rename from src/dummyThirdParty/metisDecomp/dummyMetisDecomp.C
rename to src/dummyThirdParty/metisDecomp/dummyMetisDecomp.cxx
diff --git a/src/dummyThirdParty/ptscotchDecomp/Make/files b/src/dummyThirdParty/ptscotchDecomp/Make/files
index 18402ebfca..a6920f5b7e 100644
--- a/src/dummyThirdParty/ptscotchDecomp/Make/files
+++ b/src/dummyThirdParty/ptscotchDecomp/Make/files
@@ -1,3 +1,3 @@
-dummyPtscotchDecomp.C
+dummyPtscotchDecomp.cxx
 
 LIB = $(FOAM_LIBBIN)/dummy/libptscotchDecomp
diff --git a/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C b/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.cxx
similarity index 88%
rename from src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C
rename to src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.cxx
index 63b23ce1cb..346d7d4410 100644
--- a/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C
+++ b/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.cxx
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -73,6 +73,13 @@ Foam::label Foam::ptscotchDecomp::decompose
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::ptscotchDecomp::ptscotchDecomp(const label numDomains)
+:
+    decompositionMethod(numDomains),
+    coeffsDict_()
+{}
+
+
 Foam::ptscotchDecomp::ptscotchDecomp
 (
     const dictionary& decompDict,
@@ -105,7 +112,21 @@ Foam::labelList Foam::ptscotchDecomp::decompose
     const polyMesh& mesh,
     const labelList& agglom,
     const pointField& agglomPoints,
-    const scalarField& pointWeights
+    const scalarField& agglomWeights
+) const
+{
+    FatalErrorInFunction
+        << notImplementedMessage << exit(FatalError);
+
+    return labelList();
+}
+
+
+Foam::labelList Foam::ptscotchDecomp::decompose
+(
+    const CompactListList<label>& globalCellCells,
+    const pointField&,
+    const scalarField&
 ) const
 {
     FatalErrorInFunction
diff --git a/src/dummyThirdParty/scotchDecomp/Make/files b/src/dummyThirdParty/scotchDecomp/Make/files
index f8f7b50143..adb41cf103 100644
--- a/src/dummyThirdParty/scotchDecomp/Make/files
+++ b/src/dummyThirdParty/scotchDecomp/Make/files
@@ -1,3 +1,3 @@
-dummyScotchDecomp.C
+dummyScotchDecomp.cxx
 
 LIB = $(FOAM_LIBBIN)/dummy/libscotchDecomp
diff --git a/src/dummyThirdParty/scotchDecomp/dummyScotchDecomp.C b/src/dummyThirdParty/scotchDecomp/dummyScotchDecomp.cxx
similarity index 88%
rename from src/dummyThirdParty/scotchDecomp/dummyScotchDecomp.C
rename to src/dummyThirdParty/scotchDecomp/dummyScotchDecomp.cxx
index 9259f8a85d..fab2bf49b0 100644
--- a/src/dummyThirdParty/scotchDecomp/dummyScotchDecomp.C
+++ b/src/dummyThirdParty/scotchDecomp/dummyScotchDecomp.cxx
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -74,6 +74,12 @@ Foam::label Foam::scotchDecomp::decomposeSerial
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::scotchDecomp::scotchDecomp(const label numDomains)
+:
+    metisLikeDecomp(numDomains)
+{}
+
+
 Foam::scotchDecomp::scotchDecomp
 (
     const dictionary& decompDict,
@@ -105,7 +111,21 @@ Foam::labelList Foam::scotchDecomp::decompose
     const polyMesh& mesh,
     const labelList& agglom,
     const pointField& agglomPoints,
-    const scalarField& pointWeights
+    const scalarField& agglomWeights
+) const
+{
+    FatalErrorInFunction
+        << notImplementedMessage << exit(FatalError);
+
+    return labelList();
+}
+
+
+Foam::labelList Foam::scotchDecomp::decompose
+(
+    const CompactListList<label>& globalCellCells,
+    const pointField&,
+    const scalarField&
 ) const
 {
     FatalErrorInFunction
@@ -118,8 +138,8 @@ Foam::labelList Foam::scotchDecomp::decompose
 Foam::labelList Foam::scotchDecomp::decompose
 (
     const labelListList& globalCellCells,
-    const pointField& cellCentres,
-    const scalarField& cWeights
+    const pointField&,
+    const scalarField&
 ) const
 {
     FatalErrorInFunction
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
index 5851e8e1f3..355899fe40 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -28,6 +28,7 @@ License
 
 #include "decompositionMethod.H"
 #include "globalIndex.H"
+#include "globalMeshData.H"
 #include "syncTools.H"
 #include "faceSet.H"
 #include "regionSplit.H"
@@ -331,6 +332,14 @@ const Foam::dictionary& Foam::decompositionMethod::findCoeffsDict
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::decompositionMethod::decompositionMethod(const label numDomains)
+:
+    decompDict_(dictionary::null),
+    decompRegionDict_(dictionary::null),
+    nDomains_(numDomains)
+{}
+
+
 Foam::decompositionMethod::decompositionMethod
 (
     const dictionary& decompDict,
@@ -390,18 +399,6 @@ Foam::autoPtr<Foam::decompositionMethod> Foam::decompositionMethod::New
 }
 
 
-Foam::labelList Foam::decompositionMethod::decompose
-(
-    const polyMesh& mesh,
-    const pointField& points
-) const
-{
-    scalarField weights(points.size(), scalar(1));
-
-    return decompose(mesh, points, weights);
-}
-
-
 Foam::labelList Foam::decompositionMethod::decompose
 (
     const polyMesh& mesh,
@@ -411,7 +408,7 @@ Foam::labelList Foam::decompositionMethod::decompose
 ) const
 {
     CompactListList<label> coarseCellCells;
-    calcCellCells
+    globalMeshData::calcCellCells
     (
         mesh,
         fineToCoarse,
@@ -421,56 +418,18 @@ Foam::labelList Foam::decompositionMethod::decompose
     );
 
     // Decompose based on agglomerated points
-    labelList coarseDistribution
+    labelList decomp
     (
         decompose
         (
-            coarseCellCells.unpack(),
+            coarseCellCells,
             coarsePoints,
             coarseWeights
         )
     );
 
-    // Rework back into decomposition for original mesh_
-    labelList fineDistribution(fineToCoarse.size());
-
-    forAll(fineDistribution, i)
-    {
-        fineDistribution[i] = coarseDistribution[fineToCoarse[i]];
-    }
-
-    return fineDistribution;
-}
-
-
-Foam::labelList Foam::decompositionMethod::decompose
-(
-    const polyMesh& mesh,
-    const labelList& fineToCoarse,
-    const pointField& coarsePoints
-) const
-{
-    scalarField weights(coarsePoints.size(), scalar(1));
-
-    return decompose
-    (
-        mesh,
-        fineToCoarse,
-        coarsePoints,
-        weights
-    );
-}
-
-
-Foam::labelList Foam::decompositionMethod::decompose
-(
-    const labelListList& globalCellCells,
-    const pointField& cc
-) const
-{
-    scalarField weights(cc.size(), scalar(1));
-
-    return decompose(globalCellCells, cc, weights);
+    // From coarse back to fine for original mesh
+    return labelList(decomp, fineToCoarse);
 }
 
 
@@ -483,186 +442,14 @@ void Foam::decompositionMethod::calcCellCells
     CompactListList<label>& cellCells
 )
 {
-    const labelList& faceOwner = mesh.faceOwner();
-    const labelList& faceNeighbour = mesh.faceNeighbour();
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-
-    // Create global cell numbers
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    const globalIndex globalAgglom(nLocalCoarse, UPstream::worldComm, parallel);
-
-
-    // Get agglomerate owner on other side of coupled faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelList globalNeighbour(mesh.nBoundaryFaces());
-
-    for (const polyPatch& pp : patches)
-    {
-        if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
-        {
-            label facei = pp.start();
-            label bFacei = pp.start() - mesh.nInternalFaces();
-
-            forAll(pp, i)
-            {
-                globalNeighbour[bFacei] = globalAgglom.toGlobal
-                (
-                    agglom[faceOwner[facei]]
-                );
-
-                ++facei;
-                ++bFacei;
-            }
-        }
-    }
-
-    // Get the cell on the other side of coupled patches
-    syncTools::swapBoundaryFaceList(mesh, globalNeighbour);
-
-
-    // Count number of faces (internal + coupled)
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Number of faces per coarse cell
-    labelList nFacesPerCell(nLocalCoarse, Zero);
-
-    for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
-    {
-        const label own = agglom[faceOwner[facei]];
-        const label nei = agglom[faceNeighbour[facei]];
-
-        nFacesPerCell[own]++;
-        nFacesPerCell[nei]++;
-    }
-
-    for (const polyPatch& pp : patches)
-    {
-        if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
-        {
-            label facei = pp.start();
-            label bFacei = pp.start()-mesh.nInternalFaces();
-
-            forAll(pp, i)
-            {
-                const label own = agglom[faceOwner[facei]];
-                const label globalNei = globalNeighbour[bFacei];
-
-                if
-                (
-                   !globalAgglom.isLocal(globalNei)
-                 || globalAgglom.toLocal(globalNei) != own
-                )
-                {
-                    nFacesPerCell[own]++;
-                }
-
-                ++facei;
-                ++bFacei;
-            }
-        }
-    }
-
-
-    // Fill in offset and data
-    // ~~~~~~~~~~~~~~~~~~~~~~~
-
-    cellCells.setSize(nFacesPerCell);
-
-    nFacesPerCell = 0;
-
-    labelList& m = cellCells.m();
-    const labelList& offsets = cellCells.offsets();
-
-    // 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[faceNeighbour[facei]];
-
-        m[offsets[own] + nFacesPerCell[own]++] = globalAgglom.toGlobal(nei);
-        m[offsets[nei] + nFacesPerCell[nei]++] = globalAgglom.toGlobal(own);
-    }
-
-    // For boundary faces is offsetted coupled neighbour
-    for (const polyPatch& pp : patches)
-    {
-        if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
-        {
-            label facei = pp.start();
-            label bFacei = pp.start()-mesh.nInternalFaces();
-
-            forAll(pp, i)
-            {
-                const label own = agglom[faceOwner[facei]];
-                const label globalNei = globalNeighbour[bFacei];
-
-                if
-                (
-                   !globalAgglom.isLocal(globalNei)
-                 || globalAgglom.toLocal(globalNei) != own
-                )
-                {
-                    m[offsets[own] + nFacesPerCell[own]++] = globalNei;
-                }
-
-                ++facei;
-                ++bFacei;
-            }
-        }
-    }
-
-
-    // Check for duplicates connections between cells
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Done as postprocessing step since we now have cellCells.
-
-    if (cellCells.size() == 0)
-    {
-        return;
-    }
-
-    label newIndex = 0;
-    labelHashSet nbrCells;
-
-    label startIndex = cellCells.offsets()[0];
-
-    forAll(cellCells, celli)
-    {
-        nbrCells.clear();
-        nbrCells.insert(globalAgglom.toGlobal(celli));
-
-        const label endIndex = cellCells.offsets()[celli+1];
-
-        for (label i = startIndex; i < endIndex; ++i)
-        {
-            if (nbrCells.insert(cellCells.m()[i]))
-            {
-                cellCells.m()[newIndex++] = cellCells.m()[i];
-            }
-        }
-        startIndex = endIndex;
-        cellCells.offsets()[celli+1] = newIndex;
-    }
-
-    cellCells.m().setSize(newIndex);
-
-    //forAll(cellCells, celli)
-    //{
-    //    Pout<< "Original: Coarse cell " << celli << endl;
-    //    forAll(mesh.cellCells()[celli], i)
-    //    {
-    //        Pout<< "    nbr:" << mesh.cellCells()[celli][i] << endl;
-    //    }
-    //    Pout<< "Compacted: Coarse cell " << celli << endl;
-    //    const labelUList cCells = cellCells[celli];
-    //    forAll(cCells, i)
-    //    {
-    //        Pout<< "    nbr:" << cCells[i] << endl;
-    //    }
-    //}
+    globalMeshData::calcCellCells
+    (
+        mesh,
+        agglom,
+        nLocalCoarse,
+        parallel,
+        cellCells
+    );
 }
 
 
@@ -676,184 +463,15 @@ void Foam::decompositionMethod::calcCellCells
     CompactListList<scalar>& cellCellWeights
 )
 {
-    const labelList& faceOwner = mesh.faceOwner();
-    const labelList& faceNeighbour = mesh.faceNeighbour();
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-
-    // Create global cell numbers
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    const globalIndex globalAgglom(nLocalCoarse, UPstream::worldComm, parallel);
-
-
-    // Get agglomerate owner on other side of coupled faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelList globalNeighbour(mesh.nBoundaryFaces());
-
-    for (const polyPatch& pp : patches)
-    {
-        if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
-        {
-            label facei = pp.start();
-            label bFacei = pp.start() - mesh.nInternalFaces();
-
-            forAll(pp, i)
-            {
-                globalNeighbour[bFacei] = globalAgglom.toGlobal
-                (
-                    agglom[faceOwner[facei]]
-                );
-
-                ++facei;
-                ++bFacei;
-            }
-        }
-    }
-
-    // Get the cell on the other side of coupled patches
-    syncTools::swapBoundaryFaceList(mesh, globalNeighbour);
-
-
-    // Count number of faces (internal + coupled)
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Number of faces per coarse cell
-    labelList nFacesPerCell(nLocalCoarse, Zero);
-
-    for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
-    {
-        const label own = agglom[faceOwner[facei]];
-        const label nei = agglom[faceNeighbour[facei]];
-
-        nFacesPerCell[own]++;
-        nFacesPerCell[nei]++;
-    }
-
-    for (const polyPatch& pp : patches)
-    {
-        if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
-        {
-            label facei = pp.start();
-            label bFacei = pp.start() - mesh.nInternalFaces();
-
-            forAll(pp, i)
-            {
-                const label own = agglom[faceOwner[facei]];
-                const label globalNei = globalNeighbour[bFacei];
-
-                if
-                (
-                   !globalAgglom.isLocal(globalNei)
-                 || globalAgglom.toLocal(globalNei) != own
-                )
-                {
-                    nFacesPerCell[own]++;
-                }
-
-                ++facei;
-                ++bFacei;
-            }
-        }
-    }
-
-
-    // Fill in offset and data
-    // ~~~~~~~~~~~~~~~~~~~~~~~
-
-    cellCells.setSize(nFacesPerCell);
-    cellCellWeights.setSize(nFacesPerCell);
-
-    nFacesPerCell = 0;
-
-    labelList& m = cellCells.m();
-    scalarList& w = cellCellWeights.m();
-    const labelList& offsets = cellCells.offsets();
-
-    // 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[faceNeighbour[facei]];
-
-        const label ownIndex = offsets[own] + nFacesPerCell[own]++;
-        const label neiIndex = offsets[nei] + nFacesPerCell[nei]++;
-
-        m[ownIndex] = globalAgglom.toGlobal(nei);
-        w[ownIndex] = mag(mesh.faceAreas()[facei]);
-        m[neiIndex] = globalAgglom.toGlobal(own);
-        w[ownIndex] = mag(mesh.faceAreas()[facei]);
-    }
-
-    // For boundary faces is offsetted coupled neighbour
-    for (const polyPatch& pp : patches)
-    {
-        if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
-        {
-            label facei = pp.start();
-            label bFacei = pp.start()-mesh.nInternalFaces();
-
-            forAll(pp, i)
-            {
-                const label own = agglom[faceOwner[facei]];
-                const label globalNei = globalNeighbour[bFacei];
-
-                if
-                (
-                   !globalAgglom.isLocal(globalNei)
-                 || globalAgglom.toLocal(globalNei) != own
-                )
-                {
-                    const label ownIndex = offsets[own] + nFacesPerCell[own]++;
-                    m[ownIndex] = globalNei;
-                    w[ownIndex] = mag(mesh.faceAreas()[facei]);
-                }
-
-                ++facei;
-                ++bFacei;
-            }
-        }
-    }
-
-
-    // Check for duplicates connections between cells
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Done as postprocessing step since we now have cellCells.
-
-    if (cellCells.size() == 0)
-    {
-        return;
-    }
-
-    label newIndex = 0;
-    labelHashSet nbrCells;
-
-    label startIndex = cellCells.offsets()[0];
-
-    forAll(cellCells, celli)
-    {
-        nbrCells.clear();
-        nbrCells.insert(globalAgglom.toGlobal(celli));
-
-        const label endIndex = cellCells.offsets()[celli+1];
-
-        for (label i = startIndex; i < endIndex; ++i)
-        {
-            if (nbrCells.insert(cellCells.m()[i]))
-            {
-                cellCells.m()[newIndex] = cellCells.m()[i];
-                cellCellWeights.m()[newIndex] = cellCellWeights.m()[i];
-                newIndex++;
-            }
-        }
-        startIndex = endIndex;
-        cellCells.offsets()[celli+1] = newIndex;
-        cellCellWeights.offsets()[celli+1] = newIndex;
-    }
-
-    cellCells.m().setSize(newIndex);
-    cellCellWeights.m().setSize(newIndex);
+    globalMeshData::calcCellCells
+    (
+        mesh,
+        agglom,
+        nLocalCoarse,
+        parallel,
+        cellCells,
+        cellCellWeights
+    );
 }
 
 
@@ -882,11 +500,11 @@ Foam::labelList Foam::decompositionMethod::decompose
     // Any weights specified?
     const bool hasWeights = returnReduceOr(cellWeights.size());
 
-    if (hasWeights && cellWeights.size() != mesh.nCells())
+    if (hasWeights && (cellWeights.size() != mesh.nCells()))
     {
         FatalErrorInFunction
-            << "Number of weights " << cellWeights.size()
-            << " differs from number of cells " << mesh.nCells()
+            << "Number of weights (" << cellWeights.size()
+            << ") != number of cells (" << mesh.nCells() << ")"
             << exit(FatalError);
     }
 
@@ -1219,7 +837,7 @@ void Foam::decompositionMethod::setConstraints
     List<labelPair>& explicitConnections
 ) const
 {
-    blockedFace.setSize(mesh.nFaces());
+    blockedFace.resize_nocopy(mesh.nFaces());
     blockedFace = true;
 
     specifiedProcessorFaces.clear();
@@ -1329,14 +947,4 @@ Foam::labelList Foam::decompositionMethod::decompose
 }
 
 
-Foam::labelList Foam::decompositionMethod::decompose
-(
-    const pointField& points
-) const
-{
-    NotImplemented;
-    return labelList();
-}
-
-
 // ************************************************************************* //
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H
index 343943fb24..6d0ac936f1 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2021 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef decompositionMethod_H
-#define decompositionMethod_H
+#ifndef Foam_decompositionMethod_H
+#define Foam_decompositionMethod_H
 
 #include "polyMesh.H"
 #include "CompactListList.H"
@@ -63,12 +63,6 @@ class decompositionMethod
         //- Set PtrList of constraints by reading decompDict_.
         void readConstraints();
 
-        //- No copy construct
-        decompositionMethod(const decompositionMethod&) = delete;
-
-        //- No copy assignment
-        void operator=(const decompositionMethod&) = delete;
-
 
 protected:
 
@@ -84,7 +78,7 @@ protected:
         };
 
 
-    // Protected data
+    // Protected Data
 
         //- Top-level decomposition dictionary (eg, decomposeParDict)
         const dictionary& decompDict_;
@@ -137,8 +131,24 @@ protected:
         ) const;
 
 
+    // Constructors
+
+        //- Construct with specified number of domains,
+        //- no coefficients or constraints
+        explicit decompositionMethod(const label numDomains);
+
+
 public:
 
+    // Generated Methods
+
+        //- No copy construct
+        decompositionMethod(const decompositionMethod&) = delete;
+
+        //- No copy assignment
+        void operator=(const decompositionMethod&) = delete;
+
+
     //- Runtime type information
     TypeName("decompositionMethod");
 
@@ -207,51 +217,49 @@ public:
     // Member Functions
 
         //- Number of domains
-        inline label nDomains() const noexcept
+        label nDomains() const noexcept
         {
             return nDomains_;
         }
 
+        //- True if the method is purely geometric,
+        //- often using cell centre points
+        virtual bool geometric() const { return false; }
+
         //- Is method parallel aware?
         //  (i.e. does it synchronize domains across proc boundaries)
         virtual bool parallelAware() const = 0;
 
+        // //- Is internally method parallel aware
+        // virtual bool parallelNative() const  { return false; }
+
 
     // No topology (implemented by geometric decomposers)
 
-        //- Return the wanted processor number for every coordinate.
+        //- Return the wanted processor number for every coordinate,
+        //- using uniform or specified point weights.
         virtual labelList decompose
         (
             const pointField& points,
-            const scalarField& pointWeights
+            const scalarField& pointWeights = scalarField::null()
         ) const;
 
-        //- Decompose with uniform weights on the points
-        virtual labelList decompose(const pointField& points) const;
-
 
     // Topology provided by mesh
 
-        //- Return for every coordinate the wanted processor number.
+        //- Return for every coordinate the wanted processor number,
+        //- using uniform or specified point weights.
         //  Use the mesh connectivity (if needed)
         virtual labelList decompose
         (
             const polyMesh& mesh,
             const pointField& points,
-            const scalarField& pointWeights
+            const scalarField& pointWeights = scalarField::null()
         ) const = 0;
 
-        //- Decompose with uniform weights on the points
-        virtual labelList decompose
-        (
-            const polyMesh& mesh,
-            const pointField& points
-        ) const;
-
-
-        //- Return for every coordinate the wanted processor number. Gets
-        //  passed agglomeration map (from fine to coarse cells) and coarse
-        //  cell
+        //- Return for every coordinate the wanted processor number.
+        //  Gets passed agglomeration map (from fine to coarse cells)
+        //  and coarse cell
         //  location. Can be overridden by decomposers that provide this
         //  functionality natively. Coarse cells are local to the processor
         //  (if in parallel). If you want to have coarse cells spanning
@@ -261,19 +269,25 @@ public:
             const polyMesh& mesh,
             const labelList& cellToRegion,
             const pointField& regionPoints,
-            const scalarField& regionWeights
+            const scalarField& regionWeights = scalarField::null()
         ) const;
 
-        //- Like decompose but with uniform weights on the regions
+
+    // Topology provided explicitly
+
+        //- Return for every coordinate the wanted processor number.
+        //  The connectivity is equal to mesh.cellCells() except for
+        //  - in parallel the cell numbers are global cell numbers
+        //    (starting
+        //    from 0 at processor0 and then incrementing all through the
+        //    processors)
+        //  - the connections are across coupled patches
         virtual labelList decompose
         (
-            const polyMesh& mesh,
-            const labelList& cellToRegion,
-            const pointField& regionPoints
-        ) const;
-
-
-    // Topology provided explicitly addressing
+            const CompactListList<label>& globalCellCells,
+            const pointField& cc,
+            const scalarField& cWeights = scalarField::null()
+        ) const = 0;
 
         //- Return for every coordinate the wanted processor number.
         //  The connectivity is equal to mesh.cellCells() except for
@@ -286,47 +300,12 @@ public:
         (
             const labelListList& globalCellCells,
             const pointField& cc,
-            const scalarField& cWeights
+            const scalarField& cWeights = scalarField::null()
         ) const = 0;
 
-        //- Like decompose but with uniform weights on the cells
-        virtual labelList decompose
-        (
-            const labelListList& globalCellCells,
-            const pointField& cc
-        ) const;
-
 
     // Other
 
-        //- Helper: determine (local or global) cellCells from mesh
-        //  agglomeration. Agglomeration is local to the processor.
-        //  local  : connections are in local indices. Coupled across
-        //           cyclics but not processor patches.
-        //  global : connections are in global indices. Coupled across
-        //            cyclics and processor patches.
-        static void calcCellCells
-        (
-            const polyMesh& mesh,
-            const labelList& agglom,
-            const label nLocalCoarse,
-            const bool global,
-            CompactListList<label>& cellCells
-        );
-
-        //- Helper: determine (local or global) cellCells and face weights
-        //  from mesh agglomeration.
-        //  Uses mag of faceArea as weights
-        static void calcCellCells
-        (
-            const polyMesh& mesh,
-            const labelList& agglom,
-            const label nLocalCoarse,
-            const bool parallel,
-            CompactListList<label>& cellCells,
-            CompactListList<scalar>& cellCellWeights
-        );
-
         //- Helper: extract constraints:
         //  blockedface: existing faces where owner and neighbour on same
         //               proc
@@ -390,6 +369,39 @@ public:
             const polyMesh& mesh,
             const scalarField& cWeights
         ) const;
+
+
+    // Housekeeping
+
+        //- Determine (local or global) cellCells from mesh agglomeration.
+        //  Agglomeration is local to the processor.
+        //  local  : connections are in local indices. Coupled across
+        //           cyclics but not processor patches.
+        //  global : connections are in global indices. Coupled across
+        //            cyclics and processor patches.
+        FOAM_DEPRECATED_STRICT(2023-11, "globalMeshData::calcCellCells()")
+        static void calcCellCells
+        (
+            const polyMesh& mesh,
+            const labelList& agglom,
+            const label nLocalCoarse,
+            const bool parallel,  //!< Use global cell ids in parallel
+            CompactListList<label>& cellCells
+        );
+
+        //- Determine (local or global) cellCells and face weights
+        //- from mesh agglomeration.
+        //  Uses mag of faceArea as weights
+        FOAM_DEPRECATED_STRICT(2023-11, "globalMeshData::calcCellCells()")
+        static void calcCellCells
+        (
+            const polyMesh& mesh,
+            const labelList& agglom,
+            const label nLocalCoarse,
+            const bool parallel,  //!< Use global cell ids in parallel
+            CompactListList<label>& cellCells,
+            CompactListList<scalar>& cellCellWeights
+        );
 };
 
 
diff --git a/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.C b/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.C
index 9c69a595e4..5b01f1f74e 100644
--- a/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -149,6 +149,17 @@ void Foam::geomDecomp::checkDecompositionDirections
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::geomDecomp::geomDecomp(const Vector<label>& divisions)
+:
+    decompositionMethod(divisions.x()*divisions.y()*divisions.z()),
+    delta_(0.001),
+    csys_(),
+    n_(divisions),
+    order_(0,1,2),
+    coeffsDict_(dictionary::null)
+{}
+
+
 Foam::geomDecomp::geomDecomp
 (
     const word& derivedType,
diff --git a/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.H b/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.H
index 49661c9614..6eb5cc1305 100644
--- a/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,8 +44,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef geomDecomp_H
-#define geomDecomp_H
+#ifndef Foam_geomDecomp_H
+#define Foam_geomDecomp_H
 
 #include "decompositionMethod.H"
 #include "cartesianCS.H"
@@ -105,11 +105,14 @@ protected:
         //- Check that mesh directions are compatible with decomposition
         void checkDecompositionDirections(const Vector<label>&) const;
 
-
 public:
 
     // Constructors
 
+        //- Construct with number of x/y/z division
+        //- (no coefficients or constraints)
+        explicit geomDecomp(const Vector<label>& divisions);
+
         //- Construct for derived type name and decomposition dictionary
         geomDecomp
         (
@@ -131,15 +134,15 @@ public:
 
     // Member Functions
 
+        //- Purely geometric methods
+        virtual bool geometric() const { return true; }
+
         //- Return for every coordinate the wanted processor number.
         virtual labelList decompose
         (
             const pointField& points,
-            const scalarField& pointWeights
+            const scalarField& pointWeights = scalarField::null()
         ) const = 0;
-
-        //- Decompose with uniform weights on the points
-        virtual labelList decompose(const pointField& points) const = 0;
 };
 
 
diff --git a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
index f62efd1659..fc6b72333e 100644
--- a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -674,6 +674,12 @@ Foam::label Foam::hierarchGeomDecomp::sortComponent
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::hierarchGeomDecomp::hierarchGeomDecomp(const Vector<label>& divisions)
+:
+    geomDecomp(divisions)
+{}
+
+
 Foam::hierarchGeomDecomp::hierarchGeomDecomp
 (
     const dictionary& decompDict,
@@ -686,53 +692,15 @@ Foam::hierarchGeomDecomp::hierarchGeomDecomp
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::labelList Foam::hierarchGeomDecomp::decompose
-(
-    const pointField& points
-) const
-{
-    // construct a list for the final result
-    labelList finalDecomp(points.size(), Zero);
-
-    // Start off with every point sorted onto itself.
-    labelList slice(identity(points.size()));
-
-    const pointField rotatedPoints(adjustPoints(points));
-
-    // Calculate tolerance of cell distribution. For large cases finding
-    // distribution to the cell exact would cause too many iterations so allow
-    // some slack.
-    const label allSize = returnReduce(points.size(), sumOp<label>());
-    const label sizeTol = max(1, label(1e-3*allSize/nDomains_));
-
-    // Sort recursive
-    const label nWarnings = sortComponent
-    (
-        sizeTol,
-        rotatedPoints,
-        slice,
-        0,              // Sort first component in order_
-        1,              // Offset for different x bins.
-        finalDecomp
-    );
-
-    if (nWarnings)
-    {
-        WarningInFunction
-            << "\nEncountered " << nWarnings << " occurrences where the desired"
-               " decomposition split could not be properly satisfied" << endl;
-    }
-
-    return finalDecomp;
-}
-
-
 Foam::labelList Foam::hierarchGeomDecomp::decompose
 (
     const pointField& points,
     const scalarField& weights
 ) const
 {
+    // Uniform weighting if weights are empty or poorly sized
+    const bool hasWeights = returnReduceAnd(points.size() == weights.size());
+
     // Construct a list for the final result
     labelList finalDecomp(points.size(), Zero);
 
@@ -748,16 +716,32 @@ Foam::labelList Foam::hierarchGeomDecomp::decompose
     const label sizeTol = max(1, label(1e-3*allSize/nDomains_));
 
     // Sort recursive
-    const label nWarnings = sortComponent
-    (
-        sizeTol,
-        weights,
-        rotatedPoints,
-        slice,
-        0,              // Sort first component in order_
-        1,              // Offset for different x bins.
-        finalDecomp
-    );
+    label nWarnings = 0;
+    if (hasWeights)
+    {
+        nWarnings = sortComponent
+        (
+            sizeTol,
+            weights,
+            rotatedPoints,
+            slice,
+            0,              // Sort first component in order_
+            1,              // Offset for different x bins.
+            finalDecomp
+        );
+    }
+    else
+    {
+        nWarnings = sortComponent
+        (
+            sizeTol,
+            rotatedPoints,
+            slice,
+            0,              // Sort first component in order_
+            1,              // Offset for different x bins.
+            finalDecomp
+        );
+    }
 
     if (nWarnings)
     {
diff --git a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
index 479585d022..aadd901634 100644
--- a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -62,8 +62,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef hierarchGeomDecomp_H
-#define hierarchGeomDecomp_H
+#ifndef Foam_hierarchGeomDecomp_H
+#define Foam_hierarchGeomDecomp_H
 
 #include "geomDecomp.H"
 
@@ -169,12 +169,17 @@ public:
     //- No copy assignment
     void operator=(const hierarchGeomDecomp&) = delete;
 
+
     //- Runtime type information
     TypeName("hierarchical");
 
 
     // Constructors
 
+        //- Construct with number of x/y/z division
+        //- (no coefficients or constraints)
+        explicit hierarchGeomDecomp(const Vector<label>& divisions);
+
         //- Construct given decomposition dictionary and optional region name
         explicit hierarchGeomDecomp
         (
@@ -196,42 +201,36 @@ public:
         }
 
         //- Return for every coordinate the wanted processor number.
+        //- using uniform or specified point weights.
         virtual labelList decompose
         (
-            const pointField&,
-            const scalarField& weights
+            const pointField& points,
+            const scalarField& weights = scalarField::null()
         ) const;
 
-        //- Decompose with uniform weights.
-        //  Code for weighted decomposition is a bit complex,
-        //  so kept separate for now.
-        virtual labelList decompose(const pointField&) const;
-
-
         //- Return for every coordinate the wanted processor number.
         //  Use the mesh connectivity (if needed).
+        //  With uniform or specified point weights.
         virtual labelList decompose
         (
             const polyMesh& mesh,
             const pointField& cc,
-            const scalarField& cWeights
+            const scalarField& cWeights = scalarField::null()
         ) const
         {
             checkDecompositionDirections(mesh.geometricD());
             return decompose(cc, cWeights);
         }
 
-        //- Decompose with uniform weights.
-        //  Code for weighted decomposition is a bit complex,
-        //  so kept separate for now.
+        //- Explicitly provided connectivity
         virtual labelList decompose
         (
-            const polyMesh& mesh,
-            const pointField& cc
+            const CompactListList<label>& globalCellCells,  //!< unused
+            const pointField& cc,
+            const scalarField& cWeights = scalarField::null()
         ) const
         {
-            checkDecompositionDirections(mesh.geometricD());
-            return decompose(cc);
+            return decompose(cc, cWeights);
         }
 
         //- Return for every coordinate the wanted processor number.
@@ -243,22 +242,13 @@ public:
         //  - the connections are across coupled patches
         virtual labelList decompose
         (
-            const labelListList& globalCellCells,  // unused
+            const labelListList& globalCellCells,  //!< unused
             const pointField& cc,
-            const scalarField& cWeights
+            const scalarField& cWeights = scalarField::null()
         ) const
         {
             return decompose(cc, cWeights);
         }
-
-        virtual labelList decompose
-        (
-            const labelListList& globalCellCells,  // unused
-            const pointField& cc
-        ) const
-        {
-            return decompose(cc);
-        }
 };
 
 
diff --git a/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.C b/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.C
index 0c02723cd5..a88cfd5c02 100644
--- a/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -78,7 +78,9 @@ Foam::labelList Foam::manualDecomp::decompose
                 dataFile_,
                 mesh.facesInstance(),
                 mesh.thisDb(),
-                IOobject::MUST_READ
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
             )
         )
     );
diff --git a/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.H b/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.H
index b66bf74d61..54b6804267 100644
--- a/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,8 +41,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef manualDecomp_H
-#define manualDecomp_H
+#ifndef Foam_manualDecomp_H
+#define Foam_manualDecomp_H
 
 #include "decompositionMethod.H"
 
@@ -61,8 +61,9 @@ class manualDecomp
 
         fileName dataFile_;
 
+public:
 
-    // Private Member Functions
+    // Generated Methods
 
         //- No copy construct
         manualDecomp(const manualDecomp&) = delete;
@@ -71,8 +72,6 @@ class manualDecomp
         void operator=(const manualDecomp&) = delete;
 
 
-public:
-
     //- Runtime type information
     TypeName("manual");
 
@@ -106,16 +105,29 @@ public:
         (
             const polyMesh& mesh,
             const pointField& cc,
-            const scalarField& cWeights
+            const scalarField& cWeights = scalarField::null()
         ) const;
 
+        //- Return for every coordinate the wanted processor number.
+        //  Explicitly provided connectivity - does not use mesh_.
+        virtual labelList decompose
+        (
+            const CompactListList<label>& globalCellCells,
+            const pointField& cc,
+            const scalarField& cWeights = scalarField::null()
+        ) const
+        {
+            NotImplemented;
+            return labelList();
+        }
+
         //- Return for every coordinate the wanted processor number.
         //  Explicitly provided connectivity - does not use mesh_.
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
             const pointField& cc,
-            const scalarField& cWeights
+            const scalarField& cWeights = scalarField::null()
         ) const
         {
             NotImplemented;
diff --git a/src/parallel/decompose/decompositionMethods/metisLikeDecomp/metisLikeDecomp.C b/src/parallel/decompose/decompositionMethods/metisLikeDecomp/metisLikeDecomp.C
index 08953fd02a..210dcde5f3 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-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +28,7 @@ License
 #include "metisLikeDecomp.H"
 #include "Time.H"
 #include "globalIndex.H"
+#include "globalMeshData.H"
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
@@ -39,8 +40,9 @@ Foam::label Foam::metisLikeDecomp::decomposeGeneral
     labelList& decomp
 ) const
 {
-    if (!Pstream::parRun())
+    if (!UPstream::parRun())
     {
+        // Treat zero-sized weights as uniform weighting
         return decomposeSerial
         (
             adjncy,
@@ -75,10 +77,10 @@ Foam::label Foam::metisLikeDecomp::decomposeGeneral
 
 
     List<label> allXadj;
-    if (Pstream::master())
+    if (UPstream::master())
     {
         allXadj.resize(globalCells.totalSize()+1);
-        allXadj.last() = globalAdjncy.totalSize();  // Final end offset
+        allXadj.back() = globalAdjncy.totalSize();  // Final end offset
 
         // My values - no renumbering required
         SubList<label>(allXadj, globalCells.localSize(0)) =
@@ -137,7 +139,7 @@ Foam::label Foam::metisLikeDecomp::decomposeGeneral
     }
 
     // Local to global renumbering
-    if (Pstream::master())
+    if (UPstream::master())
     {
         for (const int proci : globalCells.subProcs())
         {
@@ -147,10 +149,12 @@ Foam::label Foam::metisLikeDecomp::decomposeGeneral
         }
     }
 
-    // Ignore zero-sized weights ... and poorly sized ones too
+    // Uniform weighting if weights are empty or poorly sized
+
     List<scalar> allWeights;
     if (returnReduceAnd(cWeights.size() == globalCells.localSize()))
     {
+        // ie, hasWeights
         allWeights = globalCells.gather(cWeights);
     }
 
@@ -158,7 +162,7 @@ Foam::label Foam::metisLikeDecomp::decomposeGeneral
     // Global decomposition
     labelList allDecomp;
 
-    if (Pstream::master())
+    if (UPstream::master())
     {
         decomposeSerial
         (
@@ -183,6 +187,13 @@ Foam::label Foam::metisLikeDecomp::decomposeGeneral
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::metisLikeDecomp::metisLikeDecomp(const label numDomains)
+:
+    decompositionMethod(numDomains),
+    coeffsDict_(dictionary::null)
+{}
+
+
 Foam::metisLikeDecomp::metisLikeDecomp
 (
     const word& derivedType,
@@ -205,18 +216,16 @@ Foam::labelList Foam::metisLikeDecomp::decompose
     const scalarField& pointWeights
 ) const
 {
-    if (points.size() != mesh.nCells())
+    if (!points.empty() && (points.size() != mesh.nCells()))
     {
         FatalErrorInFunction
-            << "Can only use this decomposition method for entire mesh" << nl
-            << "and supply one coordinate (cellCentre) for every cell." << nl
-            << "The number of coordinates " << points.size() << nl
-            << "The number of cells in the mesh " << mesh.nCells() << nl
+            << "Number of cell centres (" << points.size()
+            << ") != number of cells (" << mesh.nCells() << ")"
             << exit(FatalError);
     }
 
     CompactListList<label> cellCells;
-    calcCellCells
+    globalMeshData::calcCellCells
     (
         mesh,
         identity(mesh.nCells()),
@@ -250,8 +259,8 @@ Foam::labelList Foam::metisLikeDecomp::decompose
     if (agglom.size() != mesh.nCells())
     {
         FatalErrorInFunction
-            << "Size of cell-to-coarse map " << agglom.size()
-            << " differs from number of cells in mesh " << mesh.nCells()
+            << "Agglomeration size (" << agglom.size()
+            << ") != number of cells (" << mesh.nCells() << ")"
             << exit(FatalError);
     }
 
@@ -260,7 +269,7 @@ Foam::labelList Foam::metisLikeDecomp::decompose
     //   xadj(celli) : start of information in adjncy for celli
 
     CompactListList<label> cellCells;
-    calcCellCells
+    globalMeshData::calcCellCells
     (
         mesh,
         agglom,
@@ -279,16 +288,42 @@ Foam::labelList Foam::metisLikeDecomp::decompose
         decomp
     );
 
+    // From coarse back to fine for original mesh
+    return labelList(decomp, agglom);
+}
 
-    // Rework back into decomposition for original mesh
-    labelList fineDistribution(agglom.size());
 
-    forAll(fineDistribution, i)
+Foam::labelList Foam::metisLikeDecomp::decompose
+(
+    const CompactListList<label>& globalCellCells,
+    const pointField& cellCentres,
+    const scalarField& cellWeights
+) const
+{
+    if (!cellCentres.empty() && (cellCentres.size() != globalCellCells.size()))
     {
-        fineDistribution[i] = decomp[agglom[i]];
+        FatalErrorInFunction
+            << "Number of cell centres (" << cellCentres.size()
+            << ") != number of cells (" << globalCellCells.size() << ")"
+            << exit(FatalError);
     }
 
-    return fineDistribution;
+    // CompactListList is already
+    // Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+
+    // Decompose using default weights
+    labelList decomp;
+    decomposeGeneral
+    (
+        globalCellCells.values(),
+        globalCellCells.offsets(),
+        cellWeights,
+        decomp
+    );
+
+    return decomp;
 }
 
 
@@ -299,12 +334,12 @@ Foam::labelList Foam::metisLikeDecomp::decompose
     const scalarField& cellWeights
 ) const
 {
-    if (cellCentres.size() != globalCellCells.size())
+    if (!cellCentres.empty() && (cellCentres.size() != globalCellCells.size()))
     {
         FatalErrorInFunction
-            << "Inconsistent number of cells (" << globalCellCells.size()
-            << ") and number of cell centres (" << cellCentres.size()
-            << ")." << exit(FatalError);
+            << "Number of cell centres (" << cellCentres.size()
+            << ") != number of cells (" << globalCellCells.size() << ")"
+            << exit(FatalError);
     }
 
     // Make Metis CSR (Compressed Storage Format) storage
diff --git a/src/parallel/decompose/decompositionMethods/metisLikeDecomp/metisLikeDecomp.H b/src/parallel/decompose/decompositionMethods/metisLikeDecomp/metisLikeDecomp.H
index 2960f19207..595d23fa51 100644
--- a/src/parallel/decompose/decompositionMethods/metisLikeDecomp/metisLikeDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/metisLikeDecomp/metisLikeDecomp.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,8 +37,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef metisLikeDecomp_H
-#define metisLikeDecomp_H
+#ifndef Foam_metisLikeDecomp_H
+#define Foam_metisLikeDecomp_H
 
 #include "decompositionMethod.H"
 
@@ -91,6 +91,9 @@ protected:
 
     // Constructors
 
+        //- Construct with number of domains (no coefficients or constraints)
+        explicit metisLikeDecomp(const label numDomains);
+
         //- Construct for derived type name and decomposition dictionary.
         metisLikeDecomp
         (
@@ -134,8 +137,8 @@ public:
         virtual labelList decompose
         (
             const polyMesh& mesh,
-            const pointField& points,
-            const scalarField& pointWeights
+            const pointField& points = pointField::null(),
+            const scalarField& pointWeights = scalarField::null()
         ) const;
 
         //- Return for every coordinate the wanted processor number.
@@ -147,8 +150,23 @@ public:
         (
             const polyMesh& mesh,
             const labelList& agglom,
-            const pointField& regionPoints,
-            const scalarField& regionWeights
+            const pointField& agglomPoints,
+            const scalarField& agglomWeights = scalarField::null()
+        ) const;
+
+        //- Return for every coordinate the wanted processor number.
+        //  Explicitly provided mesh connectivity.
+        //  The connectivity is equal to mesh.cellCells() except for
+        //  - in parallel the cell numbers are global cell numbers (starting
+        //    from 0 at processor0 and then incrementing all through the
+        //    processors)
+        //  - the connections are across coupled patches
+        //  See note on weights above.
+        virtual labelList decompose
+        (
+            const CompactListList<label>& globalCellCells,
+            const pointField& cellCentres = pointField::null(),
+            const scalarField& cellWeights = scalarField::null()
         ) const;
 
         //- Return for every coordinate the wanted processor number.
@@ -162,8 +180,8 @@ public:
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
-            const pointField& cellCentres,
-            const scalarField& cellWeights
+            const pointField& cellCentres = pointField::null(),
+            const scalarField& cellWeights = scalarField::null()
         ) const;
 };
 
diff --git a/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C b/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C
index f291df37ff..d3d0f3469b 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-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,6 +30,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "IFstream.H"
 #include "globalIndex.H"
+#include "globalMeshData.H"
 #include "mapDistribute.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -194,8 +195,7 @@ void Foam::multiLevelDecomp::setMethods()
 
     label nLevels = 0;
 
-    methods_.clear();
-    methods_.setSize(methodsDict_.size());
+    methods_.resize_null(methodsDict_.size());
     for (const entry& dEntry : methodsDict_)
     {
         // Dictionary entries only
@@ -211,7 +211,7 @@ void Foam::multiLevelDecomp::setMethods()
         }
     }
 
-    methods_.setSize(nLevels);
+    methods_.resize(nLevels);
 
     // Verify that nTotal is correct based on what each method delivers
 
@@ -280,7 +280,7 @@ void Foam::multiLevelDecomp::subsetGlobalCellCells
     // Now subCellCells contains indices into oldToNew which are the
     // new locations of the neighbouring cells.
 
-    cutConnections.setSize(nDomains);
+    cutConnections.resize_nocopy(nDomains);
     cutConnections = 0;
 
     forAll(subCellCells, subCelli)
@@ -399,7 +399,12 @@ void Foam::multiLevelDecomp::decompose
 
             // Subset point-wise data.
             pointField subPoints(points, domainPoints);
-            scalarField subWeights(pointWeights, domainPoints);
+            scalarField subWeights;
+            if (pointWeights.size() == points.size())
+            {
+                subWeights = scalarField(pointWeights, domainPoints);
+            }
+
             labelList subPointMap(labelUIndList(pointMap, domainPoints));
             // Subset point-point addressing (adapt global numbering)
             labelListList subPointPoints;
@@ -607,9 +612,16 @@ Foam::labelList Foam::multiLevelDecomp::decompose
 ) const
 {
     CompactListList<label> cellCells;
-    calcCellCells(mesh, identity(cc.size()), cc.size(), true, cellCells);
+    globalMeshData::calcCellCells
+    (
+        mesh,
+        identity(cc.size()),
+        cc.size(),
+        true,
+        cellCells
+    );
 
-    labelList finalDecomp(cc.size(), Zero);
+    labelList finalDecomp(cc.size(), Foam::zero{});
     labelList cellMap(identity(cc.size()));
 
     decompose
@@ -617,9 +629,9 @@ Foam::labelList Foam::multiLevelDecomp::decompose
         cellCells.unpack(),
         cc,
         cWeights,
-        cellMap,      // map back to original cells
-        0,
-        0,
+        cellMap,      // map back to original cell points
+        0,            // currLevel
+        0,            // leafOffset
 
         finalDecomp
     );
@@ -630,22 +642,48 @@ Foam::labelList Foam::multiLevelDecomp::decompose
 
 Foam::labelList Foam::multiLevelDecomp::decompose
 (
-    const labelListList& globalPointPoints,
-    const pointField& points,
-    const scalarField& pointWeights
+    const CompactListList<label>& globalCellCells,
+    const pointField& cc,
+    const scalarField& cWeights
 ) const
 {
-    labelList finalDecomp(points.size(), Zero);
-    labelList pointMap(identity(points.size()));
+    labelList finalDecomp(cc.size(), Foam::zero{});
+    labelList cellMap(identity(cc.size()));
 
     decompose
     (
-        globalPointPoints,
-        points,
-        pointWeights,
-        pointMap,       // map back to original points
-        0,
-        0,
+        globalCellCells.unpack(),
+        cc,
+        cWeights,
+        cellMap,      // map back to original cell points
+        0,            // currLevel
+        0,            // leafOffset
+
+        finalDecomp
+    );
+
+    return finalDecomp;
+}
+
+
+Foam::labelList Foam::multiLevelDecomp::decompose
+(
+    const labelListList& globalCellCells,
+    const pointField& cc,
+    const scalarField& cWeights
+) const
+{
+    labelList finalDecomp(cc.size(), Foam::zero{});
+    labelList cellMap(identity(cc.size()));
+
+    decompose
+    (
+        globalCellCells,
+        cc,
+        cWeights,
+        cellMap,      // map back to original cell points
+        0,            // currLevel
+        0,            // leafOffset
 
         finalDecomp
     );
diff --git a/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.H b/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.H
index 866a143657..468af9c45d 100644
--- a/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef multiLevelDecomp_H
-#define multiLevelDecomp_H
+#ifndef Foam_multiLevelDecomp_H
+#define Foam_multiLevelDecomp_H
 
 #include "decompositionMethod.H"
 
@@ -99,6 +99,10 @@ class multiLevelDecomp
         ) const;
 
 
+public:
+
+    // Generated Methods
+
         //- No copy construct
         multiLevelDecomp(const multiLevelDecomp&) = delete;
 
@@ -106,8 +110,6 @@ class multiLevelDecomp
         void operator=(const multiLevelDecomp&) = delete;
 
 
-public:
-
     //- Runtime type information
     TypeName("multiLevel");
 
@@ -128,8 +130,7 @@ public:
 
     // Member Functions
 
-        //- Is method parallel aware?
-        //  i.e. does it synchronize domains across proc boundaries
+        //- Is parallel aware when all sub-methods are also parallel-aware
         virtual bool parallelAware() const;
 
         //- Inherit decompose from decompositionMethod
@@ -141,7 +142,16 @@ public:
         (
             const polyMesh& mesh,
             const pointField& points,
-            const scalarField& pointWeights
+            const scalarField& pointWeights = scalarField::null()
+        ) const;
+
+        //- Return for every coordinate the wanted processor number.
+        //  Explicitly provided connectivity - does not use mesh_.
+        virtual labelList decompose
+        (
+            const CompactListList<label>& globalCellCells,
+            const pointField& cc,
+            const scalarField& cWeights = scalarField::null()
         ) const;
 
         //- Return for every coordinate the wanted processor number.
@@ -150,7 +160,7 @@ public:
         (
             const labelListList& globalCellCells,
             const pointField& cc,
-            const scalarField& cWeights
+            const scalarField& cWeights = scalarField::null()
         ) const;
 };
 
diff --git a/src/parallel/decompose/decompositionMethods/noDecomp/noDecomp.C b/src/parallel/decompose/decompositionMethods/noDecomp/noDecomp.C
index 5a8d09c86c..ab1b4fa0d3 100644
--- a/src/parallel/decompose/decompositionMethods/noDecomp/noDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/noDecomp/noDecomp.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,6 +46,12 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::noDecomp::noDecomp(const label numDomains)
+:
+    decompositionMethod(numDomains)
+{}
+
+
 Foam::noDecomp::noDecomp
 (
     const dictionary& decompDict,
@@ -56,4 +62,49 @@ Foam::noDecomp::noDecomp
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::noDecomp::decompose
+(
+    const pointField& points,
+    const scalarField&
+) const
+{
+    return labelList(points.size(), UPstream::myProcNo());
+}
+
+
+Foam::labelList Foam::noDecomp::decompose
+(
+    const polyMesh& mesh,
+    const pointField&,
+    const scalarField&
+) const
+{
+    return labelList(mesh.nCells(), UPstream::myProcNo());
+}
+
+
+Foam::labelList Foam::noDecomp::decompose
+(
+    const CompactListList<label>& globalCellCells,
+    const pointField&,
+    const scalarField&
+) const
+{
+    return labelList(globalCellCells.size(), UPstream::myProcNo());
+}
+
+
+Foam::labelList Foam::noDecomp::decompose
+(
+    const labelListList& globalCellCells,
+    const pointField&,
+    const scalarField&
+) const
+{
+    return labelList(globalCellCells.size(), UPstream::myProcNo());
+}
+
+
 // ************************************************************************* //
diff --git a/src/parallel/decompose/decompositionMethods/noDecomp/noDecomp.H b/src/parallel/decompose/decompositionMethods/noDecomp/noDecomp.H
index 8c1b08c9de..70d0980093 100644
--- a/src/parallel/decompose/decompositionMethods/noDecomp/noDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/noDecomp/noDecomp.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2015 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,8 +37,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef noDecomp_H
-#define noDecomp_H
+#ifndef Foam_noDecomp_H
+#define Foam_noDecomp_H
 
 #include "decompositionMethod.H"
 
@@ -53,7 +53,9 @@ class noDecomp
 :
     public decompositionMethod
 {
-    // Private Member Functions
+public:
+
+    // Generated Methods
 
         //- No copy construct
         noDecomp(const noDecomp&) = delete;
@@ -62,14 +64,15 @@ class noDecomp
         void operator=(const noDecomp&) = delete;
 
 
-public:
-
     //- Runtime type information
     TypeNameNoDebug("none");
 
 
     // Constructors
 
+        //- Construct with number of domains (ignored)
+        explicit noDecomp(const label numDomains);
+
         //- Construct given decomposition dictionary and optional region name
         explicit noDecomp
         (
@@ -91,33 +94,35 @@ public:
         }
 
         //- Return for every coordinate the wanted processor number.
-        //  Use the mesh connectivity (if needed)
+        virtual labelList decompose
+        (
+            const pointField& points,
+            const scalarField& pointWeights_unused = scalarField::null()
+        ) const;
+
+        //- Return for every cell the current local processor rank.
         virtual labelList decompose
         (
             const polyMesh& mesh,
-            const pointField& cc,
-            const scalarField& cWeights
-        ) const
-        {
-            return labelList(cc.size(), Pstream::myProcNo());
-        }
+            const pointField& cc_unused = pointField::null(),
+            const scalarField& cWeights_unused = scalarField::null()
+        ) const;
 
-        //- Return for every coordinate the wanted processor number.
-        //  Explicitly provided connectivity - does not use mesh_.
-        //  The connectivity is equal to mesh.cellCells() except for
-        //  - in parallel the cell numbers are global cell numbers (starting
-        //    from 0 at processor0 and then incrementing all through the
-        //    processors)
-        //  - the connections are across coupled patches
+        //- Return for every cell the current local processor rank.
+        virtual labelList decompose
+        (
+            const CompactListList<label>& globalCellCells,
+            const pointField& cc_unused = pointField::null(),
+            const scalarField& cWeights_unused = scalarField::null()
+        ) const;
+
+        //- Return for every cell the current local processor rank.
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
-            const pointField& cc,
-            const scalarField& cWeights
-        ) const
-        {
-            return labelList(globalCellCells.size(), Pstream::myProcNo());
-        }
+            const pointField& cc_unused = pointField::null(),
+            const scalarField& cWeights_unused = scalarField::null()
+        ) const;
 };
 
 
diff --git a/src/parallel/decompose/decompositionMethods/randomDecomp/randomDecomp.C b/src/parallel/decompose/decompositionMethods/randomDecomp/randomDecomp.C
index 625d2cbcb7..a9985fd4f3 100644
--- a/src/parallel/decompose/decompositionMethods/randomDecomp/randomDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/randomDecomp/randomDecomp.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -82,6 +82,13 @@ Foam::labelList Foam::randomDecomp::randomMap(const label nCells) const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::randomDecomp::randomDecomp(const label numDomains)
+:
+    decompositionMethod(numDomains),
+    agglom_(0)
+{}
+
+
 Foam::randomDecomp::randomDecomp
 (
     const dictionary& decompDict,
@@ -90,7 +97,7 @@ Foam::randomDecomp::randomDecomp
 )
 :
     decompositionMethod(decompDict, regionName),
-    agglom_(1)
+    agglom_(0)
 {
     const dictionary& coeffs = findCoeffsDict(typeName + "Coeffs", select);
 
@@ -104,16 +111,7 @@ Foam::randomDecomp::randomDecomp
 Foam::labelList Foam::randomDecomp::decompose
 (
     const pointField& points,
-    const scalarField& pointWeights
-) const
-{
-    return randomMap(points.size());
-}
-
-
-Foam::labelList Foam::randomDecomp::decompose
-(
-    const pointField& points
+    const scalarField&
 ) const
 {
     return randomMap(points.size());
@@ -127,7 +125,18 @@ Foam::labelList Foam::randomDecomp::decompose
     const scalarField&
 ) const
 {
-    return randomMap(mesh.nCells());  // or cc.size()
+    return randomMap(mesh.nCells());
+}
+
+
+Foam::labelList Foam::randomDecomp::decompose
+(
+    const CompactListList<label>& globalCellCells,
+    const pointField&,
+    const scalarField&
+) const
+{
+    return randomMap(globalCellCells.size());
 }
 
 
@@ -138,7 +147,7 @@ Foam::labelList Foam::randomDecomp::decompose
     const scalarField&
 ) const
 {
-    return randomMap(globalCellCells.size());  // or cc.size()
+    return randomMap(globalCellCells.size());
 }
 
 
diff --git a/src/parallel/decompose/decompositionMethods/randomDecomp/randomDecomp.H b/src/parallel/decompose/decompositionMethods/randomDecomp/randomDecomp.H
index 407c6e3182..5233de2621 100644
--- a/src/parallel/decompose/decompositionMethods/randomDecomp/randomDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/randomDecomp/randomDecomp.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,7 +36,7 @@ Description
     Coefficients:
     \table
         Property  | Description                             | Required | Default
-        agglom    | Number of cells to agglomerate          | no  | 1
+        agglom    | Number of cells to agglomerate          | no  | 0
     \endtable
 
 SourceFiles
@@ -62,7 +62,7 @@ class randomDecomp
 {
     // Private Data
 
-        //- Number of cells to agglomerate per random value (default: 1)
+        //- Number of cells to agglomerate per random value (default: 0)
         label agglom_;
 
 
@@ -71,14 +71,16 @@ class randomDecomp
         //- Random distribution on the 0-nCells interval
         labelList randomMap(const label nCells) const;
 
-        //- No copy construct
-        void operator=(const randomDecomp&) = delete;
+public:
 
-        //- No copy assignment
+    // Generated Methods
+
+        //- No copy construct
         randomDecomp(const randomDecomp&) = delete;
 
+        //- No copy assignment
+        void operator=(const randomDecomp&) = delete;
 
-public:
 
     //- Runtime type information
     TypeNameNoDebug("random");
@@ -86,6 +88,9 @@ public:
 
     // Constructors
 
+        //- Construct with number of domains (no coefficients or constraints)
+        explicit randomDecomp(const label numDomains);
+
         //- Construct for decomposition dictionary and optional region name
         explicit randomDecomp
         (
@@ -107,25 +112,28 @@ public:
             return true;
         }
 
-
-    // No topology (implemented by geometric decomposers)
-
         //- Return for every coordinate the wanted processor number.
         virtual labelList decompose
         (
             const pointField& points,
-            const scalarField& pointWeights
+            const scalarField& pointWeights_unused = scalarField::null()
         ) const;
 
-        //- Decompose with uniform weights on the points
-        virtual labelList decompose(const pointField& points) const;
-
         //- Return for every coordinate the wanted processor number.
         virtual labelList decompose
         (
             const polyMesh& mesh,
-            const pointField& cc,
-            const scalarField& cWeights
+            const pointField& cc_unused = pointField::null(),
+            const scalarField& cWeights_unused = scalarField::null()
+        ) const;
+
+        //- Return for every coordinate the wanted processor number.
+        //  Explicitly provided connectivity - does not use mesh_.
+        virtual labelList decompose
+        (
+            const CompactListList<label>& globalCellCells,
+            const pointField& cc_unused = pointField::null(),
+            const scalarField& cWeights_unused = scalarField::null()
         ) const;
 
         //- Return for every coordinate the wanted processor number.
@@ -133,8 +141,8 @@ public:
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
-            const pointField& cc,
-            const scalarField& cWeights
+            const pointField& cc_unused = pointField::null(),
+            const scalarField& cWeights_unused = scalarField::null()
         ) const;
 };
 
diff --git a/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C b/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
index 72edd426b7..8eeb87cc0c 100644
--- a/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -76,7 +76,7 @@ struct vectorLessOp
 } // End namespace Foam
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
 // assignToProcessorGroup : given nCells cells and nProcGroup processor
 // groups to share them, how do we share them out? Answer : each group
@@ -84,15 +84,18 @@ struct vectorLessOp
 // extra to make up the numbers. This should produce almost
 // perfect load balancing
 
-void Foam::simpleGeomDecomp::assignToProcessorGroup
+namespace Foam
+{
+
+static void assignToProcessorGroup
 (
     labelList& processorGroup,
     const label nProcGroup
 )
 {
-    label jump = processorGroup.size()/nProcGroup;
-    label jumpb = jump + 1;
-    label fstProcessorGroup = processorGroup.size() - jump*nProcGroup;
+    const label jump = processorGroup.size()/nProcGroup;
+    const label jumpb = jump + 1;
+    const label fstProcessorGroup = processorGroup.size() - jump*nProcGroup;
 
     label ind = 0;
     label j = 0;
@@ -118,7 +121,7 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
 }
 
 
-void Foam::simpleGeomDecomp::assignToProcessorGroup
+static void assignToProcessorGroup
 (
     labelList& processorGroup,
     const label nProcGroup,
@@ -139,6 +142,7 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
 
     const scalar jump = summedWeights/nProcGroup;
     const label nProcGroupM1 = nProcGroup - 1;
+
     scalar sumWeights = 0;
     label ind = 0;
     label j = 0;
@@ -160,6 +164,10 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
     }
 }
 
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::labelList Foam::simpleGeomDecomp::decomposeOneProc
 (
@@ -231,6 +239,11 @@ Foam::labelList Foam::simpleGeomDecomp::decomposeOneProc
     const scalarField& weights
 ) const
 {
+    if (weights.empty())
+    {
+        return decomposeOneProc(points);
+    }
+
     // construct a list for the final result
     labelList finalDecomp(points.size());
 
@@ -314,6 +327,12 @@ Foam::labelList Foam::simpleGeomDecomp::decomposeOneProc
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::simpleGeomDecomp::simpleGeomDecomp(const Vector<label>& divisions)
+:
+    geomDecomp(divisions)
+{}
+
+
 Foam::simpleGeomDecomp::simpleGeomDecomp
 (
     const dictionary& decompDict,
@@ -326,59 +345,51 @@ Foam::simpleGeomDecomp::simpleGeomDecomp
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::labelList Foam::simpleGeomDecomp::decompose
-(
-    const pointField& points
-) const
-{
-    if (!Pstream::parRun())
-    {
-        return decomposeOneProc(points);
-    }
-    else
-    {
-        const globalIndex globalNumbers(points.size());
-
-        pointField allPoints(globalNumbers.gather(points));
-
-        labelList allDecomp;
-        if (Pstream::master())
-        {
-            allDecomp = decomposeOneProc(allPoints);
-            allPoints.clear();  // Not needed anymore
-        }
-
-        return globalNumbers.scatter(allDecomp);
-    }
-}
-
-
 Foam::labelList Foam::simpleGeomDecomp::decompose
 (
     const pointField& points,
     const scalarField& weights
 ) const
 {
-    if (returnReduceOr(points.size() != weights.size()))
+    // Uniform weighting if weights are empty or poorly sized
+    const bool hasWeights = returnReduceAnd(points.size() == weights.size());
+
+    if (!UPstream::parRun())
     {
-        // Ignore zero-sized weights ... and poorly sized ones too
-        return decompose(points);
-    }
-    else if (!Pstream::parRun())
-    {
-        return decomposeOneProc(points, weights);
+        if (hasWeights)
+        {
+            return decomposeOneProc(points, weights);
+        }
+        else
+        {
+            return decomposeOneProc(points);
+        }
     }
     else
     {
+        // Parallel
         const globalIndex globalNumbers(points.size());
 
-        pointField allPoints(globalNumbers.gather(points));
-        scalarField allWeights(globalNumbers.gather(weights));
-
         labelList allDecomp;
-        if (Pstream::master())
+        pointField allPoints(globalNumbers.gather(points));
+        scalarField allWeights;
+
+        if (hasWeights)
         {
-            allDecomp = decomposeOneProc(allPoints, allWeights);
+            // Non-uniform weighting
+            allWeights = globalNumbers.gather(weights);
+        }
+
+        if (UPstream::master())
+        {
+            if (hasWeights)
+            {
+                allDecomp = decomposeOneProc(allPoints, allWeights);
+            }
+            else
+            {
+                allDecomp = decomposeOneProc(allPoints);
+            }
             allPoints.clear();  // Not needed anymore
             allWeights.clear(); // ...
         }
diff --git a/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H b/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
index 3e9bd76ab5..618fe2d599 100644
--- a/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,8 +44,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef simpleGeomDecomp_H
-#define simpleGeomDecomp_H
+#ifndef Foam_simpleGeomDecomp_H
+#define Foam_simpleGeomDecomp_H
 
 #include "geomDecomp.H"
 
@@ -62,21 +62,6 @@ class simpleGeomDecomp
 {
     // Private Member Functions
 
-        static void assignToProcessorGroup
-        (
-            labelList& processorGroup,
-            const label nProcGroup
-        );
-
-        static void assignToProcessorGroup
-        (
-            labelList& processorGroup,
-            const label nProcGroup,
-            const labelList& indices,
-            const scalarField& weights,
-            const scalar summedWeights
-        );
-
         labelList decomposeOneProc(const pointField& points) const;
 
         labelList decomposeOneProc
@@ -88,11 +73,14 @@ class simpleGeomDecomp
 
 public:
 
-    //- No copy construct
-    void operator=(const simpleGeomDecomp&) = delete;
+    // Generated Methods
+
+        //- No copy construct
+        simpleGeomDecomp(const simpleGeomDecomp&) = delete;
+
+        //- No copy assignment
+        void operator=(const simpleGeomDecomp&) = delete;
 
-    //- No copy assignment
-    simpleGeomDecomp(const simpleGeomDecomp&) = delete;
 
     //- Runtime type information
     TypeName("simple");
@@ -100,6 +88,10 @@ public:
 
     // Constructors
 
+        //- Construct with number of x/y/z division
+        //- (no coefficients or constraints)
+        explicit simpleGeomDecomp(const Vector<label>& divisions);
+
         //- Construct given decomposition dictionary and optional region name
         explicit simpleGeomDecomp
         (
@@ -114,51 +106,49 @@ public:
 
     // Member Functions
 
-        //- Simple sends all points to the master for decomposition.
+        //- Sends all points to the master for decomposition.
         virtual bool parallelAware() const
         {
             return true;
         }
 
-        //- Decompose with uniform weights.
-        virtual labelList decompose(const pointField& points) const;
-
         //- Return for every coordinate the wanted processor number.
+        //- using uniform or specified point weights.
         virtual labelList decompose
         (
             const pointField& points,
-            const scalarField& weights
+            const scalarField& weights = scalarField::null()
         ) const;
 
-        //- Decompose with uniform weights.
+        //- Decompose with uniform or specified point weights.
         virtual labelList decompose
         (
-            const polyMesh& mesh,
-            const pointField& points
-        ) const
-        {
-            checkDecompositionDirections(mesh.geometricD());
-            return decompose(points);
-        }
-
-        //- Return for every coordinate the wanted processor number.
-        virtual labelList decompose
-        (
-            const polyMesh& mesh,
+            const polyMesh& mesh,  //!< To check mesh dimensions
             const pointField& points,
-            const scalarField& weights
+            const scalarField& pointWeights = scalarField::null()
         ) const
         {
             checkDecompositionDirections(mesh.geometricD());
-            return decompose(points, weights);
+            return decompose(points, pointWeights);
         }
 
         //- Explicitly provided connectivity
         virtual labelList decompose
         (
-            const labelListList& globalCellCells,  // unused
+            const CompactListList<label>& globalCellCells,  //!< unused
             const pointField& cc,
-            const scalarField& cWeights
+            const scalarField& cWeights = scalarField::null()
+        ) const
+        {
+            return decompose(cc, cWeights);
+        }
+
+        //- Explicitly provided connectivity
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,  //!< unused
+            const pointField& cc,
+            const scalarField& cWeights = scalarField::null()
         ) const
         {
             return decompose(cc, cWeights);
diff --git a/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C b/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C
index 5c6d7e3825..e33541667a 100644
--- a/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -102,7 +102,11 @@ Foam::labelList Foam::structuredDecomp::decompose
     );
     const fvMesh& subMesh = subsetter.subMesh();
     pointField subCc(cc, subsetter.cellMap());
-    scalarField subWeights(cWeights, subsetter.cellMap());
+    scalarField subWeights;
+    if (cWeights.size() == cc.size())
+    {
+        subWeights = scalarField(cWeights, subsetter.cellMap());
+    }
 
     // Decompose the layer of cells
     labelList subDecomp(method_().decompose(subMesh, subCc, subWeights));
diff --git a/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.H b/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.H
index 29b9eafcc2..c9aeee8e37 100644
--- a/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef structuredDecomp_H
-#define structuredDecomp_H
+#ifndef Foam_structuredDecomp_H
+#define Foam_structuredDecomp_H
 
 #include "decompositionMethod.H"
 #include "wordRes.H"
@@ -61,7 +61,9 @@ class structuredDecomp
         autoPtr<decompositionMethod> method_;
 
 
-    // Private Member Functions
+public:
+
+    // Generated Methods
 
         //- No copy construct
         structuredDecomp(const structuredDecomp&) = delete;
@@ -70,8 +72,6 @@ class structuredDecomp
         void operator=(const structuredDecomp&) = delete;
 
 
-public:
-
     //- Runtime type information
     TypeName("structured");
 
@@ -102,16 +102,29 @@ public:
         (
             const polyMesh& mesh,
             const pointField& points,
-            const scalarField& pointWeights
+            const scalarField& pointWeights = scalarField::null()
         ) const;
 
+        //- Return for every coordinate the wanted processor number.
+        //  Explicitly provided connectivity - does not use mesh_.
+        virtual labelList decompose
+        (
+            const CompactListList<label>& globalCellCells,
+            const pointField& cc,
+            const scalarField& cWeights = scalarField::null()
+        ) const
+        {
+            NotImplemented;
+            return labelList();
+        }
+
         //- Return for every coordinate the wanted processor number.
         //  Explicitly provided connectivity - does not use mesh_.
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
             const pointField& cc,
-            const scalarField& cWeights
+            const scalarField& cWeights = scalarField::null()
         ) const
         {
             NotImplemented;
diff --git a/src/parallel/decompose/kahipDecomp/kahipDecomp.C b/src/parallel/decompose/kahipDecomp/kahipDecomp.C
index 8d29fdd6ac..b889efa020 100644
--- a/src/parallel/decompose/kahipDecomp/kahipDecomp.C
+++ b/src/parallel/decompose/kahipDecomp/kahipDecomp.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -108,7 +108,7 @@ Foam::label Foam::kahipDecomp::decomposeSerial
     ConstPrecisionAdaptor<int, label, List> xadj_param(xadj);
 
     // Output: cell -> processor addressing
-    decomp.resize(numCells);
+    decomp.resize_nocopy(numCells);
     decomp = 0;
     PrecisionAdaptor<int, label, List> decomp_param(decomp, false);
 
@@ -149,8 +149,8 @@ Foam::label Foam::kahipDecomp::decomposeSerial
     else if (hasWeights && (cWeights.size() != numCells))
     {
         FatalErrorInFunction
-            << "Number of cell weights " << cWeights.size()
-            << " does not equal number of cells " << numCells
+            << "Number of weights (" << cWeights.size()
+            << ") != number of cells (" << numCells << ")"
             << exit(FatalError);
     }
 
@@ -160,7 +160,7 @@ Foam::label Foam::kahipDecomp::decomposeSerial
     if (hasWeights)
     {
         // Convert to integers.
-        cellWeights.resize(cWeights.size());
+        cellWeights.resize_nocopy(cWeights.size());
         forAll(cellWeights, i)
         {
             cellWeights[i] = static_cast<int>
@@ -299,6 +299,12 @@ Foam::label Foam::kahipDecomp::decomposeSerial
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::kahipDecomp::kahipDecomp(const label numDomains)
+:
+    metisLikeDecomp(numDomains)
+{}
+
+
 Foam::kahipDecomp::kahipDecomp
 (
     const dictionary& decompDict,
diff --git a/src/parallel/decompose/kahipDecomp/kahipDecomp.H b/src/parallel/decompose/kahipDecomp/kahipDecomp.H
index 0b83b93616..dece23b282 100644
--- a/src/parallel/decompose/kahipDecomp/kahipDecomp.H
+++ b/src/parallel/decompose/kahipDecomp/kahipDecomp.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,8 +59,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef kahipDecomp_H
-#define kahipDecomp_H
+#ifndef Foam_kahipDecomp_H
+#define Foam_kahipDecomp_H
 
 #include "metisLikeDecomp.H"
 #include "Enum.H"
@@ -90,13 +90,6 @@ protected:
         ) const;
 
 
-        //- No copy construct
-        kahipDecomp(const kahipDecomp&) = delete;
-
-        //- No copy assignment
-        void operator=(const kahipDecomp&) = delete;
-
-
 public:
 
     //- The predefined KaHIP configuration types
@@ -115,12 +108,24 @@ public:
     static const Enum<configs> configNames;
 
 
+    // Generated Methods
+
+        //- No copy construct
+        kahipDecomp(const kahipDecomp&) = delete;
+
+        //- No copy assignment
+        void operator=(const kahipDecomp&) = delete;
+
+
     //- Runtime type information
     TypeName("kahip");
 
 
     // Constructors
 
+        //- Construct with number of domains (no coefficients or constraints)
+        explicit kahipDecomp(const label numDomains);
+
         //- Construct given decomposition dictionary and optional region name
         explicit kahipDecomp
         (
@@ -135,6 +140,7 @@ public:
 
     // Member Functions
 
+        //- Method is parallel aware
         virtual bool parallelAware() const
         {
             return true;
diff --git a/src/parallel/decompose/metisDecomp/metisDecomp.C b/src/parallel/decompose/metisDecomp/metisDecomp.C
index efa706b58b..d62307cb48 100644
--- a/src/parallel/decompose/metisDecomp/metisDecomp.C
+++ b/src/parallel/decompose/metisDecomp/metisDecomp.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -101,29 +101,31 @@ Foam::label Foam::metisDecomp::decomposeSerial
     // Face weights (so on the edges of the dual)
     List<idx_t> faceWeights;
 
-    // Check for externally provided cellweights and if so initialise weights
+    bool hasWeights = !cWeights.empty();
+
     // Note: min, not gMin since routine runs on master only.
-    const scalar minWeights = min(cWeights);
+    const scalar minWeights = hasWeights ? min(cWeights) : scalar(1);
 
-    if (!cWeights.empty())
+    if (minWeights <= 0)
     {
-        if (minWeights <= 0)
-        {
-            WarningInFunction
-                << "Illegal minimum weight " << minWeights
-                << endl;
-        }
-
-        if (cWeights.size() != numCells)
-        {
-            FatalErrorInFunction
-                << "Number of cell weights " << cWeights.size()
-                << " does not equal number of cells " << numCells
-                << exit(FatalError);
-        }
+        hasWeights = false;
+        WarningInFunction
+            << "Illegal minimum weight " << minWeights
+            << " ... ignoring"
+            << endl;
+    }
+    else if (hasWeights && (cWeights.size() != numCells))
+    {
+        FatalErrorInFunction
+            << "Number of weights (" << cWeights.size()
+            << ") != number of cells (" << numCells << ")"
+            << exit(FatalError);
+    }
 
+    if (hasWeights)
+    {
         // Convert to integers.
-        cellWeights.setSize(cWeights.size());
+        cellWeights.resize_nocopy(cWeights.size());
         forAll(cellWeights, i)
         {
             cellWeights[i] = idx_t(cWeights[i]/minWeights);
@@ -261,6 +263,12 @@ Foam::label Foam::metisDecomp::decomposeSerial
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::metisDecomp::metisDecomp(const label numDomains)
+:
+    metisLikeDecomp(numDomains)
+{}
+
+
 Foam::metisDecomp::metisDecomp
 (
     const dictionary& decompDict,
diff --git a/src/parallel/decompose/metisDecomp/metisDecomp.H b/src/parallel/decompose/metisDecomp/metisDecomp.H
index 0578000263..4a2feb0a46 100644
--- a/src/parallel/decompose/metisDecomp/metisDecomp.H
+++ b/src/parallel/decompose/metisDecomp/metisDecomp.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -60,8 +60,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef metisDecomp_H
-#define metisDecomp_H
+#ifndef Foam_metisDecomp_H
+#define Foam_metisDecomp_H
 
 #include "metisLikeDecomp.H"
 
@@ -90,6 +90,10 @@ protected:
         ) const;
 
 
+public:
+
+    // Generated Methods
+
         //- No copy construct
         metisDecomp(const metisDecomp&) = delete;
 
@@ -97,14 +101,15 @@ protected:
         void operator=(const metisDecomp&) = delete;
 
 
-public:
-
     //- Runtime type information
     TypeName("metis");
 
 
     // Constructors
 
+        //- Construct with number of domains (no coefficients or constraints)
+        explicit metisDecomp(const label numDomains);
+
         //- Construct given decomposition dictionary and optional region name
         explicit metisDecomp
         (
@@ -119,6 +124,7 @@ public:
 
     // Member Functions
 
+        //- Knows about coupled boundaries
         virtual bool parallelAware() const
         {
             return true;
diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
index cbe191fa8f..7a2ec4948d 100644
--- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
+++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
@@ -29,6 +29,7 @@ License
 #include "ptscotchDecomp.H"
 #include "addToRunTimeSelectionTable.H"
 #include "floatScalar.H"
+#include "globalMeshData.H"
 #include "Time.H"
 #include "PrecisionAdaptor.H"
 #include "OFstream.H"
@@ -123,7 +124,7 @@ Foam::label Foam::ptscotchDecomp::decompose
     ConstPrecisionAdaptor<SCOTCH_Num, label, List> xadj_param(xadj);
 
     // Output: cell -> processor addressing
-    decomp.resize(numCells);
+    decomp.resize_nocopy(numCells);
     decomp = 0;
     PrecisionAdaptor<SCOTCH_Num, label, List> decomp_param(decomp, false);
 
@@ -497,6 +498,13 @@ Foam::label Foam::ptscotchDecomp::decompose
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::ptscotchDecomp::ptscotchDecomp(const label numDomains)
+:
+    decompositionMethod(numDomains),
+    coeffsDict_(dictionary::null)
+{}
+
+
 Foam::ptscotchDecomp::ptscotchDecomp
 (
     const dictionary& decompDict,
@@ -520,23 +528,20 @@ Foam::labelList Foam::ptscotchDecomp::decompose
     // Where to write graph
     graphPath_ = getGraphPathBase(mesh);
 
-    if (points.size() != mesh.nCells())
+    if (!points.empty() && (points.size() != mesh.nCells()))
     {
         FatalErrorInFunction
-            << "Can only use this decomposition method for entire mesh" << nl
-            << "and supply one coordinate (cellCentre) for every cell." << nl
-            << "The number of coordinates " << points.size() << nl
-            << "The number of cells in the mesh " << mesh.nCells() << nl
+            << "Number of cell centres (" << points.size()
+            << ") != number of cells (" << mesh.nCells() << ")"
             << exit(FatalError);
     }
 
-
     // 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;
-    calcCellCells
+    globalMeshData::calcCellCells
     (
         mesh,
         identity(mesh.nCells()),
@@ -564,7 +569,7 @@ Foam::labelList Foam::ptscotchDecomp::decompose
     const polyMesh& mesh,
     const labelList& agglom,
     const pointField& agglomPoints,
-    const scalarField& pointWeights
+    const scalarField& agglomWeights
 ) const
 {
     // Where to write graph
@@ -573,8 +578,8 @@ Foam::labelList Foam::ptscotchDecomp::decompose
     if (agglom.size() != mesh.nCells())
     {
         FatalErrorInFunction
-            << "Size of cell-to-coarse map " << agglom.size()
-            << " differs from number of cells in mesh " << mesh.nCells()
+            << "Agglomeration size (" << agglom.size()
+            << ") != number of cells (" << mesh.nCells() << ")"
             << exit(FatalError);
     }
 
@@ -583,7 +588,7 @@ Foam::labelList Foam::ptscotchDecomp::decompose
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
     CompactListList<label> cellCells;
-    calcCellCells
+    globalMeshData::calcCellCells
     (
         mesh,
         agglom,
@@ -598,19 +603,49 @@ Foam::labelList Foam::ptscotchDecomp::decompose
     (
         cellCells.values(),
         cellCells.offsets(),
-        pointWeights,
+        agglomWeights,
         decomp
     );
 
-    // Rework back into decomposition for original mesh
-    labelList fineDistribution(agglom.size());
+    // From coarse back to fine for original mesh
+    return labelList(decomp, agglom);
+}
 
-    forAll(fineDistribution, i)
+
+Foam::labelList Foam::ptscotchDecomp::decompose
+(
+    const CompactListList<label>& globalCellCells,
+    const pointField& cellCentres,
+    const scalarField& cWeights
+) const
+{
+    // Where to write graph
+    graphPath_ = "ptscotch";
+
+    if (!cellCentres.empty() && (cellCentres.size() != globalCellCells.size()))
     {
-        fineDistribution[i] = decomp[agglom[i]];
+        FatalErrorInFunction
+            << "Number of cell centres (" << cellCentres.size()
+            << ") != number of cells (" << globalCellCells.size() << ")"
+            << exit(FatalError);
     }
 
-    return fineDistribution;
+    // CompactListList is already
+    // Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+
+    // Decompose using default weights
+    labelList decomp;
+    decompose
+    (
+        globalCellCells.values(),
+        globalCellCells.offsets(),
+        cWeights,
+        decomp
+    );
+
+    return decomp;
 }
 
 
@@ -624,22 +659,21 @@ Foam::labelList Foam::ptscotchDecomp::decompose
     // Where to write graph
     graphPath_ = "ptscotch";
 
-    if (cellCentres.size() != globalCellCells.size())
+    if (!cellCentres.empty() && (cellCentres.size() != globalCellCells.size()))
     {
         FatalErrorInFunction
-            << "Inconsistent number of cells (" << globalCellCells.size()
-            << ") and number of cell centres (" << cellCentres.size()
-            << ")." << exit(FatalError);
+            << "Number of cell centres (" << cellCentres.size()
+            << ") != number of cells (" << globalCellCells.size() << ")"
+            << exit(FatalError);
     }
 
-
     // Make Metis CSR (Compressed Storage Format) storage
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
 
     auto cellCells(CompactListList<label>::pack(globalCellCells));
 
-    // Decompose using weights
+    // Decompose using default weights
     labelList decomp;
     decompose
     (
diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H
index 6467fd86e0..956b2aee95 100644
--- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H
+++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,8 +54,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef ptscotchDecomp_H
-#define ptscotchDecomp_H
+#ifndef Foam_ptscotchDecomp_H
+#define Foam_ptscotchDecomp_H
 
 #include "decompositionMethod.H"
 
@@ -90,6 +90,10 @@ class ptscotchDecomp
             labelList& decomp
         ) const;
 
+public:
+
+    // Generated Methods
+
         //- No copy construct
         ptscotchDecomp(const ptscotchDecomp&) = delete;
 
@@ -97,14 +101,15 @@ class ptscotchDecomp
         void operator=(const ptscotchDecomp&) = delete;
 
 
-public:
-
     //- Runtime type information
     TypeName("ptscotch");
 
 
     // Constructors
 
+        //- Construct with number of domains (no coefficients or constraints)
+        explicit ptscotchDecomp(const label numDomains);
+
         //- Construct given decomposition dictionary and optional region name
         explicit ptscotchDecomp
         (
@@ -119,6 +124,7 @@ public:
 
     // Member Functions
 
+        //- Knows about coupled boundaries
         virtual bool parallelAware() const
         {
             return true;
@@ -133,8 +139,8 @@ public:
         virtual labelList decompose
         (
             const polyMesh& mesh,
-            const pointField& points,
-            const scalarField& pointWeights
+            const pointField& points = pointField::null(),
+            const scalarField& pointWeights = scalarField::null()
         ) const;
 
         //- Return for every coordinate the wanted processor number.
@@ -146,8 +152,16 @@ public:
         (
             const polyMesh& mesh,
             const labelList& agglom,
-            const pointField& regionPoints,
-            const scalarField& regionWeights
+            const pointField& agglomPoints,
+            const scalarField& agglomWeights = scalarField::null()
+        ) const;
+
+        //- Return for every coordinate the wanted processor number.
+        virtual labelList decompose
+        (
+            const CompactListList<label>& globalCellCells,
+            const pointField& cc = pointField::null(),
+            const scalarField& cWeights = scalarField::null()
         ) const;
 
         //- Return for every coordinate the wanted processor number.
@@ -161,8 +175,8 @@ public:
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
-            const pointField& cc,
-            const scalarField& cWeights
+            const pointField& cc = pointField::null(),
+            const scalarField& cWeights = scalarField::null()
         ) const;
 };
 
diff --git a/src/parallel/decompose/scotchDecomp/scotchDecomp.C b/src/parallel/decompose/scotchDecomp/scotchDecomp.C
index 1768b7169e..e12252f806 100644
--- a/src/parallel/decompose/scotchDecomp/scotchDecomp.C
+++ b/src/parallel/decompose/scotchDecomp/scotchDecomp.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2021,2023 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -120,7 +120,7 @@ Foam::label Foam::scotchDecomp::decomposeSerial
     ConstPrecisionAdaptor<SCOTCH_Num, label, List> xadj_param(xadj);
 
     // Output: cell -> processor addressing
-    decomp.resize(numCells);
+    decomp.resize_nocopy(numCells);
     decomp = 0;
     PrecisionAdaptor<SCOTCH_Num, label, List> decomp_param(decomp, false);
 
@@ -512,6 +512,12 @@ Foam::label Foam::scotchDecomp::decomposeSerial
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::scotchDecomp::scotchDecomp(const label numDomains)
+:
+    metisLikeDecomp(numDomains)
+{}
+
+
 Foam::scotchDecomp::scotchDecomp
 (
     const dictionary& decompDict,
@@ -548,7 +554,7 @@ Foam::labelList Foam::scotchDecomp::decompose
     const polyMesh& mesh,
     const labelList& agglom,
     const pointField& agglomPoints,
-    const scalarField& pointWeights
+    const scalarField& agglomWeights
 ) const
 {
     // Where to write graph
@@ -559,7 +565,26 @@ Foam::labelList Foam::scotchDecomp::decompose
         mesh,
         agglom,
         agglomPoints,
-        pointWeights
+        agglomWeights
+    );
+}
+
+
+Foam::labelList Foam::scotchDecomp::decompose
+(
+    const CompactListList<label>& globalCellCells,
+    const pointField& cellCentres,
+    const scalarField& cWeights
+) const
+{
+    // Where to write graph
+    graphPath_ = "scotch.grf";
+
+    return metisLikeDecomp::decompose
+    (
+        globalCellCells,
+        cellCentres,
+        cWeights
     );
 }
 
diff --git a/src/parallel/decompose/scotchDecomp/scotchDecomp.H b/src/parallel/decompose/scotchDecomp/scotchDecomp.H
index 6ce9588c97..b8b5dd4cee 100644
--- a/src/parallel/decompose/scotchDecomp/scotchDecomp.H
+++ b/src/parallel/decompose/scotchDecomp/scotchDecomp.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2021,2023 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -249,8 +249,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef scotchDecomp_H
-#define scotchDecomp_H
+#ifndef Foam_scotchDecomp_H
+#define Foam_scotchDecomp_H
 
 #include "metisLikeDecomp.H"
 
@@ -284,6 +284,9 @@ protected:
             labelList& decomp
         ) const;
 
+public:
+
+    // Generated Methods
 
         //- No copy construct
         scotchDecomp(const scotchDecomp&) = delete;
@@ -292,14 +295,15 @@ protected:
         void operator=(const scotchDecomp&) = delete;
 
 
-public:
-
     //- Runtime type information
     TypeName("scotch");
 
 
     // Constructors
 
+        //- Construct with number of domains (no coefficients or constraints)
+        explicit scotchDecomp(const label numDomains);
+
         //- Construct given decomposition dictionary and optional region name
         explicit scotchDecomp
         (
@@ -327,25 +331,33 @@ public:
         virtual labelList decompose
         (
             const polyMesh& mesh,
-            const pointField& points,
-            const scalarField& pointWeights
+            const pointField& points = pointField::null(),
+            const scalarField& pointWeights = scalarField::null()
         ) const;
 
         //- Return for every coordinate the wanted processor number.
         virtual labelList decompose
         (
             const polyMesh& mesh,
-            const labelList& agglom,
-            const pointField& regionPoints,
-            const scalarField& regionWeights
+            const labelList& agglom,  //!< agglomeration: fine-to-coarse
+            const pointField& agglomPoints,
+            const scalarField& agglomWeights = scalarField::null()
         ) const;
 
-        //- Return for every coordinate the wanted processor number.
+        //- Return the wanted processor number for every cell.
+        virtual labelList decompose
+        (
+            const CompactListList<label>& globalCellCells,
+            const pointField& cc = pointField::null(),
+            const scalarField& cWeights = scalarField::null()
+        ) const;
+
+        //- Return the wanted processor number for every cell.
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
-            const pointField& cc,
-            const scalarField& cWeights
+            const pointField& cc = pointField::null(),
+            const scalarField& cWeights = scalarField::null()
         ) const;
 };
 
diff --git a/src/renumber/SloanRenumber/Make/options b/src/renumber/SloanRenumber/Make/options
index 5af9c713ee..c47a6d48da 100644
--- a/src/renumber/SloanRenumber/Make/options
+++ b/src/renumber/SloanRenumber/Make/options
@@ -2,11 +2,9 @@ EXE_INC = \
     $(c++LESSWARN) \
     -I$(BOOST_INC_DIR) \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/renumber/renumberMethods/lnInclude
 
 LIB_LIBS = \
     -L$(BOOST_LIB_DIR) -lboost_system \
     -lmeshTools \
-    -ldecompositionMethods \
     -lrenumberMethods
diff --git a/src/renumber/SloanRenumber/SloanRenumber.C b/src/renumber/SloanRenumber/SloanRenumber.C
index 04b8f265ca..f3dea4f762 100644
--- a/src/renumber/SloanRenumber/SloanRenumber.C
+++ b/src/renumber/SloanRenumber/SloanRenumber.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2017 OpenFOAM Foundation
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,7 +30,7 @@ License
 
 #include "SloanRenumber.H"
 #include "addToRunTimeSelectionTable.H"
-#include "decompositionMethod.H"
+#include "globalMeshData.H"
 #include "processorPolyPatch.H"
 #include "syncTools.H"
 
diff --git a/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C b/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
index fe8bcc2936..b2662f90f6 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-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,7 +29,7 @@ License
 #include "CuthillMcKeeRenumber.H"
 #include "addToRunTimeSelectionTable.H"
 #include "bandCompression.H"
-#include "decompositionMethod.H"
+#include "globalMeshData.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -68,7 +68,7 @@ Foam::labelList Foam::CuthillMcKeeRenumber::renumber
 ) const
 {
     CompactListList<label> cellCells;
-    decompositionMethod::calcCellCells
+    globalMeshData::calcCellCells
     (
         mesh,
         identity(mesh.nCells()),
diff --git a/src/renumber/renumberMethods/Make/options b/src/renumber/renumberMethods/Make/options
index 0cd60f1b7a..a3ae8da833 100644
--- a/src/renumber/renumberMethods/Make/options
+++ b/src/renumber/renumberMethods/Make/options
@@ -1,9 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
-    -lmeshTools \
-    -ldecompositionMethods
+    -lmeshTools
diff --git a/src/renumber/renumberMethods/renumberMethod/renumberMethod.C b/src/renumber/renumberMethods/renumberMethod/renumberMethod.C
index 3c2b1acf34..2ea1cf0c7e 100644
--- a/src/renumber/renumberMethods/renumberMethod/renumberMethod.C
+++ b/src/renumber/renumberMethods/renumberMethod/renumberMethod.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +27,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "renumberMethod.H"
-#include "decompositionMethod.H"
+#include "globalMeshData.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -74,7 +74,7 @@ Foam::labelList Foam::renumberMethod::renumber
 ) const
 {
     CompactListList<label> cellCells;
-    decompositionMethod::calcCellCells
+    globalMeshData::calcCellCells
     (
         mesh,
         identity(mesh.nCells()),
@@ -117,7 +117,7 @@ Foam::labelList Foam::renumberMethod::renumber
 ) const
 {
     CompactListList<label> coarseCellCells;
-    decompositionMethod::calcCellCells
+    globalMeshData::calcCellCells
     (
         mesh,
         fineToCoarse,
@@ -132,15 +132,8 @@ Foam::labelList Foam::renumberMethod::renumber
         renumber(coarseCellCells, coarsePoints)
     );
 
-    // Rework back into renumbering for original mesh_
-    labelList fineDistribution(fineToCoarse.size());
-
-    forAll(fineDistribution, i)
-    {
-        fineDistribution[i] = coarseDistribution[fineToCoarse[i]];
-    }
-
-    return fineDistribution;
+    // From coarse back to fine for original mesh
+    return labelList(coarseDistribution, fineToCoarse);
 }
 
 
diff --git a/src/renumber/renumberMethods/springRenumber/springRenumber.C b/src/renumber/renumberMethods/springRenumber/springRenumber.C
index ce6370649f..f7b79dfc93 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-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,8 +27,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "springRenumber.H"
+#include "globalMeshData.H"
 #include "addToRunTimeSelectionTable.H"
-#include "decompositionMethod.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -164,7 +164,7 @@ Foam::labelList Foam::springRenumber::renumber
 ) const
 {
     CompactListList<label> cellCells;
-    decompositionMethod::calcCellCells
+    globalMeshData::calcCellCells
     (
         mesh,
         identity(mesh.nCells()),