diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 7ced14726d..c48abfbc39 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -334,8 +334,7 @@ int main(int argc, char *argv[])
     (
         decompositionMethod::New
         (
-            decomposeDict,
-            mesh
+            decomposeDict
         )
     );
     decompositionMethod& decomposer = decomposerPtr();
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
index 19a80bac5f..3394b928b2 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
@@ -489,11 +489,17 @@ int main(int argc, char *argv[])
         );
         autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
         (
-            decomposeDict,
-            mesh
+            decomposeDict
         );
 
-        labelList cellToRegion(decomposePtr().decompose(mesh.cellCentres()));
+        labelList cellToRegion
+        (
+            decomposePtr().decompose
+            (
+                mesh,
+                mesh.cellCentres()
+            )
+        );
 
         // For debugging: write out region
         {
diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index e9c20f1afb..7e3cd14588 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -381,7 +381,8 @@ void subsetVolFields
     const fvMesh& mesh,
     const fvMesh& subMesh,
     const labelList& cellMap,
-    const labelList& faceMap
+    const labelList& faceMap,
+    const labelHashSet& addedPatches
 )
 {
     const labelList patchMap(identity(mesh.boundaryMesh().size()));
@@ -412,14 +413,7 @@ void subsetVolFields
         //       get initialised.
         forAll(tSubFld().boundaryField(), patchI)
         {
-            const fvPatchField<typename GeoField::value_type>& pfld =
-                tSubFld().boundaryField()[patchI];
-
-            if
-            (
-                isA<calculatedFvPatchField<typename GeoField::value_type> >
-                (pfld)
-            )
+            if (addedPatches.found(patchI))
             {
                 tSubFld().boundaryField()[patchI] ==
                     pTraits<typename GeoField::value_type>::zero;
@@ -440,7 +434,8 @@ void subsetSurfaceFields
 (
     const fvMesh& mesh,
     const fvMesh& subMesh,
-    const labelList& faceMap
+    const labelList& faceMap,
+    const labelHashSet& addedPatches
 )
 {
     const labelList patchMap(identity(mesh.boundaryMesh().size()));
@@ -470,14 +465,7 @@ void subsetSurfaceFields
         //       get initialised.
         forAll(tSubFld().boundaryField(), patchI)
         {
-            const fvsPatchField<typename GeoField::value_type>& pfld =
-                tSubFld().boundaryField()[patchI];
-
-            if
-            (
-                isA<calculatedFvsPatchField<typename GeoField::value_type> >
-                (pfld)
-            )
+            if (addedPatches.found(patchI))
             {
                 tSubFld().boundaryField()[patchI] ==
                     pTraits<typename GeoField::value_type>::zero;
@@ -852,6 +840,15 @@ void createAndWriteRegion
         newMesh
     );
 
+
+    // Make map of all added patches
+    labelHashSet addedPatches(2*interfaceToPatch.size());
+    forAllConstIter(EdgeMap<label>, interfaceToPatch, iter)
+    {
+        addedPatches.insert(iter());
+        addedPatches.insert(iter()+1);
+    }
+
     Info<< "Mapping fields" << endl;
 
     // Map existing fields
@@ -863,66 +860,76 @@ void createAndWriteRegion
         mesh,
         newMesh(),
         map().cellMap(),
-        map().faceMap()
+        map().faceMap(),
+        addedPatches
     );
     subsetVolFields<volVectorField>
     (
         mesh,
         newMesh(),
         map().cellMap(),
-        map().faceMap()
+        map().faceMap(),
+        addedPatches
     );
     subsetVolFields<volSphericalTensorField>
     (
         mesh,
         newMesh(),
         map().cellMap(),
-        map().faceMap()
+        map().faceMap(),
+        addedPatches
     );
     subsetVolFields<volSymmTensorField>
     (
         mesh,
         newMesh(),
         map().cellMap(),
-        map().faceMap()
+        map().faceMap(),
+        addedPatches
     );
     subsetVolFields<volTensorField>
     (
         mesh,
         newMesh(),
         map().cellMap(),
-        map().faceMap()
+        map().faceMap(),
+        addedPatches
     );
 
     subsetSurfaceFields<surfaceScalarField>
     (
         mesh,
         newMesh(),
-        map().faceMap()
+        map().faceMap(),
+        addedPatches
     );
     subsetSurfaceFields<surfaceVectorField>
     (
         mesh,
         newMesh(),
-        map().faceMap()
+        map().faceMap(),
+        addedPatches
     );
     subsetSurfaceFields<surfaceSphericalTensorField>
     (
         mesh,
         newMesh(),
-        map().faceMap()
+        map().faceMap(),
+        addedPatches
     );
     subsetSurfaceFields<surfaceSymmTensorField>
     (
         mesh,
         newMesh(),
-        map().faceMap()
+        map().faceMap(),
+        addedPatches
     );
     subsetSurfaceFields<surfaceTensorField>
     (
         mesh,
         newMesh(),
-        map().faceMap()
+        map().faceMap(),
+        addedPatches
     );
 
 
diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposeParDict b/applications/utilities/parallelProcessing/decomposePar/decomposeParDict
index a183ebd0cb..0f09fca353 100644
--- a/applications/utilities/parallelProcessing/decomposePar/decomposeParDict
+++ b/applications/utilities/parallelProcessing/decomposePar/decomposeParDict
@@ -19,6 +19,7 @@ FoamFile
 
 numberOfSubdomains  4;
 
+
 //- Keep owner and neighbour on same processor for faces in zones:
 // preserveFaceZones (heater solid1 solid3);
 
@@ -28,15 +29,40 @@ numberOfSubdomains  4;
 //preservePatches (cyclic_half0 cyclic_half1);
 
 
-method          scotch;
+// method          scotch;
 // method          hierarchical;
 // method          simple;
 // method          metis;
 // method          manual;
+method          multiLevel;
+
+multiLevelCoeffs
+{
+    level0
+    {
+        numberOfSubdomains  2;
+        method simple;
+        simpleCoeffs
+        {
+            n           (2 1 1);
+            delta       0.001;
+        }
+    }
+    level1
+    {
+        numberOfSubdomains  2;
+        method scotch;
+    }
+}
+
+
+// Desired output
+
+
 
 simpleCoeffs
 {
-    n           (2 2 1);
+    n           (2 1 1);
     delta       0.001;
 }
 
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C
index 6adf063962..877395e38d 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C
@@ -109,13 +109,12 @@ void Foam::domainDecomposition::distributeCells()
 
     autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
     (
-        decompositionDict_,
-        *this
+        decompositionDict_
     );
 
     if (sameProcFaces.empty())
     {
-        cellToProc_ = decomposePtr().decompose(cellCentres());
+        cellToProc_ = decomposePtr().decompose(*this, cellCentres());
     }
     else
     {
@@ -174,7 +173,12 @@ void Foam::domainDecomposition::distributeCells()
 
         // Do decomposition on agglomeration
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        cellToProc_ = decomposePtr().decompose(globalRegion, regionCentres);
+        cellToProc_ = decomposePtr().decompose
+        (
+            *this,
+            globalRegion,
+            regionCentres
+        );
     }
 
     Info<< "\nFinished decomposition in "
diff --git a/applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C b/applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C
index 5df7310f8e..173d6314df 100644
--- a/applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C
+++ b/applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C
@@ -594,8 +594,7 @@ int main(int argc, char *argv[])
         (
             decompositionMethod::New
             (
-                decompositionDict,
-                mesh
+                decompositionDict
             )
         );
 
@@ -612,7 +611,7 @@ int main(int argc, char *argv[])
                 << endl;
         }
 
-        finalDecomp = decomposer().decompose(mesh.cellCentres());
+        finalDecomp = decomposer().decompose(mesh, mesh.cellCentres());
     }
 
     // Dump decomposition to volScalarField
diff --git a/src/dummyThirdParty/metisDecomp/dummyMetisDecomp.C b/src/dummyThirdParty/metisDecomp/dummyMetisDecomp.C
index e7b61fd28d..4be4cd995b 100644
--- a/src/dummyThirdParty/metisDecomp/dummyMetisDecomp.C
+++ b/src/dummyThirdParty/metisDecomp/dummyMetisDecomp.C
@@ -46,7 +46,7 @@ namespace Foam
     (
         decompositionMethod,
         metisDecomp,
-        dictionaryMesh
+        dictionary
     );
 }
 
@@ -80,12 +80,10 @@ Foam::label Foam::metisDecomp::decompose
 
 Foam::metisDecomp::metisDecomp
 (
-    const dictionary& decompositionDict,
-    const polyMesh& mesh
+    const dictionary& decompositionDict
 )
 :
-    decompositionMethod(decompositionDict),
-    mesh_(mesh)
+    decompositionMethod(decompositionDict)
 {}
 
 
@@ -93,6 +91,7 @@ Foam::metisDecomp::metisDecomp
 
 Foam::labelList Foam::metisDecomp::decompose
 (
+    const polyMesh& mesh,
     const pointField& points,
     const scalarField& pointWeights
 )
@@ -112,6 +111,7 @@ Foam::labelList Foam::metisDecomp::decompose
 
 Foam::labelList Foam::metisDecomp::decompose
 (
+    const polyMesh& mesh,
     const labelList& agglom,
     const pointField& agglomPoints,
     const scalarField& agglomWeights
diff --git a/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C b/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C
index 2e877de4a9..9e59239fd1 100644
--- a/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C
+++ b/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C
@@ -48,7 +48,7 @@ namespace Foam
     (
         decompositionMethod,
         ptscotchDecomp,
-        dictionaryMesh
+        dictionary
     );
 }
 
@@ -85,12 +85,10 @@ Foam::label Foam::ptscotchDecomp::decompose
 
 Foam::ptscotchDecomp::ptscotchDecomp
 (
-    const dictionary& decompositionDict,
-    const polyMesh& mesh
+    const dictionary& decompositionDict
 )
 :
-    decompositionMethod(decompositionDict),
-    mesh_(mesh)
+    decompositionMethod(decompositionDict)
 {}
 
 
@@ -98,6 +96,7 @@ Foam::ptscotchDecomp::ptscotchDecomp
 
 Foam::labelList Foam::ptscotchDecomp::decompose
 (
+    const polyMesh& mesh,
     const pointField& points,
     const scalarField& pointWeights
 )
@@ -117,6 +116,7 @@ Foam::labelList Foam::ptscotchDecomp::decompose
 
 Foam::labelList Foam::ptscotchDecomp::decompose
 (
+    const polyMesh& mesh,
     const labelList& agglom,
     const pointField& agglomPoints,
     const scalarField& pointWeights
diff --git a/src/dummyThirdParty/scotchDecomp/dummyScotchDecomp.C b/src/dummyThirdParty/scotchDecomp/dummyScotchDecomp.C
index 3c84633932..007f708a5d 100644
--- a/src/dummyThirdParty/scotchDecomp/dummyScotchDecomp.C
+++ b/src/dummyThirdParty/scotchDecomp/dummyScotchDecomp.C
@@ -47,7 +47,7 @@ namespace Foam
     (
         decompositionMethod,
         scotchDecomp,
-        dictionaryMesh
+        dictionary
     );
 }
 
@@ -59,6 +59,7 @@ void Foam::scotchDecomp::check(const int retVal, const char* str)
 
 Foam::label Foam::scotchDecomp::decompose
 (
+    const fileName& meshPath,
     const List<int>& adjncy,
     const List<int>& xadj,
     const scalarField& cWeights,
@@ -68,13 +69,14 @@ Foam::label Foam::scotchDecomp::decompose
 {
     FatalErrorIn
     (
-        "label scotchDecomp::decompose"
-        "("
-            "const List<int>&, "
-            "const List<int>&, "
-            "const scalarField&, "
-            "List<int>&"
-        ")"
+        "label scotchDecomp::decompose\n"
+        "(\n"
+            "const fileName& meshPath,\n"
+            "const List<int>&,\n"
+            "const List<int>&,\n"
+            "const scalarField&,\n"
+            "List<int>&\n"
+        ")\n"
     )   << notImplementedMessage << exit(FatalError);
 
     return -1;
@@ -85,12 +87,10 @@ Foam::label Foam::scotchDecomp::decompose
 
 Foam::scotchDecomp::scotchDecomp
 (
-    const dictionary& decompositionDict,
-    const polyMesh& mesh
+    const dictionary& decompositionDict
 )
 :
-    decompositionMethod(decompositionDict),
-    mesh_(mesh)
+    decompositionMethod(decompositionDict)
 {}
 
 
@@ -98,6 +98,7 @@ Foam::scotchDecomp::scotchDecomp
 
 Foam::labelList Foam::scotchDecomp::decompose
 (
+    const polyMesh& mesh,
     const pointField& points,
     const scalarField& pointWeights
 )
@@ -117,6 +118,7 @@ Foam::labelList Foam::scotchDecomp::decompose
 
 Foam::labelList Foam::scotchDecomp::decompose
 (
+    const polyMesh& mesh,
     const labelList& agglom,
     const pointField& agglomPoints,
     const scalarField& pointWeights
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
index 1fd66dfda8..b5085a465a 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
@@ -415,11 +415,7 @@ Foam::autoHexMeshDriver::autoHexMeshDriver
 
     {
         // Decomposition
-        decomposerPtr_ = decompositionMethod::New
-        (
-            decomposeDict,
-            mesh_
-        );
+        decomposerPtr_ = decompositionMethod::New(decomposeDict);
         decompositionMethod& decomposer = decomposerPtr_();
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index 0622273a4e..43d1c74060 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -1214,6 +1214,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                 // Normal decomposition
                 distribution = decomposer.decompose
                 (
+                    mesh_,
                     mesh_.cellCentres(),
                     cellWeights
                 );
@@ -1224,6 +1225,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
             // Normal decomposition
             distribution = decomposer.decompose
             (
+                mesh_,
                 mesh_.cellCentres(),
                 cellWeights
             );
diff --git a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
index a046f05050..e736ec8945 100644
--- a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
@@ -34,6 +34,7 @@ License
 
 #include "IFstream.H"
 #include "decompositionMethod.H"
+#include "geomDecomp.H"
 #include "vectorList.H"
 #include "PackedBoolList.H"
 
@@ -855,6 +856,19 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
                 << " does not decompose in parallel."
                 << " Please choose one that does." << exit(FatalError);
         }
+
+        if (!isA<geomDecomp>(decomposer_()))
+        {
+            FatalErrorIn
+            (
+                "distributedTriSurfaceMesh::independentlyDistributedBbs"
+                "(const triSurface&)"
+            )   << "The decomposition method " << decomposer_().typeName
+                << " is not a geometric decomposition method." << endl
+                << "Only geometric decomposition methods are currently"
+                << " supported."
+                << exit(FatalError);
+        }
     }
 
     // Do decomposition according to triangle centre
@@ -864,8 +878,11 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
         triCentres[triI] = s[triI].centre(s.points());
     }
 
+
+    geomDecomp& decomposer = refCast<geomDecomp>(decomposer_());
+
     // Do the actual decomposition
-    labelList distribution(decomposer_->decompose(triCentres));
+    labelList distribution(decomposer.decompose(triCentres));
 
     // Find bounding box for all triangles on new distribution.
 
diff --git a/src/parallel/decompose/decompositionMethods/Make/files b/src/parallel/decompose/decompositionMethods/Make/files
index 43bb9a25c1..fccb8b7c8f 100644
--- a/src/parallel/decompose/decompositionMethods/Make/files
+++ b/src/parallel/decompose/decompositionMethods/Make/files
@@ -3,5 +3,6 @@ geomDecomp/geomDecomp.C
 simpleGeomDecomp/simpleGeomDecomp.C
 hierarchGeomDecomp/hierarchGeomDecomp.C
 manualDecomp/manualDecomp.C
+multiLevelDecomp/multiLevelDecomp.C
 
 LIB = $(FOAM_LIBBIN)/libdecompositionMethods
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
index 4e026f9e48..e68e31c7b9 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -37,7 +37,6 @@ namespace Foam
 {
     defineTypeNameAndDebug(decompositionMethod, 0);
     defineRunTimeSelectionTable(decompositionMethod, dictionary);
-    defineRunTimeSelectionTable(decompositionMethod, dictionaryMesh);
 }
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -71,57 +70,45 @@ Foam::autoPtr<Foam::decompositionMethod> Foam::decompositionMethod::New
 }
 
 
-Foam::autoPtr<Foam::decompositionMethod> Foam::decompositionMethod::New
-(
-    const dictionary& decompositionDict,
-    const polyMesh& mesh
-)
-{
-    const word methodType(decompositionDict.lookup("method"));
-
-    Info<< "Selecting decompositionMethod " << methodType << endl;
-
-    dictionaryMeshConstructorTable::iterator cstrIter =
-        dictionaryMeshConstructorTablePtr_->find(methodType);
-
-    if (cstrIter == dictionaryMeshConstructorTablePtr_->end())
-    {
-        FatalErrorIn
-        (
-            "decompositionMethod::New"
-            "(const dictionary& decompositionDict, "
-            "const polyMesh& mesh)"
-        )   << "Unknown decompositionMethod "
-            << methodType << nl << nl
-            << "Valid decompositionMethods are : " << endl
-            << dictionaryMeshConstructorTablePtr_->sortedToc()
-            << exit(FatalError);
-    }
-
-    return autoPtr<decompositionMethod>(cstrIter()(decompositionDict, mesh));
-}
-
-
 Foam::labelList Foam::decompositionMethod::decompose
 (
+    const polyMesh& mesh,
     const pointField& points
 )
 {
-    scalarField weights(0);
+    scalarField weights(points.size(), 1.0);
 
-    return decompose(points, weights);
+    return decompose(mesh, points, weights);
 }
 
 
 Foam::labelList Foam::decompositionMethod::decompose
 (
+    const polyMesh& mesh,
     const labelList& fineToCoarse,
     const pointField& coarsePoints,
     const scalarField& coarseWeights
 )
 {
+    CompactListList<label> coarseCellCells;
+    calcCellCells
+    (
+        mesh,
+        fineToCoarse,
+        coarsePoints.size(),
+        coarseCellCells
+    );
+
     // Decompose based on agglomerated points
-    labelList coarseDistribution(decompose(coarsePoints, coarseWeights));
+    labelList coarseDistribution
+    (
+        decompose
+        (
+            coarseCellCells(),
+            coarsePoints,
+            coarseWeights
+        )
+    );
 
     // Rework back into decomposition for original mesh_
     labelList fineDistribution(fineToCoarse.size());
@@ -137,22 +124,20 @@ Foam::labelList Foam::decompositionMethod::decompose
 
 Foam::labelList Foam::decompositionMethod::decompose
 (
+    const polyMesh& mesh,
     const labelList& fineToCoarse,
     const pointField& coarsePoints
 )
 {
-    // Decompose based on agglomerated points
-    labelList coarseDistribution(decompose(coarsePoints));
+    scalarField cWeights(coarsePoints.size(), 1.0);
 
-    // Rework back into decomposition for original mesh_
-    labelList fineDistribution(fineToCoarse.size());
-
-    forAll(fineDistribution, i)
-    {
-        fineDistribution[i] = coarseDistribution[fineToCoarse[i]];
-    }
-
-    return fineDistribution;
+    return decompose
+    (
+        mesh,
+        fineToCoarse,
+        coarsePoints,
+        cWeights
+    );
 }
 
 
@@ -162,57 +147,12 @@ Foam::labelList Foam::decompositionMethod::decompose
     const pointField& cc
 )
 {
-    scalarField cWeights(0);
+    scalarField cWeights(cc.size(), 1.0);
 
     return decompose(globalCellCells, cc, cWeights);
 }
 
 
-void Foam::decompositionMethod::calcCellCells
-(
-    const polyMesh& mesh,
-    const labelList& fineToCoarse,
-    const label nCoarse,
-    labelListList& cellCells
-)
-{
-    if (fineToCoarse.size() != mesh.nCells())
-    {
-        FatalErrorIn
-        (
-            "decompositionMethod::calcCellCells"
-            "(const labelList&, labelListList&) const"
-        )   << "Only valid for mesh agglomeration." << exit(FatalError);
-    }
-
-    List<DynamicList<label> > dynCellCells(nCoarse);
-
-    forAll(mesh.faceNeighbour(), faceI)
-    {
-        label own = fineToCoarse[mesh.faceOwner()[faceI]];
-        label nei = fineToCoarse[mesh.faceNeighbour()[faceI]];
-
-        if (own != nei)
-        {
-            if (findIndex(dynCellCells[own], nei) == -1)
-            {
-                dynCellCells[own].append(nei);
-            }
-            if (findIndex(dynCellCells[nei], own) == -1)
-            {
-                dynCellCells[nei].append(own);
-            }
-        }
-    }
-
-    cellCells.setSize(dynCellCells.size());
-    forAll(dynCellCells, coarseI)
-    {
-        cellCells[coarseI].transfer(dynCellCells[coarseI]);
-    }
-}
-
-
 // Return the minimum face between two cells. Only relevant for
 // cells with multiple faces inbetween.
 Foam::label Foam::decompositionMethod::masterFace
@@ -250,208 +190,209 @@ Foam::label Foam::decompositionMethod::masterFace
 }
 
 
-void Foam::decompositionMethod::calcCSR
+//void Foam::decompositionMethod::calcCSR
+//(
+//    const polyMesh& mesh,
+//    List<int>& adjncy,
+//    List<int>& xadj
+//)
+//{
+//    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+//
+//    // Make Metis CSR (Compressed Storage Format) storage
+//    //   adjncy      : contains neighbours (= edges in graph)
+//    //   xadj(celli) : start of information in adjncy for celli
+//
+//
+//    // Count unique faces between cells
+//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+//
+//    labelList nFacesPerCell(mesh.nCells(), 0);
+//
+//    // Internal faces
+//    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+//    {
+//        label own = mesh.faceOwner()[faceI];
+//        label nei = mesh.faceNeighbour()[faceI];
+//
+//        if (faceI == masterFace(mesh, own, nei))
+//        {
+//            nFacesPerCell[own]++;
+//            nFacesPerCell[nei]++;
+//        }
+//    }
+//
+//    // Coupled faces. Only cyclics done.
+//    HashSet<edge, Hash<edge> > cellPair(mesh.nFaces()-mesh.nInternalFaces());
+//
+//    forAll(pbm, patchI)
+//    {
+//        if
+//        (
+//            isA<cyclicPolyPatch>(pbm[patchI])
+//         && refCast<const cyclicPolyPatch>(pbm[patchI]).owner()
+//        )
+//        {
+//            const cyclicPolyPatch& cycPatch = refCast<const cyclicPolyPatch>
+//            (
+//                pbm[patchI]
+//            );
+//
+//            const unallocLabelList& faceCells = cycPatch.faceCells();
+//            const unallocLabelList& nbrCells =
+//                cycPatch.neighbPatch().faceCells();
+//
+//            forAll(faceCells, facei)
+//            {
+//                label own = faceCells[facei];
+//                label nei = nbrCells[facei];
+//
+//                if (cellPair.insert(edge(own, nei)))
+//                {
+//                    nFacesPerCell[own]++;
+//                    nFacesPerCell[nei]++;
+//                }
+//            }
+//        }
+//    }
+//
+//
+//    // Size tables
+//    // ~~~~~~~~~~~
+//
+//    // Sum nFacesPerCell
+//    xadj.setSize(mesh.nCells()+1);
+//
+//    label nConnections = 0;
+//
+//    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
+//    {
+//        xadj[cellI] = nConnections;
+//        nConnections += nFacesPerCell[cellI];
+//    }
+//    xadj[mesh.nCells()] = nConnections;
+//    adjncy.setSize(nConnections);
+//
+//
+//
+//    // Fill tables
+//    // ~~~~~~~~~~~
+//
+//    nFacesPerCell = 0;
+//
+//    // Internal faces
+//    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+//    {
+//        label own = mesh.faceOwner()[faceI];
+//        label nei = mesh.faceNeighbour()[faceI];
+//
+//        if (faceI == masterFace(mesh, own, nei))
+//        {
+//            adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
+//            adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+//        }
+//    }
+//
+//    // Coupled faces. Only cyclics done.
+//    cellPair.clear();
+//    forAll(pbm, patchI)
+//    {
+//        if
+//        (
+//            isA<cyclicPolyPatch>(pbm[patchI])
+//         && refCast<const cyclicPolyPatch>(pbm[patchI]).owner()
+//        )
+//        {
+//            const cyclicPolyPatch& cycPatch = refCast<const cyclicPolyPatch>
+//            (
+//                pbm[patchI]
+//            );
+//
+//            const unallocLabelList& faceCells = cycPatch.faceCells();
+//            const unallocLabelList& nbrCells =
+//                cycPatch.neighbPatch().faceCells();
+//
+//            forAll(faceCells, facei)
+//            {
+//                label own = faceCells[facei];
+//                label nei = nbrCells[facei];
+//
+//                if (cellPair.insert(edge(own, nei)))
+//                {
+//                    adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
+//                    adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+//                }
+//            }
+//        }
+//    }
+//}
+//
+//
+//// From cell-cell connections to Metis format (like CompactListList)
+//void Foam::decompositionMethod::calcCSR
+//(
+//    const labelListList& cellCells,
+//    List<int>& adjncy,
+//    List<int>& xadj
+//)
+//{
+//    labelHashSet nbrCells;
+//
+//    // Count number of internal faces
+//    label nConnections = 0;
+//
+//    forAll(cellCells, coarseI)
+//    {
+//        nbrCells.clear();
+//
+//        const labelList& cCells = cellCells[coarseI];
+//
+//        forAll(cCells, i)
+//        {
+//            if (nbrCells.insert(cCells[i]))
+//            {
+//                nConnections++;
+//            }
+//        }
+//    }
+//
+//    // Create the adjncy array as twice the size of the total number of
+//    // internal faces
+//    adjncy.setSize(nConnections);
+//
+//    xadj.setSize(cellCells.size()+1);
+//
+//
+//    // Fill in xadj
+//    // ~~~~~~~~~~~~
+//    label freeAdj = 0;
+//
+//    forAll(cellCells, coarseI)
+//    {
+//        xadj[coarseI] = freeAdj;
+//
+//        nbrCells.clear();
+//
+//        const labelList& cCells = cellCells[coarseI];
+//
+//        forAll(cCells, i)
+//        {
+//            if (nbrCells.insert(cCells[i]))
+//            {
+//                adjncy[freeAdj++] = cCells[i];
+//            }
+//        }
+//    }
+//    xadj[cellCells.size()] = freeAdj;
+//}
+
+
+void Foam::decompositionMethod::calcCellCells
 (
     const polyMesh& mesh,
-    List<int>& adjncy,
-    List<int>& xadj
-)
-{
-    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
-
-    // Make Metis CSR (Compressed Storage Format) storage
-    //   adjncy      : contains neighbours (= edges in graph)
-    //   xadj(celli) : start of information in adjncy for celli
-
-
-    // Count unique faces between cells
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelList nFacesPerCell(mesh.nCells(), 0);
-
-    // Internal faces
-    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
-    {
-        label own = mesh.faceOwner()[faceI];
-        label nei = mesh.faceNeighbour()[faceI];
-
-        if (faceI == masterFace(mesh, own, nei))
-        {
-            nFacesPerCell[own]++;
-            nFacesPerCell[nei]++;
-        }
-    }
-
-    // Coupled faces. Only cyclics done.
-    HashSet<edge, Hash<edge> > cellPair(mesh.nFaces()-mesh.nInternalFaces());
-
-    forAll(pbm, patchI)
-    {
-        if
-        (
-            isA<cyclicPolyPatch>(pbm[patchI])
-         && refCast<const cyclicPolyPatch>(pbm[patchI]).owner()
-        )
-        {
-            const cyclicPolyPatch& cycPatch = refCast<const cyclicPolyPatch>
-            (
-                pbm[patchI]
-            );
-
-            const unallocLabelList& faceCells = cycPatch.faceCells();
-            const unallocLabelList& nbrCells =
-                cycPatch.neighbPatch().faceCells();
-
-            forAll(faceCells, facei)
-            {
-                label own = faceCells[facei];
-                label nei = nbrCells[facei];
-
-                if (cellPair.insert(edge(own, nei)))
-                {
-                    nFacesPerCell[own]++;
-                    nFacesPerCell[nei]++;
-                }
-            }
-        }
-    }
-
-
-    // Size tables
-    // ~~~~~~~~~~~
-
-    // Sum nFacesPerCell
-    xadj.setSize(mesh.nCells()+1);
-
-    label nConnections = 0;
-
-    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
-    {
-        xadj[cellI] = nConnections;
-        nConnections += nFacesPerCell[cellI];
-    }
-    xadj[mesh.nCells()] = nConnections;
-    adjncy.setSize(nConnections);
-
-
-
-    // Fill tables
-    // ~~~~~~~~~~~
-
-    nFacesPerCell = 0;
-
-    // Internal faces
-    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
-    {
-        label own = mesh.faceOwner()[faceI];
-        label nei = mesh.faceNeighbour()[faceI];
-
-        if (faceI == masterFace(mesh, own, nei))
-        {
-            adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
-            adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
-        }
-    }
-
-    // Coupled faces. Only cyclics done.
-    cellPair.clear();
-    forAll(pbm, patchI)
-    {
-        if
-        (
-            isA<cyclicPolyPatch>(pbm[patchI])
-         && refCast<const cyclicPolyPatch>(pbm[patchI]).owner()
-        )
-        {
-            const cyclicPolyPatch& cycPatch = refCast<const cyclicPolyPatch>
-            (
-                pbm[patchI]
-            );
-
-            const unallocLabelList& faceCells = cycPatch.faceCells();
-            const unallocLabelList& nbrCells =
-                cycPatch.neighbPatch().faceCells();
-
-            forAll(faceCells, facei)
-            {
-                label own = faceCells[facei];
-                label nei = nbrCells[facei];
-
-                if (cellPair.insert(edge(own, nei)))
-                {
-                    adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
-                    adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
-                }
-            }
-        }
-    }
-}
-
-
-// From cell-cell connections to Metis format (like CompactListList)
-void Foam::decompositionMethod::calcCSR
-(
-    const labelListList& cellCells,
-    List<int>& adjncy,
-    List<int>& xadj
-)
-{
-    labelHashSet nbrCells;
-
-    // Count number of internal faces
-    label nConnections = 0;
-
-    forAll(cellCells, coarseI)
-    {
-        nbrCells.clear();
-
-        const labelList& cCells = cellCells[coarseI];
-
-        forAll(cCells, i)
-        {
-            if (nbrCells.insert(cCells[i]))
-            {
-                nConnections++;
-            }
-        }
-    }
-
-    // Create the adjncy array as twice the size of the total number of
-    // internal faces
-    adjncy.setSize(nConnections);
-
-    xadj.setSize(cellCells.size()+1);
-
-
-    // Fill in xadj
-    // ~~~~~~~~~~~~
-    label freeAdj = 0;
-
-    forAll(cellCells, coarseI)
-    {
-        xadj[coarseI] = freeAdj;
-
-        nbrCells.clear();
-
-        const labelList& cCells = cellCells[coarseI];
-
-        forAll(cCells, i)
-        {
-            if (nbrCells.insert(cCells[i]))
-            {
-                adjncy[freeAdj++] = cCells[i];
-            }
-        }
-    }
-    xadj[cellCells.size()] = freeAdj;
-}
-
-
-void Foam::decompositionMethod::calcDistributedCSR
-(
-    const polyMesh& mesh,
-    List<int>& adjncy,
-    List<int>& xadj
+    const labelList& agglom,
+    const label nCoarse,
+    CompactListList<label>& cellCells
 )
 {
     // Create global cell numbers
@@ -459,14 +400,6 @@ void Foam::decompositionMethod::calcDistributedCSR
 
     globalIndex globalCells(mesh.nCells());
 
-
-    //
-    // Make Metis Distributed CSR (Compressed Storage Format) storage
-    //   adjncy      : contains cellCells (= edges in graph)
-    //   xadj(celli) : start of information in adjncy for celli
-    //
-
-
     const labelList& faceOwner = mesh.faceOwner();
     const labelList& faceNeighbour = mesh.faceNeighbour();
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
@@ -475,7 +408,7 @@ void Foam::decompositionMethod::calcDistributedCSR
     // Get renumbered owner on other side of coupled faces
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    List<int> globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
+    labelList globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
 
     forAll(patches, patchI)
     {
@@ -504,18 +437,15 @@ void Foam::decompositionMethod::calcDistributedCSR
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Number of faces per cell
-    List<int> nFacesPerCell(mesh.nCells(), 0);
+    labelList nFacesPerCell(mesh.nCells(), 0);
 
     for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
     {
-        label own = faceOwner[faceI];
-        label nei = faceNeighbour[faceI];
+        label own = agglom[faceOwner[faceI]];
+        label nei = agglom[faceNeighbour[faceI]];
 
-        if (faceI == masterFace(mesh, own, nei))
-        {
-            nFacesPerCell[own]++;
-            nFacesPerCell[nei]++;
-        }
+        nFacesPerCell[own]++;
+        nFacesPerCell[nei]++;
     }
 
     // Handle coupled faces
@@ -532,7 +462,7 @@ void Foam::decompositionMethod::calcDistributedCSR
 
             forAll(pp, i)
             {
-                label own =  faceOwner[faceI];
+                label own = agglom[faceOwner[faceI]];
                 label globalNei = globalNeighbour[bFaceI];
                 if (cellPair.insert(edge(own, globalNei)))
                 {
@@ -545,43 +475,24 @@ void Foam::decompositionMethod::calcDistributedCSR
     }
 
 
-    // Fill in xadj
-    // ~~~~~~~~~~~~
+    // Fill in offset and data
+    // ~~~~~~~~~~~~~~~~~~~~~~~
 
-    xadj.setSize(mesh.nCells()+1);
-
-    int freeAdj = 0;
-
-    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
-    {
-        xadj[cellI] = freeAdj;
-
-        freeAdj += nFacesPerCell[cellI];
-    }
-    xadj[mesh.nCells()] = freeAdj;
-
-
-
-    // Fill in adjncy
-    // ~~~~~~~~~~~~~~
-
-    adjncy.setSize(freeAdj);
+    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++)
     {
-        label own = faceOwner[faceI];
-        label nei = faceNeighbour[faceI];
+        label own = agglom[faceOwner[faceI]];
+        label nei = agglom[faceNeighbour[faceI]];
 
-        if (faceI == masterFace(mesh, own, nei))
-        {
-            adjncy[xadj[own] + nFacesPerCell[own]++] =
-                globalCells.toGlobal(nei);
-            adjncy[xadj[nei] + nFacesPerCell[nei]++] =
-                globalCells.toGlobal(own);
-        }
+        m[offsets[own] + nFacesPerCell[own]++] = globalCells.toGlobal(nei);
+        m[offsets[nei] + nFacesPerCell[nei]++] = globalCells.toGlobal(own);
     }
 
     // For boundary faces is offsetted coupled neighbour
@@ -597,17 +508,41 @@ void Foam::decompositionMethod::calcDistributedCSR
 
             forAll(pp, i)
             {
-                label own = faceOwner[faceI];
+                label own = agglom[faceOwner[faceI]];
                 label globalNei = globalNeighbour[bFaceI];
                 if (cellPair.insert(edge(own, globalNei)))
                 {
-                    adjncy[xadj[own] + nFacesPerCell[own]++] = globalNei;
+                    m[offsets[own] + nFacesPerCell[own]++] = globalNei;
                 }
                 faceI++;
                 bFaceI++;
             }
         }
     }
+
+
+    // Check for duplicates connections between cells as a postprocessing step
+    // (since quite rare)
+    label startIndex = 0;
+    label newIndex = 0;
+    labelHashSet nbrCells;
+    forAll(cellCells, cellI)
+    {
+        nbrCells.clear();
+
+        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;
+        endIndex = newIndex;
+    }
 }
 
 
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H
index dabca144c8..26ad6cc5e7 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H
@@ -37,6 +37,7 @@ SourceFiles
 
 #include "polyMesh.H"
 #include "pointField.H"
+#include "CompactListList.H"
 
 namespace Foam
 {
@@ -56,43 +57,37 @@ protected:
         label nProcessors_;
 
 
-        //- Helper: determine (non-parallel) cellCells from mesh agglomeration.
+        //- Helper: determine (global) cellCells from mesh agglomeration.
         static void calcCellCells
         (
             const polyMesh& mesh,
             const labelList& agglom,
             const label nCoarse,
-            labelListList& cellCells
+            CompactListList<label>& cellCells
         );
 
         //- Calculate the minimum face between two neighbouring cells
         //  (usually there is only one)
         static label masterFace(const polyMesh&, const label, const label);
 
-        // From mesh to compact row storage format
-        // (like CompactListList)
-        static void calcCSR
-        (
-            const polyMesh& mesh,
-            List<int>& adjncy,
-            List<int>& xadj
-        );
+//        // From mesh to compact row storage format
+//        // (like CompactListList)
+//        static void calcCSR
+//        (
+//            const polyMesh& mesh,
+//            List<int>& adjncy,
+//            List<int>& xadj
+//        );
+//
+//        // From cell-cell connections to compact row storage format
+//        // (like CompactListList)
+//        static void calcCSR
+//        (
+//            const labelListList& cellCells,
+//            List<int>& adjncy,
+//            List<int>& xadj
+//        );
 
-        // From cell-cell connections to compact row storage format
-        // (like CompactListList)
-        static void calcCSR
-        (
-            const labelListList& cellCells,
-            List<int>& adjncy,
-            List<int>& xadj
-        );
-
-        static void calcDistributedCSR
-        (
-            const polyMesh& mesh,
-            List<int>& adjncy,
-            List<int>& xadj
-        );
 
 private:
 
@@ -122,18 +117,6 @@ public:
             (decompositionDict)
         );
 
-        declareRunTimeSelectionTable
-        (
-            autoPtr,
-            decompositionMethod,
-            dictionaryMesh,
-            (
-                const dictionary& decompositionDict,
-                const polyMesh& mesh
-            ),
-            (decompositionDict, mesh)
-        );
-
 
     // Selectors
 
@@ -143,13 +126,6 @@ public:
             const dictionary& decompositionDict
         );
 
-        //- Return a reference to the selected decomposition method
-        static autoPtr<decompositionMethod> New
-        (
-            const dictionary& decompositionDict,
-            const polyMesh& mesh
-        );
-
 
     // Constructors
 
@@ -171,63 +147,105 @@ public:
 
     // Member Functions
 
+        label nDomains() const
+        {
+            return nProcessors_;
+        }
+
         //- Is method parallel aware (i.e. does it synchronize domains across
         //  proc boundaries)
         virtual bool parallelAware() const = 0;
 
-        //- Return for every coordinate the wanted processor number. Use the
-        //  mesh connectivity (if needed)
-        virtual labelList decompose
-        (
-            const pointField& points,
-            const scalarField& pointWeights
-        ) = 0;
 
-        //- Like decompose but with uniform weights on the points
-        virtual labelList decompose(const pointField&);
+        // No topology (implemented by geometric decomposers)
+
+            //- Return for every coordinate the wanted processor number.
+            virtual labelList decompose
+            (
+                const pointField& points,
+                const scalarField& pointWeights
+            )
+            {
+                notImplemented
+                (
+                    "decompositionMethod:decompose(const pointField&"
+                    ", const scalarField&)"
+                );
+                return labelList(0);
+            }
+
+            //- Like decompose but with uniform weights on the points
+            virtual labelList decompose(const pointField&)
+            {
+                notImplemented
+                (
+                    "decompositionMethod:decompose(const pointField&)"
+                );
+                return labelList(0);
+            }
 
 
-        //- 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
-        //  processors use the globalCellCells instead.
-        virtual labelList decompose
-        (
-            const labelList& cellToRegion,
-            const pointField& regionPoints,
-            const scalarField& regionWeights
-        );
+        // Topology provided by mesh
 
-        //- Like decompose but with uniform weights on the regions
-        virtual labelList decompose
-        (
-            const labelList& cellToRegion,
-            const pointField& regionPoints
-        );
+            //- Return for every coordinate the wanted processor number. Use the
+            //  mesh connectivity (if needed)
+            virtual labelList decompose
+            (
+                const polyMesh& mesh,
+                const pointField& points,
+                const scalarField& pointWeights
+            ) = 0;
+
+            //- Like decompose but with uniform weights on the points
+            virtual labelList decompose(const polyMesh&, const pointField&);
 
 
-        //- 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
-        virtual labelList decompose
-        (
-            const labelListList& globalCellCells,
-            const pointField& cc,
-            const scalarField& cWeights
-        ) = 0;
+            //- 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
+            //  processors use the globalCellCells instead.
+            virtual labelList decompose
+            (
+                const polyMesh& mesh,
+                const labelList& cellToRegion,
+                const pointField& regionPoints,
+                const scalarField& regionWeights
+            );
 
-        //- Like decompose but with uniform weights on the cells
-        virtual labelList decompose
-        (
-            const labelListList& globalCellCells,
-            const pointField& cc
-        );
+            //- Like decompose but with uniform weights on the regions
+            virtual labelList decompose
+            (
+                const polyMesh& mesh,
+                const labelList& cellToRegion,
+                const pointField& regionPoints
+            );
+
+
+        // Topology provided explicitly addressing
+
+            //- 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 labelListList& globalCellCells,
+                const pointField& cc,
+                const scalarField& cWeights
+            ) = 0;
+
+            //- Like decompose but with uniform weights on the cells
+            virtual labelList decompose
+            (
+                const labelListList& globalCellCells,
+                const pointField& cc
+            );
 
 };
 
diff --git a/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.H b/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.H
index 14dfc5a1c1..f10b5f7c6e 100644
--- a/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.H
@@ -71,6 +71,17 @@ public:
             const dictionary& decompositionDict,
             const word& derivedType
         );
+
+        //- Return for every coordinate the wanted processor number.
+        virtual labelList decompose
+        (
+            const pointField& points,
+            const scalarField& pointWeights
+        ) = 0;
+
+        //- Like decompose but with uniform weights on the points
+        virtual labelList decompose(const pointField&) = 0;
+
 };
 
 
diff --git a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
index 98f8081608..3761d11b7f 100644
--- a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
@@ -40,13 +40,6 @@ namespace Foam
         hierarchGeomDecomp,
         dictionary
     );
-
-    addToRunTimeSelectionTable
-    (
-        decompositionMethod,
-        hierarchGeomDecomp,
-        dictionaryMesh
-    );
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -685,35 +678,6 @@ void Foam::hierarchGeomDecomp::sortComponent
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::hierarchGeomDecomp::hierarchGeomDecomp
-(
-    const dictionary& decompositionDict
-)
-:
-    geomDecomp(decompositionDict, typeName),
-    decompOrder_()
-{
-    setDecompOrder();
-}
-
-
-Foam::hierarchGeomDecomp::hierarchGeomDecomp
-(
-    const dictionary& decompositionDict,
-    const polyMesh&
-)
-:
-    geomDecomp(decompositionDict, hierarchGeomDecomp::typeName),
-    decompOrder_()
-{
-    setDecompOrder();
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
 Foam::labelList Foam::hierarchGeomDecomp::decompose
 (
     const pointField& points
@@ -797,4 +761,20 @@ Foam::labelList Foam::hierarchGeomDecomp::decompose
 }
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::hierarchGeomDecomp::hierarchGeomDecomp
+(
+    const dictionary& decompositionDict
+)
+:
+    geomDecomp(decompositionDict, typeName),
+    decompOrder_()
+{
+    setDecompOrder();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 // ************************************************************************* //
diff --git a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
index 62eef72ae3..b82b9890da 100644
--- a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
@@ -154,7 +154,6 @@ class hierarchGeomDecomp
             labelList& finalDecomp          // overall decomposition
         );
 
-
         //- Disallow default bitwise copy construct and assignment
         void operator=(const hierarchGeomDecomp&);
         hierarchGeomDecomp(const hierarchGeomDecomp&);
@@ -171,13 +170,6 @@ public:
         //- Construct given the decomposition dictionary
         hierarchGeomDecomp(const dictionary& decompositionDict);
 
-        //- Construct given the decomposition dictionary and mesh
-        hierarchGeomDecomp
-        (
-            const dictionary& decompositionDict,
-            const polyMesh& mesh
-        );
-
 
     //- Destructor
     virtual ~hierarchGeomDecomp()
@@ -192,8 +184,7 @@ public:
             return true;
         }
 
-        //- Return for every coordinate the wanted processor number. Use the
-        //  mesh connectivity (if needed)
+        //- Return for every coordinate the wanted processor number.
         virtual labelList decompose
         (
             const pointField&,
@@ -204,6 +195,26 @@ public:
         //  so kept separate for now.
         virtual labelList decompose(const pointField&);
 
+
+        //- Return for every coordinate the wanted processor number. Use the
+        //  mesh connectivity (if needed)
+        virtual labelList decompose
+        (
+            const polyMesh& mesh,
+            const pointField& cc,
+            const scalarField& cWeights
+        )
+        {
+            return decompose(cc, cWeights);
+        }
+
+        //- Without weights. Code for weighted decomposition is a bit complex
+        //  so kept separate for now.
+        virtual labelList decompose(const polyMesh& mesh, const pointField& cc)
+        {
+            return decompose(cc);
+        }
+
         //- Return for every coordinate the wanted processor number. Explicitly
         //  provided connectivity - does not use mesh_.
         //  The connectivity is equal to mesh.cellCells() except for
diff --git a/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.C b/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.C
index 4c51c2a591..b0a1377e44 100644
--- a/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.C
@@ -43,34 +43,14 @@ namespace Foam
         manualDecomp,
         dictionary
     );
-
-    addToRunTimeSelectionTable
-    (
-        decompositionMethod,
-        manualDecomp,
-        dictionaryMesh
-    );
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::manualDecomp::manualDecomp(const dictionary& decompositionDict)
-:
-    decompositionMethod(decompositionDict)
-{
-    notImplemented("manualDecomp(const dictionary&)");
-}
-
-
-Foam::manualDecomp::manualDecomp
-(
-    const dictionary& decompositionDict,
-    const polyMesh& mesh
-)
 :
     decompositionMethod(decompositionDict),
-    meshPtr_(&mesh),
     decompDataFile_
     (
         decompositionDict.subDict(word(decompositionDict.lookup("method"))
@@ -83,12 +63,11 @@ Foam::manualDecomp::manualDecomp
 
 Foam::labelList Foam::manualDecomp::decompose
 (
+    const polyMesh& mesh,
     const pointField& points,
     const scalarField& pointWeights
 )
 {
-    const polyMesh& mesh = *meshPtr_;
-
     labelIOList finalDecomp
     (
         IOobject
diff --git a/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.H b/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.H
index e3da10ccac..9976b9c08d 100644
--- a/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/manualDecomp/manualDecomp.H
@@ -50,8 +50,6 @@ class manualDecomp
 {
     // Private data
 
-        const polyMesh* meshPtr_;
-
         fileName decompDataFile_;
 
 
@@ -73,14 +71,6 @@ public:
         //- Construct given the decomposition dictionary
         manualDecomp(const dictionary& decompositionDict);
 
-        //- Construct given the decomposition dictionary and mesh
-        manualDecomp
-        (
-            const dictionary& decompositionDict,
-            const polyMesh& mesh
-        );
-
-
     //- Destructor
     virtual ~manualDecomp()
     {}
@@ -99,8 +89,9 @@ public:
         //  mesh connectivity (if needed)
         virtual labelList decompose
         (
-            const pointField& points,
-            const scalarField& pointWeights
+            const polyMesh& mesh,
+            const pointField& cc,
+            const scalarField& cWeights
         );
 
         //- Return for every coordinate the wanted processor number. Explicitly
@@ -117,8 +108,14 @@ public:
             const scalarField& cWeights
         )
         {
-            return decompose(cc, cWeights);
+            notImplemented
+            (
+                "decompose(const labelListList&, const pointField&"
+                ", const scalarField)"
+            );
+            return labelList(0);
         }
+
 };
 
 
diff --git a/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C b/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C
new file mode 100644
index 0000000000..e9a61c729a
--- /dev/null
+++ b/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C
@@ -0,0 +1,293 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "multiLevelDecomp.H"
+#include "addToRunTimeSelectionTable.H"
+#include "IFstream.H"
+#include "globalIndex.H"
+#include "mapDistribute.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(multiLevelDecomp, 0);
+
+    addToRunTimeSelectionTable
+    (
+        decompositionMethod,
+        multiLevelDecomp,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Given a subset of cells determine the new global indices. The problem
+// is in the cells from neighbouring processors which need to be renumbered.
+void Foam::multiLevelDecomp::subsetGlobalCellCells
+(
+    const labelListList& cellCells,
+    const labelList& set,
+    labelListList& subCellCells,
+    label& cutConnections
+) const
+{
+    // Determine new index for cells by inverting subset
+    labelList oldToNew(invert(cellCells.size(), set));
+
+    globalIndex globalCells(cellCells.size());
+
+    // Subset locally the elements for which I have data
+    subCellCells = UIndirectList<labelList>(cellCells, set);
+
+    // Get new indices for neighbouring processors
+    List<Map<label> > compactMap;
+    mapDistribute map(globalCells, subCellCells, compactMap);
+    map.distribute(oldToNew);
+
+    // Now subCellCells contains indices into oldToNew which are the
+    // new locations of the neighbouring cells.
+
+    cutConnections = 0;
+
+    forAll(subCellCells, subCellI)
+    {
+        labelList& cCells = subCellCells[subCellI];
+
+        // Keep the connections to valid mapped cells
+        label newI = 0;
+        forAll(cCells, i)
+        {
+            label subCellI = oldToNew[cCells[i]];
+            if (subCellI == -1)
+            {
+                cutConnections++;
+            }
+            else
+            {
+                cCells[newI++] = subCellI;
+            }
+        }
+        cCells.setSize(newI);
+    }
+}
+
+
+void Foam::multiLevelDecomp::decompose
+(
+    const labelListList& pointPoints,
+    const pointField& points,
+    const scalarField& pointWeights,
+    const labelList& pointMap,      // map back to original points
+    const label levelI,
+
+    labelField& finalDecomp
+)
+{
+    labelList dist
+    (
+        methods_[levelI].decompose
+        (
+            pointPoints,
+            points,
+            pointWeights
+        )
+    );
+
+//Pout<< "At level " << levelI << endl;
+//forAll(dist, i)
+//{
+//    Pout<< "    " << i << " at:" << points[i]
+//        << " original:" << pointMap[i] << "  dist:" << dist[i]
+//        << " connected to:" << pointField(points, pointPoints[i])
+//        << endl;
+//}
+//Pout<< endl;
+
+    forAll(pointMap, i)
+    {
+        label orig = pointMap[i];
+        finalDecomp[orig] += dist[i];
+    }
+
+    if (levelI != methods_.size()-1)
+    {
+        // Recurse
+
+        // Determine points per domain
+        label n = methods_[levelI].nDomains();
+        labelListList domainToPoints(invertOneToMany(n, dist));
+
+        // 'Make space' for new levels of decomposition
+        finalDecomp *= methods_[levelI+1].nDomains();
+
+        // Extract processor+local index from point-point addressing
+
+        forAll(domainToPoints, domainI)
+        {
+            const labelList& myPoints = domainToPoints[domainI];
+
+            // Subset point-wise data.
+            pointField subPoints(points, myPoints);
+            scalarField subWeights(pointWeights, myPoints);
+            labelList subPointMap(UIndirectList<label>(pointMap, myPoints));
+            // Subset point-point addressing (adapt global numbering)
+            labelListList subPointPoints;
+            label nOutsideConnections;
+            subsetGlobalCellCells
+            (
+                pointPoints,
+                myPoints,
+                subPointPoints,
+                nOutsideConnections
+            );
+
+            //Info<< "At level " << levelI << "  domain " << domainI
+            //    << " have connections to other domains "
+            //    << returnReduce(nOutsideConnections, sumOp<label>())
+            //    << endl;
+
+            decompose
+            (
+                subPointPoints,
+                subPoints,
+                subWeights,
+                subPointMap,
+                levelI+1,
+
+                finalDecomp
+            );
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::multiLevelDecomp::multiLevelDecomp(const dictionary& decompositionDict)
+:
+    decompositionMethod(decompositionDict)
+{
+    const dictionary& myDict = decompositionDict_.subDict(typeName + "Coeffs");
+
+    methods_.setSize(myDict.size());
+    label i = 0;
+    forAllConstIter(dictionary, myDict, iter)
+    {
+        methods_.set(i++, decompositionMethod::New(iter().dict()));
+    }
+
+    label n = 1;
+    Info<< "decompositionMethod " << type() << " :" << endl;
+    forAll(methods_, i)
+    {
+        Info<< "    level " << i << " decomposing with " << methods_[i].type()
+            << " into " << methods_[i].nDomains() << " subdomains." << endl;
+
+        n *= methods_[i].nDomains();
+    }
+
+    if (n != nDomains())
+    {
+        FatalErrorIn("multiLevelDecomp::multiLevelDecomp(const dictionary&)")
+            << "Top level decomposition specifies " << nDomains()
+            << " domains which is not equal to the product of"
+            << " all sub domains " << n
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::multiLevelDecomp::parallelAware() const
+{
+    forAll(methods_, i)
+    {
+        if (!methods_[i].parallelAware())
+        {
+            return false;
+        }
+    }
+    return true;
+}
+
+
+Foam::labelList Foam::multiLevelDecomp::decompose
+(
+    const polyMesh& mesh,
+    const pointField& cc,
+    const scalarField& cWeights
+)
+{
+    CompactListList<label> cellCells;
+    calcCellCells(mesh, identity(cc.size()), cc.size(), cellCells);
+
+    labelField finalDecomp(cc.size(), 0);
+    labelList cellMap(identity(cc.size()));
+
+    decompose
+    (
+        cellCells(),
+        cc,
+        cWeights,
+        cellMap,      // map back to original cells
+        0,
+
+        finalDecomp
+    );
+
+    return finalDecomp;
+}
+
+
+Foam::labelList Foam::multiLevelDecomp::decompose
+(
+    const labelListList& globalPointPoints,
+    const pointField& points,
+    const scalarField& pointWeights
+)
+{
+    labelField finalDecomp(points.size(), 0);
+    labelList pointMap(identity(points.size()));
+
+    decompose
+    (
+        globalPointPoints,
+        points,
+        pointWeights,
+        pointMap,       // map back to original points
+        0,
+
+        finalDecomp
+    );
+
+    return finalDecomp;
+}
+
+
+// ************************************************************************* //
diff --git a/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.H b/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.H
new file mode 100644
index 0000000000..a45c01f28e
--- /dev/null
+++ b/src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.H
@@ -0,0 +1,136 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::multiLevelDecomp
+
+Description
+    Decomposition given using consecutive application of decomposers.
+
+SourceFiles
+    multiLevelDecomp.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiLevelDecomp_H
+#define multiLevelDecomp_H
+
+#include "decompositionMethod.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class multiLevelDecomp Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiLevelDecomp
+:
+    public decompositionMethod
+{
+    // Private data
+
+        PtrList<decompositionMethod> methods_;
+
+
+    // Private Member Functions
+
+        //- Given connectivity across processors work out connectivity
+        //  for a (consistent) subset
+        void subsetGlobalCellCells
+        (
+            const labelListList& cellCells,
+            const labelList& set,
+            labelListList& subCellCells,
+            label& cutConnections
+        ) const;
+
+        //- Decompose level methodI without addressing
+        void decompose
+        (
+            const labelListList& pointPoints,
+            const pointField& points,
+            const scalarField& pointWeights,
+            const labelList& pointMap,  // map back to original points
+            const label levelI,
+
+            labelField& finalDecomp
+        );
+
+        //- Disallow default bitwise copy construct and assignment
+        void operator=(const multiLevelDecomp&);
+        multiLevelDecomp(const multiLevelDecomp&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiLevel");
+
+
+    // Constructors
+
+        //- Construct given the decomposition dictionary
+        multiLevelDecomp(const dictionary& decompositionDict);
+
+
+    //- Destructor
+    virtual ~multiLevelDecomp()
+    {}
+
+
+    // Member Functions
+
+        //- Is method parallel aware (i.e. does it synchronize domains across
+        //  proc boundaries)
+        virtual bool parallelAware() const;
+
+        //- Return for every coordinate the wanted processor number. Use the
+        //  mesh connectivity (if needed)
+        virtual labelList decompose
+        (
+            const polyMesh& mesh,
+            const pointField& points,
+            const scalarField& pointWeights
+        );
+
+        //- 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
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C b/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
index 6447c4fa74..f4963e2e5a 100644
--- a/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
@@ -39,13 +39,6 @@ namespace Foam
         simpleGeomDecomp,
         dictionary
     );
-
-    addToRunTimeSelectionTable
-    (
-        decompositionMethod,
-        simpleGeomDecomp,
-        dictionaryMesh
-    );
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -133,26 +126,6 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::simpleGeomDecomp::simpleGeomDecomp(const dictionary& decompositionDict)
-:
-    geomDecomp(decompositionDict, typeName)
-{}
-
-
-Foam::simpleGeomDecomp::simpleGeomDecomp
-(
-    const dictionary& decompositionDict,
-    const polyMesh&
-)
-:
-    geomDecomp(decompositionDict, typeName)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
 Foam::labelList Foam::simpleGeomDecomp::decompose(const pointField& points)
 {
     // construct a list for the final result
@@ -315,4 +288,16 @@ Foam::labelList Foam::simpleGeomDecomp::decompose
 
     return finalDecomp;
 }
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::simpleGeomDecomp::simpleGeomDecomp(const dictionary& decompositionDict)
+:
+    geomDecomp(decompositionDict, typeName)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 // ************************************************************************* //
diff --git a/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H b/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
index 714916f62d..7336b3b85d 100644
--- a/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
@@ -76,13 +76,6 @@ public:
         //- Construct given the decomposition dictionary
         simpleGeomDecomp(const dictionary& decompositionDict);
 
-        //- Construct given the decomposition dictionary and mesh
-        simpleGeomDecomp
-        (
-            const dictionary& decompositionDict,
-            const polyMesh& mesh
-        );
-
 
     //- Destructor
     virtual ~simpleGeomDecomp()
@@ -98,16 +91,24 @@ public:
             return false;
         }
 
-        virtual labelList decompose
-        (
-            const pointField& points
-        );
+        virtual labelList decompose(const pointField&);
+
+        virtual labelList decompose(const pointField&, const scalarField&);
+
+        virtual labelList decompose(const polyMesh&, const pointField& points)
+        {
+            return decompose(points);
+        }
 
         virtual labelList decompose
         (
+            const polyMesh&,
             const pointField& points,
             const scalarField& pointWeights
-        );
+        )
+        {
+            return decompose(points, pointWeights);
+        }
 
         //- Explicitly provided connectivity
         virtual labelList decompose
diff --git a/src/parallel/decompose/metisDecomp/metisDecomp.C b/src/parallel/decompose/metisDecomp/metisDecomp.C
index b198c88415..0c24fc6ac7 100644
--- a/src/parallel/decompose/metisDecomp/metisDecomp.C
+++ b/src/parallel/decompose/metisDecomp/metisDecomp.C
@@ -41,12 +41,7 @@ namespace Foam
 {
     defineTypeNameAndDebug(metisDecomp, 0);
 
-    addToRunTimeSelectionTable
-    (
-        decompositionMethod,
-        metisDecomp,
-        dictionaryMesh
-    );
+    addToRunTimeSelectionTable(decompositionMethod, metisDecomp, dictionary);
 }
 
 
@@ -170,31 +165,31 @@ Foam::label Foam::metisDecomp::decompose
             }
         }
 
-        if (metisCoeffs.readIfPresent("cellWeightsFile", weightsFile))
-        {
-            Info<< "metisDecomp : Using cell-based weights." << endl;
-
-            IOList<int> cellIOWeights
-            (
-                IOobject
-                (
-                    weightsFile,
-                    mesh_.time().timeName(),
-                    mesh_,
-                    IOobject::MUST_READ,
-                    IOobject::AUTO_WRITE
-                )
-            );
-            cellWeights.transfer(cellIOWeights);
-
-            if (cellWeights.size() != xadj.size()-1)
-            {
-                FatalErrorIn("metisDecomp::decompose(const pointField&)")
-                    << "Number of cell weights " << cellWeights.size()
-                    << " does not equal number of cells " << xadj.size()-1
-                    << exit(FatalError);
-            }
-        }
+        //if (metisCoeffs.readIfPresent("cellWeightsFile", weightsFile))
+        //{
+        //    Info<< "metisDecomp : Using cell-based weights." << endl;
+        //
+        //    IOList<int> cellIOWeights
+        //    (
+        //        IOobject
+        //        (
+        //            weightsFile,
+        //            mesh_.time().timeName(),
+        //            mesh_,
+        //            IOobject::MUST_READ,
+        //            IOobject::AUTO_WRITE
+        //        )
+        //    );
+        //    cellWeights.transfer(cellIOWeights);
+        //
+        //    if (cellWeights.size() != xadj.size()-1)
+        //    {
+        //        FatalErrorIn("metisDecomp::decompose(const pointField&)")
+        //            << "Number of cell weights " << cellWeights.size()
+        //            << " does not equal number of cells " << xadj.size()-1
+        //            << exit(FatalError);
+        //    }
+        //}
     }
 
     int nProcs = nProcessors_;
@@ -304,14 +299,9 @@ Foam::label Foam::metisDecomp::decompose
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::metisDecomp::metisDecomp
-(
-    const dictionary& decompositionDict,
-    const polyMesh& mesh
-)
+Foam::metisDecomp::metisDecomp(const dictionary& decompositionDict)
 :
-    decompositionMethod(decompositionDict),
-    mesh_(mesh)
+    decompositionMethod(decompositionDict)
 {}
 
 
@@ -319,11 +309,12 @@ Foam::metisDecomp::metisDecomp
 
 Foam::labelList Foam::metisDecomp::decompose
 (
+    const polyMesh& mesh,
     const pointField& points,
     const scalarField& pointWeights
 )
 {
-    if (points.size() != mesh_.nCells())
+    if (points.size() != mesh.nCells())
     {
         FatalErrorIn
         (
@@ -332,71 +323,53 @@ Foam::labelList Foam::metisDecomp::decompose
             << endl
             << "and supply one coordinate (cellCentre) for every cell." << endl
             << "The number of coordinates " << points.size() << endl
-            << "The number of cells in the mesh " << mesh_.nCells()
+            << "The number of cells in the mesh " << mesh.nCells()
             << exit(FatalError);
     }
 
-    List<int> adjncy;
-    List<int> xadj;
-    calcCSR(mesh_, adjncy, xadj);
+    CompactListList<label> cellCells;
+    calcCellCells(mesh, identity(mesh.nCells()), mesh.nCells(), cellCells);
 
     // Decompose using default weights
-    List<int> finalDecomp;
-    decompose(adjncy, xadj, pointWeights, finalDecomp);
+    labelList decomp;
+    decompose(cellCells.m(), cellCells.offsets(), pointWeights, decomp);
 
-    // Copy back to labelList
-    labelList decomp(finalDecomp.size());
-    forAll(decomp, i)
-    {
-        decomp[i] = finalDecomp[i];
-    }
     return decomp;
 }
 
 
 Foam::labelList Foam::metisDecomp::decompose
 (
+    const polyMesh& mesh,
     const labelList& agglom,
     const pointField& agglomPoints,
     const scalarField& agglomWeights
 )
 {
-    if (agglom.size() != mesh_.nCells())
+    if (agglom.size() != mesh.nCells())
     {
         FatalErrorIn
         (
             "metisDecomp::decompose"
             "(const labelList&, const pointField&, const scalarField&)"
         )   << "Size of cell-to-coarse map " << agglom.size()
-            << " differs from number of cells in mesh " << mesh_.nCells()
+            << " differs from number of cells in mesh " << 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
-    List<int> adjncy;
-    List<int> xadj;
-    {
-        // Get cellCells on coarse mesh.
-        labelListList cellCells;
-        calcCellCells
-        (
-            mesh_,
-            agglom,
-            agglomPoints.size(),
-            cellCells
-        );
 
-        calcCSR(cellCells, adjncy, xadj);
-    }
+    CompactListList<label> cellCells;
+    calcCellCells(mesh, agglom, agglomPoints.size(), cellCells);
 
     // Decompose using default weights
-    List<int> finalDecomp;
-    decompose(adjncy, xadj, agglomWeights, finalDecomp);
+    labelList finalDecomp;
+    decompose(cellCells.m(), cellCells.offsets(), agglomWeights, finalDecomp);
 
 
-    // Rework back into decomposition for original mesh_
+    // Rework back into decomposition for original mesh
     labelList fineDistribution(agglom.size());
 
     forAll(fineDistribution, i)
@@ -404,7 +377,7 @@ Foam::labelList Foam::metisDecomp::decompose
         fineDistribution[i] = finalDecomp[agglom[i]];
     }
 
-    return fineDistribution;
+    return finalDecomp;
 }
 
 
@@ -431,21 +404,12 @@ Foam::labelList Foam::metisDecomp::decompose
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
 
-    List<int> adjncy;
-    List<int> xadj;
-    calcCSR(globalCellCells, adjncy, xadj);
-
+    CompactListList<label> cellCells(globalCellCells);
 
     // Decompose using default weights
-    List<int> finalDecomp;
-    decompose(adjncy, xadj, cellWeights, finalDecomp);
+    labelList decomp;
+    decompose(cellCells.m(), cellCells.offsets(), cellWeights, decomp);
 
-    // Copy back to labelList
-    labelList decomp(finalDecomp.size());
-    forAll(decomp, i)
-    {
-        decomp[i] = finalDecomp[i];
-    }
     return decomp;
 }
 
diff --git a/src/parallel/decompose/metisDecomp/metisDecomp.H b/src/parallel/decompose/metisDecomp/metisDecomp.H
index f18ba8bd0f..d6c7e26e41 100644
--- a/src/parallel/decompose/metisDecomp/metisDecomp.H
+++ b/src/parallel/decompose/metisDecomp/metisDecomp.H
@@ -48,10 +48,6 @@ class metisDecomp
 :
     public decompositionMethod
 {
-    // Private data
-
-        const polyMesh& mesh_;
-
 
     // Private Member Functions
 
@@ -76,12 +72,8 @@ public:
 
     // Constructors
 
-        //- Construct given the decomposition dictionary and mesh
-        metisDecomp
-        (
-            const dictionary& decompositionDict,
-            const polyMesh& mesh
-        );
+        //- Construct given the decomposition dictionary
+        metisDecomp(const dictionary&);
 
 
     //- Destructor
@@ -104,6 +96,7 @@ public:
         //  value. The overall sum of weights might otherwise overflow.
         virtual labelList decompose
         (
+            const polyMesh& mesh,
             const pointField& points,
             const scalarField& pointWeights
         );
@@ -115,6 +108,7 @@ public:
         //  See note on weights above.
         virtual labelList decompose
         (
+            const polyMesh& mesh,
             const labelList& agglom,
             const pointField& regionPoints,
             const scalarField& regionWeights
diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
index cded55bcd7..b4aae52f38 100644
--- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
+++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
@@ -141,12 +141,7 @@ namespace Foam
 {
     defineTypeNameAndDebug(ptscotchDecomp, 0);
 
-    addToRunTimeSelectionTable
-    (
-        decompositionMethod,
-        ptscotchDecomp,
-        dictionaryMesh
-    );
+    addToRunTimeSelectionTable(decompositionMethod, ptscotchDecomp, dictionary);
 }
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -414,14 +409,9 @@ Foam::label Foam::ptscotchDecomp::decompose
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::ptscotchDecomp::ptscotchDecomp
-(
-    const dictionary& decompositionDict,
-    const polyMesh& mesh
-)
+Foam::ptscotchDecomp::ptscotchDecomp(const dictionary& decompositionDict)
 :
-    decompositionMethod(decompositionDict),
-    mesh_(mesh)
+    decompositionMethod(decompositionDict)
 {}
 
 
@@ -429,11 +419,12 @@ Foam::ptscotchDecomp::ptscotchDecomp
 
 Foam::labelList Foam::ptscotchDecomp::decompose
 (
+    const polyMesh& mesh,
     const pointField& points,
     const scalarField& pointWeights
 )
 {
-    if (points.size() != mesh_.nCells())
+    if (points.size() != mesh.nCells())
     {
         FatalErrorIn
         (
@@ -443,7 +434,7 @@ Foam::labelList Foam::ptscotchDecomp::decompose
             << endl
             << "and supply one coordinate (cellCentre) for every cell." << endl
             << "The number of coordinates " << points.size() << endl
-            << "The number of cells in the mesh " << mesh_.nCells()
+            << "The number of cells in the mesh " << mesh.nCells()
             << exit(FatalError);
     }
 
@@ -457,20 +448,14 @@ Foam::labelList Foam::ptscotchDecomp::decompose
     // Make Metis CSR (Compressed Storage Format) storage
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
-    // Connections
-    Field<int> adjncy;
-    // Offsets into adjncy
-    Field<int> xadj;
-    calcDistributedCSR
-    (
-        mesh_,
-        adjncy,
-        xadj
-    );
+
+
+    CompactListList<label> cellCells;
+    calcCellCells(mesh, identity(mesh.nCells()), mesh.nCells(), cellCells);
 
     // Decompose using default weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, pointWeights, finalDecomp);
+    decompose(cellCells.m(), cellCells.offsets(), pointWeights, finalDecomp);
 
     // Copy back to labelList
     labelList decomp(finalDecomp.size());
@@ -484,52 +469,40 @@ Foam::labelList Foam::ptscotchDecomp::decompose
 
 Foam::labelList Foam::ptscotchDecomp::decompose
 (
+    const polyMesh& mesh,
     const labelList& agglom,
     const pointField& agglomPoints,
     const scalarField& pointWeights
 )
 {
-    if (agglom.size() != mesh_.nCells())
+    if (agglom.size() != mesh.nCells())
     {
         FatalErrorIn
         (
             "ptscotchDecomp::decompose(const labelList&, const pointField&)"
         )   << "Size of cell-to-coarse map " << agglom.size()
-            << " differs from number of cells in mesh " << mesh_.nCells()
+            << " differs from number of cells in mesh " << mesh.nCells()
             << exit(FatalError);
     }
 
 //    // For running sequential ...
 //    if (Pstream::nProcs() <= 1)
 //    {
-//        return scotchDecomp(decompositionDict_, mesh_)
+//        return scotchDecomp(decompositionDict_, mesh)
 //            .decompose(agglom, agglomPoints, pointWeights);
 //    }
 
     // Make Metis CSR (Compressed Storage Format) storage
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
-    List<int> adjncy;
-    List<int> xadj;
-    {
-        // Get cellCells on coarse mesh.
-        labelListList cellCells;
-        calcCellCells
-        (
-            mesh_,
-            agglom,
-            agglomPoints.size(),
-            cellCells
-        );
-
-        calcCSR(cellCells, adjncy, xadj);
-    }
+    CompactListList<label> cellCells;
+    calcCellCells(mesh, agglom, agglomPoints.size(), cellCells);
 
     // Decompose using weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, pointWeights, finalDecomp);
+    decompose(cellCells.m(), cellCells.offsets(), pointWeights, finalDecomp);
 
-    // Rework back into decomposition for original mesh_
+    // Rework back into decomposition for original mesh
     labelList fineDistribution(agglom.size());
 
     forAll(fineDistribution, i)
@@ -561,7 +534,7 @@ Foam::labelList Foam::ptscotchDecomp::decompose
 //    // For running sequential ...
 //    if (Pstream::nProcs() <= 1)
 //    {
-//        return scotchDecomp(decompositionDict_, mesh_)
+//        return scotchDecomp(decompositionDict_, mesh)
 //            .decompose(globalCellCells, cellCentres, cWeights);
 //    }
 
@@ -570,13 +543,11 @@ Foam::labelList Foam::ptscotchDecomp::decompose
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
 
-    List<int> adjncy;
-    List<int> xadj;
-    calcCSR(globalCellCells, adjncy, xadj);
+    CompactListList<label> cellCells(globalCellCells);
 
     // Decompose using weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, cWeights, finalDecomp);
+    decompose(cellCells.m(), cellCells.offsets(), cWeights, finalDecomp);
 
     // Copy back to labelList
     labelList decomp(finalDecomp.size());
diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H
index 45649a0b81..d276137d4f 100644
--- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H
+++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H
@@ -48,11 +48,6 @@ class ptscotchDecomp
 :
     public decompositionMethod
 {
-    // Private data
-
-        const polyMesh& mesh_;
-
-
     // Private Member Functions
 
         //- Check and print error message
@@ -80,11 +75,7 @@ public:
     // Constructors
 
         //- Construct given the decomposition dictionary and mesh
-        ptscotchDecomp
-        (
-            const dictionary& decompositionDict,
-            const polyMesh& mesh
-        );
+        ptscotchDecomp(const dictionary& decompositionDict);
 
 
     //- Destructor
@@ -104,6 +95,7 @@ public:
         //  mesh connectivity (if needed)
         virtual labelList decompose
         (
+            const polyMesh& mesh,
             const pointField& points,
             const scalarField& pointWeights
         );
@@ -114,6 +106,7 @@ public:
         //  functionality natively.
         virtual labelList decompose
         (
+            const polyMesh& mesh,
             const labelList& agglom,
             const pointField& regionPoints,
             const scalarField& regionWeights
diff --git a/src/parallel/decompose/scotchDecomp/scotchDecomp.C b/src/parallel/decompose/scotchDecomp/scotchDecomp.C
index 534044646f..dc870abacd 100644
--- a/src/parallel/decompose/scotchDecomp/scotchDecomp.C
+++ b/src/parallel/decompose/scotchDecomp/scotchDecomp.C
@@ -145,7 +145,7 @@ namespace Foam
     (
         decompositionMethod,
         scotchDecomp,
-        dictionaryMesh
+        dictionary
     );
 }
 
@@ -165,6 +165,7 @@ void Foam::scotchDecomp::check(const int retVal, const char* str)
 // Call scotch with options from dictionary.
 Foam::label Foam::scotchDecomp::decompose
 (
+    const fileName& meshPath,
     const List<int>& adjncy,
     const List<int>& xadj,
     const scalarField& cWeights,
@@ -184,7 +185,7 @@ Foam::label Foam::scotchDecomp::decompose
 
             if (writeGraph)
             {
-                OFstream str(mesh_.time().path() / mesh_.name() + ".grf");
+                OFstream str(meshPath + ".grf");
 
                 Info<< "Dumping Scotch graph file to " << str.name() << endl
                     << "Use this in combination with gpart." << endl;
@@ -408,14 +409,9 @@ Foam::label Foam::scotchDecomp::decompose
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::scotchDecomp::scotchDecomp
-(
-    const dictionary& decompositionDict,
-    const polyMesh& mesh
-)
+Foam::scotchDecomp::scotchDecomp(const dictionary& decompositionDict)
 :
-    decompositionMethod(decompositionDict),
-    mesh_(mesh)
+    decompositionMethod(decompositionDict)
 {}
 
 
@@ -423,11 +419,12 @@ Foam::scotchDecomp::scotchDecomp
 
 Foam::labelList Foam::scotchDecomp::decompose
 (
+    const polyMesh& mesh,
     const pointField& points,
     const scalarField& pointWeights
 )
 {
-    if (points.size() != mesh_.nCells())
+    if (points.size() != mesh.nCells())
     {
         FatalErrorIn
         (
@@ -437,20 +434,24 @@ Foam::labelList Foam::scotchDecomp::decompose
             << endl
             << "and supply one coordinate (cellCentre) for every cell." << endl
             << "The number of coordinates " << points.size() << endl
-            << "The number of cells in the mesh " << mesh_.nCells()
+            << "The number of cells in the mesh " << 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
-    List<int> adjncy;
-    List<int> xadj;
-    calcCSR(mesh_, adjncy, xadj);
+
+    CompactListList<label> cellCells;
+    calcCellCells(mesh, identity(mesh.nCells()), mesh.nCells(), cellCells);
 
     // Decompose using default weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, pointWeights, finalDecomp);
+    decompose
+    (
+        mesh.time().path()/mesh.name(),
+        cellCells.m(),
+        cellCells.offsets(),
+        pointWeights,
+        finalDecomp
+    );
 
     // Copy back to labelList
     labelList decomp(finalDecomp.size());
@@ -464,43 +465,35 @@ Foam::labelList Foam::scotchDecomp::decompose
 
 Foam::labelList Foam::scotchDecomp::decompose
 (
+    const polyMesh& mesh,
     const labelList& agglom,
     const pointField& agglomPoints,
     const scalarField& pointWeights
 )
 {
-    if (agglom.size() != mesh_.nCells())
+    if (agglom.size() != mesh.nCells())
     {
         FatalErrorIn
         (
             "scotchDecomp::decompose(const labelList&, const pointField&)"
         )   << "Size of cell-to-coarse map " << agglom.size()
-            << " differs from number of cells in mesh " << mesh_.nCells()
+            << " differs from number of cells in mesh " << 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
-    List<int> adjncy;
-    List<int> xadj;
-    {
-        // Get cellCells on coarse mesh.
-        labelListList cellCells;
-        calcCellCells
-        (
-            mesh_,
-            agglom,
-            agglomPoints.size(),
-            cellCells
-        );
-
-        calcCSR(cellCells, adjncy, xadj);
-    }
+    CompactListList<label> cellCells;
+    calcCellCells(mesh, agglom, agglomPoints.size(), cellCells);
 
     // Decompose using weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, pointWeights, finalDecomp);
+    decompose
+    (
+        mesh.time().path()/mesh.name(),
+        cellCells.m(),
+        cellCells.offsets(),
+        pointWeights,
+        finalDecomp
+    );
 
     // Rework back into decomposition for original mesh_
     labelList fineDistribution(agglom.size());
@@ -537,13 +530,11 @@ Foam::labelList Foam::scotchDecomp::decompose
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
 
-    List<int> adjncy;
-    List<int> xadj;
-    calcCSR(globalCellCells, adjncy, xadj);
+    CompactListList<label> cellCells(globalCellCells);
 
     // Decompose using weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, cWeights, finalDecomp);
+    decompose(".", cellCells.m(), cellCells.offsets(), cWeights, finalDecomp);
 
     // Copy back to labelList
     labelList decomp(finalDecomp.size());
diff --git a/src/parallel/decompose/scotchDecomp/scotchDecomp.H b/src/parallel/decompose/scotchDecomp/scotchDecomp.H
index 249b191f6c..2cd376074c 100644
--- a/src/parallel/decompose/scotchDecomp/scotchDecomp.H
+++ b/src/parallel/decompose/scotchDecomp/scotchDecomp.H
@@ -48,11 +48,6 @@ class scotchDecomp
 :
     public decompositionMethod
 {
-    // Private data
-
-        const polyMesh& mesh_;
-
-
     // Private Member Functions
 
         //- Check and print error message
@@ -60,6 +55,7 @@ class scotchDecomp
 
         label decompose
         (
+            const fileName& meshPath,
             const List<int>& adjncy,
             const List<int>& xadj,
             const scalarField& cWeights,
@@ -80,11 +76,7 @@ public:
     // Constructors
 
         //- Construct given the decomposition dictionary and mesh
-        scotchDecomp
-        (
-            const dictionary& decompositionDict,
-            const polyMesh& mesh
-        );
+        scotchDecomp(const dictionary& decompositionDict);
 
 
     //- Destructor
@@ -107,6 +99,7 @@ public:
         //  might otherwise overflow.
         virtual labelList decompose
         (
+            const polyMesh& mesh,
             const pointField& points,
             const scalarField& pointWeights
         );
@@ -118,6 +111,7 @@ public:
         //  See note on weights above.
         virtual labelList decompose
         (
+            const polyMesh& mesh,
             const labelList& agglom,
             const pointField& regionPoints,
             const scalarField& regionWeights