diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C
index 2cfcc5060f..13078ffe45 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C
@@ -241,12 +241,12 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
     }
 
 
-    // Allocate a communicator for the coarse level
-    label coarseComm = UPstream::allocateCommunicator
-    (
-        fineMesh.comm(),
-        identity(UPstream::nProcs(fineMesh.comm()))  //TBD
-    );
+//    // Allocate a communicator for the coarse level
+//    label coarseComm = UPstream::allocateCommunicator
+//    (
+//        fineMesh.comm(),
+//        identity(UPstream::nProcs(fineMesh.comm()))  //TBD
+//    );
 
 
     // Add the coarse level
@@ -269,7 +269,7 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
                         restrictMap
                     ),
                     fineLevelIndex,
-                    coarseComm
+                    fineMesh.comm() //coarseComm
                 ).ptr()
             );
 
@@ -278,6 +278,18 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
     }
 
 
+    if (debug)
+    {
+        Pout<< "GAMGAgglomeration :"
+            << " agglomerated level " << fineLevelIndex
+            << " from nCells:" << fineMeshAddr.size()
+            << " nFaces:" << upperAddr.size()
+            << " to nCells:" << nCoarseCells
+            << " nFaces:" << coarseOwner.size()
+            << endl;
+    }
+
+
     // Set the coarse ldu addressing onto the list
     meshLevels_.set
     (
@@ -290,11 +302,99 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
             coarseInterfaceAddr,
             coarseInterfaces,
             fineMeshAddr.patchSchedule(),
-            coarseComm,
+            fineMesh.comm(),    //coarseComm,
             true
         )
     );
 }
 
 
+void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
+(
+    const labelList& procAgglomMap,
+    const labelList& procIDs,
+    const label allMeshComm,
+
+    const label levelIndex
+) const
+{
+    if (procCellOffsets_.empty())
+    {
+        Pout<< "GAMGAgglomeration : sizing to " << size()
+            << " levels" << endl;
+        procAgglomMap_.setSize(size());
+        agglomProcIDs_.setSize(size());
+        procMeshLevels_.setSize(size());
+        procCommunicator_.setSize(size(), -1);
+        procCellOffsets_.setSize(size());
+        procFaceMap_.setSize(size());
+        procBoundaryMap_.setSize(size());
+        procBoundaryFaceMap_.setSize(size());
+    }
+
+
+    const lduMesh& myMesh = meshLevels_[levelIndex];
+
+    label meshComm = myMesh.comm();
+
+    label oldWarn = UPstream::warnComm;
+    UPstream::warnComm = meshComm;
+
+    Pout<< "GAMGAgglomeration : gathering myMesh (level="
+        << levelIndex
+        << ") using communicator " << meshComm << endl;
+
+    PtrList<lduMesh> otherMeshes;
+    lduPrimitiveMesh::gather(myMesh, procIDs, otherMeshes);
+
+    Pout<< "** Own Mesh " << myMesh.info() << endl;
+    forAll(otherMeshes, i)
+    {
+        Pout<< "** otherMesh " << i << " "
+            << otherMeshes[i].info()
+            << endl;
+    }
+    Pout<< endl;
+
+    procAgglomMap_.set(levelIndex, new labelList(procAgglomMap));
+    agglomProcIDs_.set(levelIndex, new labelList(procIDs));
+    procCommunicator_[levelIndex] = allMeshComm;
+
+    if (Pstream::myProcNo(meshComm) == procIDs[0])
+    {
+        // Agglomerate all addressing
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        procCellOffsets_.set(levelIndex, new labelList(0));
+        procFaceMap_.set(levelIndex, new labelListList(0));
+        procBoundaryMap_.set(levelIndex, new labelListList(0));
+        procBoundaryFaceMap_.set(levelIndex, new labelListListList(0));
+
+        procMeshLevels_.set
+        (
+            levelIndex,
+            new lduPrimitiveMesh
+            (
+                procCommunicator_[levelIndex],
+                procAgglomMap,
+
+                procIDs,
+                myMesh,
+                otherMeshes,
+
+                procCellOffsets_[levelIndex],
+                procFaceMap_[levelIndex],
+                procBoundaryMap_[levelIndex],
+                procBoundaryFaceMap_[levelIndex]
+            )
+        );
+
+        Pout<< "** Agglomerated Mesh " << procMeshLevels_[levelIndex].info()
+            << endl;
+    }
+
+    UPstream::warnComm = oldWarn;
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
index 5ed00450ab..313018c9c9 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
@@ -197,10 +197,17 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
 Foam::GAMGAgglomeration::~GAMGAgglomeration()
 {
     // Temporary store the user-defined communicators so we can delete them
-    labelList communicators(meshLevels_.size());
+    labelHashSet communicators(meshLevels_.size()+procMeshLevels_.size());
     forAll(meshLevels_, leveli)
     {
-        communicators[leveli] = meshLevels_[leveli].comm();
+        communicators.insert(meshLevels_[leveli].comm());
+    }
+    forAll(procMeshLevels_, leveli)
+    {
+        if (procMeshLevels_.set(leveli))
+        {
+            communicators.insert(procMeshLevels_[leveli].comm());
+        }
     }
 
     Pout<< "~GAMGAgglomeration() : current communicators:" << communicators
@@ -221,9 +228,9 @@ Foam::GAMGAgglomeration::~GAMGAgglomeration()
         }
     }
 
-    forAll(communicators, i)
+    forAllConstIter(labelHashSet, communicators, iter)
     {
-        UPstream::freeCommunicator(communicators[i]);
+        UPstream::freeCommunicator(iter.key());
     }
 }
 
@@ -255,112 +262,76 @@ const Foam::lduInterfacePtrsList& Foam::GAMGAgglomeration::interfaceLevel
 }
 
 
-void Foam::GAMGAgglomeration::gatherMeshes
+const Foam::labelList& Foam::GAMGAgglomeration::procAgglomMap
 (
-    const labelList& procAgglomMap,
-    const labelList& procIDs,
-    const label allMeshComm
+    const label leveli
 ) const
 {
-    if (allMeshPtr_.valid())
-    {
-        FatalErrorIn("GAMGAgglomeration::gatherMeshes(..)")
-            << "Processor-agglomerated mesh already constructed"
-            << exit(FatalError);
-    }
-
-    const lduMesh& coarsestMesh = meshLevels_.last();
-
-    label coarseComm = coarsestMesh.comm();
-
-    label oldWarn = UPstream::warnComm;
-    UPstream::warnComm = coarseComm;
-
-    Pout<< "GAMGAgglomeration : gathering coarsestmesh (level="
-        << meshLevels_.size()-1
-        << ") using communicator " << coarseComm << endl;
-
-    lduPrimitiveMesh::gather(coarsestMesh, procIDs, otherMeshes_);
-
-    Pout<< "** Own Mesh " << coarsestMesh.info() << endl;
-    forAll(otherMeshes_, i)
-    {
-        Pout<< "** otherMesh " << i << " "
-            << otherMeshes_[i].info()
-            << endl;
-    }
-    Pout<< endl;
-
-    if (Pstream::myProcNo(coarseComm) == procIDs[0])
-    {
-        // Agglomerate all addressing
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        allMeshPtr_.reset
-        (
-            new lduPrimitiveMesh
-            (
-                allMeshComm,
-                procAgglomMap,
-
-                procIDs,
-                coarsestMesh,
-                otherMeshes_,
-
-                cellOffsets_,
-                faceMap_,
-                boundaryMap_,
-                boundaryFaceMap_
-            )
-        );
-
-        Pout<< "** Agglomerated Mesh " << allMeshPtr_().info() << endl;
-    }
-
-    UPstream::warnComm = oldWarn;
+    return procAgglomMap_[leveli];
 }
 
 
-const Foam::lduPrimitiveMesh& Foam::GAMGAgglomeration::allMesh() const
+const Foam::labelList& Foam::GAMGAgglomeration::agglomProcIDs
+(
+    const label leveli
+) const
 {
-    if (!allMeshPtr_.valid())
-    {
-        FatalErrorIn("GAMGAgglomeration::allMesh() const")
-            << "Processor-agglomerated mesh not constructed. Did you call"
-            << " gatherMeshes?"
-            << exit(FatalError);
-    }
-    return allMeshPtr_();
+    return agglomProcIDs_[leveli];
 }
 
 
-const Foam::labelList& Foam::GAMGAgglomeration::cellOffsets() const
+bool Foam::GAMGAgglomeration::hasProcMesh(const label leveli) const
 {
-    return cellOffsets_;
+    return procMeshLevels_.set(leveli);
 }
 
 
-const Foam::labelListList& Foam::GAMGAgglomeration::faceMap() const
+const Foam::lduMesh& Foam::GAMGAgglomeration::procMeshLevel(const label leveli)
+const
 {
-    return faceMap_;
+    return procMeshLevels_[leveli];
 }
 
 
-const Foam::labelListList& Foam::GAMGAgglomeration::boundaryMap() const
+Foam::label Foam::GAMGAgglomeration::procCommunicator(const label leveli) const
 {
-    return boundaryMap_;
+    return procCommunicator_[leveli];
 }
 
 
-const Foam::labelListListList& Foam::GAMGAgglomeration::boundaryFaceMap() const
+const Foam::labelList& Foam::GAMGAgglomeration::cellOffsets
+(
+    const label leveli
+) const
 {
-    return boundaryFaceMap_;
+    return procCellOffsets_[leveli];
 }
 
 
-const Foam::PtrList<Foam::lduMesh>& Foam::GAMGAgglomeration::otherMeshes() const
+const Foam::labelListList& Foam::GAMGAgglomeration::faceMap
+(
+    const label leveli
+) const
 {
-    return otherMeshes_;
+    return procFaceMap_[leveli];
+}
+
+
+const Foam::labelListList& Foam::GAMGAgglomeration::boundaryMap
+(
+    const label leveli
+) const
+{
+    return procBoundaryMap_[leveli];
+}
+
+
+const Foam::labelListListList& Foam::GAMGAgglomeration::boundaryFaceMap
+(
+    const label leveli
+) const
+{
+    return procBoundaryFaceMap_[leveli];
 }
 
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H
index 7d218956ea..fd281965dc 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H
@@ -30,7 +30,6 @@ Description
 SourceFiles
     GAMGAgglomeration.C
     GAMGAgglomerationTemplates.C
-    GAMGAgglomerate.C
     GAMGAgglomerateLduAddressing.C
 
 \*---------------------------------------------------------------------------*/
@@ -102,25 +101,30 @@ protected:
 
         //- Processor agglomeration
 
-            //- Combined coarsest meshes for all processors
-            mutable autoPtr<lduPrimitiveMesh> allMeshPtr_;
+            //- Per level, per processor the processor it agglomerates into
+            mutable PtrList<labelList> procAgglomMap_;
 
-            //- Per remote processor the local equivalent. TBD: hack. Can be
-            //  deleted once lduMatrices get received and merged into single
-            //  one without constructing local ones
-            mutable PtrList<lduMesh> otherMeshes_;
+            //- Per level the set of processors to agglomerate. Element 0 is
+            //  the 'master' of the cluster.
+            mutable PtrList<labelList> agglomProcIDs_;
 
-            //- Mapping from processor to coarsestAllMesh cells
-            mutable labelList cellOffsets_;
+            //- Combined mesh for all processors
+            mutable PtrList<lduPrimitiveMesh> procMeshLevels_;
 
-            //- Mapping from processor to coarsestAllMesh face
-            mutable labelListList faceMap_;
+            //- Communicator for given level
+            mutable labelList procCommunicator_;
 
-            //- Mapping from processor to coarsestAllMesh boundary
-            mutable labelListList boundaryMap_;
+            //- Mapping from processor to procMeshLevel cells
+            mutable PtrList<labelList> procCellOffsets_;
 
-            //- Mapping from processor to coarsestAllMesh boundary face
-            mutable labelListListList boundaryFaceMap_;
+            //- Mapping from processor to procMeshLevel face
+            mutable PtrList<labelListList> procFaceMap_;
+
+            //- Mapping from processor to procMeshLevel boundary
+            mutable PtrList<labelListList> procBoundaryMap_;
+
+            //- Mapping from processor to procMeshLevel boundary face
+            mutable PtrList<labelListListList> procBoundaryFaceMap_;
 
 
         //- Assemble coarse mesh addressing
@@ -284,32 +288,45 @@ public:
             //                    (rank in allMeshComm!). Global information.
             //  - procIDs       : local information: same for all in
             //                    agglomerated processor.
-            void gatherMeshes
+            void procAgglomerateLduAddressing
             (
                 const labelList& procAgglomMap,
                 const labelList& procIDs,
-                const label allMeshComm
+                const label allMeshComm,
+                const label levelIndex
             ) const;
 
+
+            //- Mapping from processor to agglomerated processor (global, all
+            //  processors have the same information)
+            const labelList& procAgglomMap(const label leveli) const;
+
+            //- Set of processors to agglomerate. Element 0 is the
+            //  master processor. (local, same only on those processor that
+            //  agglomerate)
+            const labelList& agglomProcIDs(const label leveli) const;
+
+            //- Check that level has combined mesh
+            bool hasProcMesh(const label leveli) const;
+
             //- Combined mesh
-            const lduPrimitiveMesh& allMesh() const;
+            const lduMesh& procMeshLevel(const label leveli) const;
 
-            //- Mapping from processor to coarsestAllMesh cells
-            const labelList& cellOffsets() const;
+            //- Communicator for procMesh (stored separately so also works
+            //  even if no mesh is stored on this processor)
+            label procCommunicator(const label leveli) const;
 
-            //- Mapping from processor to coarsestAllMesh face
-            const labelListList& faceMap() const;
+            //- Mapping from processor to procMesh cells
+            const labelList& cellOffsets(const label leveli) const;
 
-            //- Mapping from processor to coarsestAllMesh boundary
-            const labelListList& boundaryMap() const;
+            //- Mapping from processor to procMesh face
+            const labelListList& faceMap(const label leveli) const;
 
-            //- Mapping from processor to coarsestAllMesh boundary face
-            const labelListListList& boundaryFaceMap() const;
+            //- Mapping from processor to procMesh boundary
+            const labelListList& boundaryMap(const label leveli) const;
 
-            //- Per remote processor the local equivalent. TBD: hack. Can be
-            //  deleted once lduMatrices get received and merged into single
-            //  one without constructing local ones
-            const PtrList<lduMesh>& otherMeshes() const;
+            //- Mapping from processor to procMesh boundary face
+            const labelListListList& boundaryFaceMap(const label leveli) const;
 
 };
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C
index 2568ad56fe..3662f1c354 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C
@@ -91,10 +91,13 @@ Foam::GAMGSolver::GAMGSolver
 
     if (matrixLevels_.size())
     {
-        const label coarsestLevel = matrixLevels_.size() - 1;
-
         if (directSolveCoarsest_)
         {
+            const label coarsestLevel = matrixLevels_.size() - 1;
+
+            Pout<< "GAMGSolver :"
+                << " coarsestLevel:" << coarsestLevel << endl;
+
             const lduMesh& coarsestMesh = matrixLevels_[coarsestLevel].mesh();
 
             label coarseComm = coarsestMesh.comm();
@@ -118,73 +121,63 @@ Foam::GAMGSolver::GAMGSolver
         }
         else if (processorAgglomerate_)
         {
-            const lduMatrix& coarsestMatrix = matrixLevels_[coarsestLevel];
-            const lduInterfaceFieldPtrsList& coarsestInterfaces =
-                interfaceLevels_[coarsestLevel];
-            const FieldField<Field, scalar>& coarsestBouCoeffs =
-                interfaceLevelsBouCoeffs_[coarsestLevel];
-            const FieldField<Field, scalar>& coarsestIntCoeffs =
-                interfaceLevelsIntCoeffs_[coarsestLevel];
-            const lduMesh& coarsestMesh = coarsestMatrix.mesh();
+            // Pick a level to processor agglomerate
+            label agglomLevel = matrixLevels_.size() - 1;//1;
 
 
-            label coarseComm = coarsestMesh.comm();
+            // Get mesh and matrix at this level
+            const lduMatrix& levelMatrix = matrixLevels_[agglomLevel];
+            const lduMesh& levelMesh = levelMatrix.mesh();
+
+
+            label levelComm = levelMesh.comm();
             label oldWarn = UPstream::warnComm;
-            UPstream::warnComm = coarseComm;
+            UPstream::warnComm = levelComm;
 
-            Pout<< "Solve generic on coarsestmesh (level=" << coarsestLevel
-                << ") using communicator " << coarseComm << endl;
-
-            //const List<int>& procIDs = UPstream::procID(coarseComm);
+            Pout<< "Solve generic on mesh (level=" << agglomLevel
+                << ") using communicator " << levelComm << endl;
 
             // Processor restriction map: per processor the coarse processor
-            labelList procAgglomMap(UPstream::nProcs(coarseComm));
-            procAgglomMap[0] = 0;
-            procAgglomMap[1] = 0;
-            procAgglomMap[2] = 1;
-            procAgglomMap[3] = 1;
-
-            // Determine the master processors
-            Map<label> agglomToMaster(procAgglomMap.size());
-
-            forAll(procAgglomMap, procI)
-            {
-                label coarseI = procAgglomMap[procI];
-
-                Map<label>::iterator fnd = agglomToMaster.find(coarseI);
-                if (fnd == agglomToMaster.end())
-                {
-                    agglomToMaster.insert(coarseI, procI);
-                }
-                else
-                {
-                    fnd() = max(fnd(), procI);
-                }
-            }
-
-            labelList masterProcs(agglomToMaster.size());
-            forAllConstIter(Map<label>, agglomToMaster, iter)
-            {
-                masterProcs[iter.key()] = iter();
-            }
-
-            // Allocate a communicator for the linear solver
-
-            label masterComm = UPstream::allocateCommunicator
-            (
-                coarseComm,
-                masterProcs
-            );
-            Pout<< "** Allocated communicator " << masterComm
-                << " for indices " << masterProcs
-                << " in processor list " << UPstream::procID(coarseComm)
-                << endl;
-
-
+            labelList procAgglomMap(UPstream::nProcs(levelComm));
+            // Master processor
+            labelList masterProcs;
+            // Local processors that agglomerate. agglomProcIDs[0] is in
+            // masterProc.
             List<int> agglomProcIDs;
-            // Collect all the processors in my agglomeration
+
             {
-                label myProcID = Pstream::myProcNo(coarseComm);
+                procAgglomMap[0] = 0;
+                procAgglomMap[1] = 0;
+                procAgglomMap[2] = 1;
+                procAgglomMap[3] = 1;
+
+                // Determine the master processors
+                Map<label> agglomToMaster(procAgglomMap.size());
+
+                forAll(procAgglomMap, procI)
+                {
+                    label coarseI = procAgglomMap[procI];
+
+                    Map<label>::iterator fnd = agglomToMaster.find(coarseI);
+                    if (fnd == agglomToMaster.end())
+                    {
+                        agglomToMaster.insert(coarseI, procI);
+                    }
+                    else
+                    {
+                        fnd() = max(fnd(), procI);
+                    }
+                }
+
+                masterProcs.setSize(agglomToMaster.size());
+                forAllConstIter(Map<label>, agglomToMaster, iter)
+                {
+                    masterProcs[iter.key()] = iter();
+                }
+
+
+                // Collect all the processors in my agglomeration
+                label myProcID = Pstream::myProcNo(levelComm);
                 label myAgglom = procAgglomMap[myProcID];
 
                 // Get all processors agglomerating to the same coarse processor
@@ -198,26 +191,34 @@ Foam::GAMGSolver::GAMGSolver
                 Swap(agglomProcIDs[0], agglomProcIDs[index]);
             }
 
+
             Pout<< "procAgglomMap:" << procAgglomMap << endl;
             Pout<< "agglomProcIDs:" << agglomProcIDs << endl;
 
+            // Allocate a communicator for the processor-agglomerated matrix
+            label procAgglomComm = UPstream::allocateCommunicator
+            (
+                levelComm,
+                masterProcs
+            );
+            Pout<< "** Allocated communicator " << procAgglomComm
+                << " for indices " << masterProcs
+                << " in processor list " << UPstream::procID(levelComm)
+                << endl;
+
+
 
             // Gather matrix and mesh onto agglomProcIDs[0]
             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
             procAgglomerateMatrix
             (
-                // Current mesh and  matrix
-                coarsestMesh,
-                coarsestMatrix,
-                coarsestInterfaces,
-                coarsestBouCoeffs,
-                coarsestIntCoeffs,
-
                 // Agglomeration information
                 procAgglomMap,
                 agglomProcIDs,
-                masterComm,
+                procAgglomComm,
+
+                agglomLevel,  // level (coarse, not fine level!)
 
                 // Resulting matrix
                 allMatrixPtr_,
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
index 9b13c25459..955a798318 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
@@ -192,18 +192,13 @@ class GAMGSolver
         //- Agglomerate processor matrices
         void procAgglomerateMatrix
         (
-            // Current mesh and  matrix
-            const lduMesh& coarsestMesh,
-            const lduMatrix& coarsestMatrix,
-            const lduInterfaceFieldPtrsList& coarsestInterfaces,
-            const FieldField<Field, scalar>& coarsestBouCoeffs,
-            const FieldField<Field, scalar>& coarsestIntCoeffs,
-
             // Agglomeration information
             const labelList& procAgglomMap,
             const List<int>& agglomProcIDs,
             label masterComm,
 
+            const label levelI,
+
             // Resulting matrix
             autoPtr<lduMatrix>& allMatrixPtr,
             FieldField<Field, scalar>& allInterfaceBouCoeffs,
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C
index 66ede61db7..6bd97630c9 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C
@@ -352,18 +352,13 @@ void Foam::GAMGSolver::gatherMatrices
 
 void Foam::GAMGSolver::procAgglomerateMatrix
 (
-    // Current mesh and  matrix
-    const lduMesh& coarsestMesh,
-    const lduMatrix& coarsestMatrix,
-    const lduInterfaceFieldPtrsList& coarsestInterfaces,
-    const FieldField<Field, scalar>& coarsestBouCoeffs,
-    const FieldField<Field, scalar>& coarsestIntCoeffs,
-
     // Agglomeration information
     const labelList& procAgglomMap,
     const List<int>& agglomProcIDs,
     label masterComm,
 
+    const label levelI,
+
     // Resulting matrix
     autoPtr<lduMatrix>& allMatrixPtr,
     FieldField<Field, scalar>& allInterfaceBouCoeffs,
@@ -372,6 +367,16 @@ void Foam::GAMGSolver::procAgglomerateMatrix
     lduInterfaceFieldPtrsList& allInterfaces
 )
 {
+    const lduMatrix& coarsestMatrix = matrixLevels_[levelI];
+    const lduInterfaceFieldPtrsList& coarsestInterfaces =
+        interfaceLevels_[levelI];
+    const FieldField<Field, scalar>& coarsestBouCoeffs =
+        interfaceLevelsBouCoeffs_[levelI];
+    const FieldField<Field, scalar>& coarsestIntCoeffs =
+        interfaceLevelsIntCoeffs_[levelI];
+    const lduMesh& coarsestMesh = coarsestMatrix.mesh();
+
+
     label coarseComm = coarsestMesh.comm();
 
     label oldWarn = UPstream::warnComm;
@@ -379,7 +384,16 @@ void Foam::GAMGSolver::procAgglomerateMatrix
 
 
     // Construct (on the agglomeration) a complete mesh with mapping
-    agglomeration_.gatherMeshes(procAgglomMap, agglomProcIDs, masterComm);
+    Pout<< "procAgglomerateMatrix :" << " level:" << levelI << endl;
+
+
+    agglomeration_.procAgglomerateLduAddressing
+    (
+        procAgglomMap,
+        agglomProcIDs,
+        masterComm,
+        levelI     // coarsest level
+    );
 
 
     // Gather all matrix coefficients onto agglomProcIDs[0]
@@ -425,7 +439,13 @@ void Foam::GAMGSolver::procAgglomerateMatrix
         Pout<< endl;
 
 
-        const lduMesh& allMesh = agglomeration_.allMesh();
+        const lduMesh& allMesh = agglomeration_.procMeshLevel(levelI);
+        const labelList& cellOffsets = agglomeration_.cellOffsets(levelI);
+        const labelListList& faceMap = agglomeration_.faceMap(levelI);
+        const labelListList& boundaryMap = agglomeration_.boundaryMap(levelI);
+        const labelListListList& boundaryFaceMap =
+            agglomeration_.boundaryFaceMap(levelI);
+
 
         allMatrixPtr.reset(new lduMatrix(allMesh));
         lduMatrix& allMatrix = allMatrixPtr();
@@ -447,7 +467,7 @@ void Foam::GAMGSolver::procAgglomerateMatrix
                 (
                     allDiag,
                     otherMats[i].diag().size(),
-                    agglomeration_.cellOffsets()[i+1]
+                    cellOffsets[i+1]
                 ).assign
                 (
                     otherMats[i].diag()
@@ -460,14 +480,14 @@ void Foam::GAMGSolver::procAgglomerateMatrix
             UIndirectList<scalar>
             (
                 allLower,
-                agglomeration_.faceMap()[0]
+                faceMap[0]
             ) = coarsestMatrix.lower();
             forAll(otherMats, i)
             {
                 UIndirectList<scalar>
                 (
                     allLower,
-                    agglomeration_.faceMap()[i+1]
+                    faceMap[i+1]
                 ) = otherMats[i].lower();
             }
         }
@@ -477,14 +497,14 @@ void Foam::GAMGSolver::procAgglomerateMatrix
             UIndirectList<scalar>
             (
                 allUpper,
-                agglomeration_.faceMap()[0]
+                faceMap[0]
             ) = coarsestMatrix.upper();
             forAll(otherMats, i)
             {
                 UIndirectList<scalar>
                 (
                     allUpper,
-                    agglomeration_.faceMap()[i+1]
+                    faceMap[i+1]
                 ) = otherMats[i].upper();
             }
         }
@@ -510,7 +530,7 @@ void Foam::GAMGSolver::procAgglomerateMatrix
         }
 
         labelList nBounFaces(allMeshInterfaces.size());
-        forAll(agglomeration_.boundaryMap(), procI)
+        forAll(boundaryMap, procI)
         {
             const FieldField<Field, scalar>& procBouCoeffs
             (
@@ -525,7 +545,7 @@ void Foam::GAMGSolver::procAgglomerateMatrix
               : otherIntCoeffs[procI-1]
             );
 
-            const labelList& bMap = agglomeration_.boundaryMap()[procI];
+            const labelList& bMap = boundaryMap[procI];
             forAll(bMap, procIntI)
             {
                 label allIntI = bMap[procIntI];
@@ -587,8 +607,7 @@ void Foam::GAMGSolver::procAgglomerateMatrix
                     scalarField& allBou = allInterfaceBouCoeffs[allIntI];
                     scalarField& allInt = allInterfaceIntCoeffs[allIntI];
 
-                    const labelList& map =
-                        agglomeration_.boundaryFaceMap()[procI][procIntI];
+                    const labelList& map = boundaryFaceMap[procI][procIntI];
 
 Pout<< "    from proc:" << procI << " interface:" << procIntI
 << " mapped to faces:" << map << endl;
@@ -612,12 +631,11 @@ Pout<< "    from proc:" << procI << " interface:" << procIntI
                 {
                     // Boundary has become internal face
 
-                    const labelList& map =
-                        agglomeration_.boundaryFaceMap()[procI][procIntI];
+                    const labelList& map = boundaryFaceMap[procI][procIntI];
 
-                    const labelList& allU = allMesh.lduAddr().upperAddr();
-                    const labelList& allL = allMesh.lduAddr().lowerAddr();
-                    const label off = agglomeration_.cellOffsets()[procI];
+                    //const labelList& allU = allMesh.lduAddr().upperAddr();
+                    //const labelList& allL = allMesh.lduAddr().lowerAddr();
+                    //const label off = cellOffsets[procI];
 
                     const scalarField& procBou = procBouCoeffs[procIntI];
                     const scalarField& procInt = procIntCoeffs[procIntI];
@@ -652,41 +670,41 @@ Pout<< "    from proc:" << procI << " interface:" << procIntI
                         }
 
 
-                        // Simple check
-                        label allFaceI =
-                        (
-                            map[i] >= 0
-                          ? map[i]
-                          : -map[i]-1
-                        );
-
-                        const labelList& fCells =
-                            lduPrimitiveMesh::mesh
-                            (
-                                coarsestMesh,
-                                agglomeration_.otherMeshes(),
-                                procI
-                            ).lduAddr().patchAddr(procIntI);
-
-                        label allCellI = off + fCells[i];
-
-                        if
-                        (
-                            allCellI != allL[allFaceI]
-                         && allCellI != allU[allFaceI]
-                        )
-                        {
-                            FatalErrorIn
-                            (
-                                "GAMGSolver::GAMGSolver()"
-                            )   << "problem."
-                                << " allFaceI:" << allFaceI
-                                << " local cellI:" << fCells[i]
-                                << " allCellI:" << allCellI
-                                << " allLower:" << allL[allFaceI]
-                                << " allUpper:" << allU[allFaceI]
-                                << abort(FatalError);
-                        }
+                        //// Simple check
+                        //label allFaceI =
+                        //(
+                        //    map[i] >= 0
+                        //  ? map[i]
+                        //  : -map[i]-1
+                        //);
+                        //
+                        //const labelList& fCells =
+                        //    lduPrimitiveMesh::mesh
+                        //    (
+                        //        coarsestMesh,
+                        //        agglomeration_.otherMeshes(),
+                        //        procI
+                        //    ).lduAddr().patchAddr(procIntI);
+                        //
+                        //label allCellI = off + fCells[i];
+                        //
+                        //if
+                        //(
+                        //    allCellI != allL[allFaceI]
+                        // && allCellI != allU[allFaceI]
+                        //)
+                        //{
+                        //    FatalErrorIn
+                        //    (
+                        //        "GAMGSolver::GAMGSolver()"
+                        //    )   << "problem."
+                        //        << " allFaceI:" << allFaceI
+                        //        << " local cellI:" << fCells[i]
+                        //        << " allCellI:" << allCellI
+                        //        << " allLower:" << allL[allFaceI]
+                        //        << " allUpper:" << allU[allFaceI]
+                        //        << abort(FatalError);
+                        //}
                     }
                 }
             }
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
index 8b04036c47..f62a4ee769 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
@@ -442,6 +442,11 @@ void Foam::GAMGSolver::solveCoarsestLevel
 ) const
 {
     const label coarsestLevel = matrixLevels_.size() - 1;
+
+    Pout<< "solveCoarsestLevel :"
+        << " coarsestLevel:" << coarsestLevel << endl;
+
+
     label coarseComm = matrixLevels_[coarsestLevel].mesh().comm();
     label oldWarn = UPstream::warnComm;
     UPstream::warnComm = coarseComm;
@@ -453,53 +458,21 @@ void Foam::GAMGSolver::solveCoarsestLevel
     }
     else if (processorAgglomerate_)
     {
+        const labelList& agglomProcIDs = agglomeration_.agglomProcIDs
+        (
+            coarsestLevel
+        );
 
-//XXXX
-        labelList procAgglomMap(UPstream::nProcs(coarseComm));
-        procAgglomMap[0] = 0;
-        procAgglomMap[1] = 0;
-        procAgglomMap[2] = 1;
-        procAgglomMap[3] = 1;
 
-        // Determine the master processors
-        Map<label> agglomToMaster(procAgglomMap.size());
-
-        forAll(procAgglomMap, procI)
-        {
-            label coarseI = procAgglomMap[procI];
-
-            Map<label>::iterator fnd = agglomToMaster.find(coarseI);
-            if (fnd == agglomToMaster.end())
-            {
-                agglomToMaster.insert(coarseI, procI);
-            }
-            else
-            {
-                fnd() = max(fnd(), procI);
-            }
-        }
-
-        List<int> agglomProcIDs;
-        // Collect all the processors in my agglomeration
-        {
-            label myProcID = Pstream::myProcNo(coarseComm);
-            label myAgglom = procAgglomMap[myProcID];
-
-            // Get all processors agglomerating to the same coarse processor
-            agglomProcIDs = findIndices(procAgglomMap, myAgglom);
-            // Make sure the master is the first element.
-            label index = findIndex
-            (
-                agglomProcIDs,
-                agglomToMaster[myAgglom]
-            );
-            Swap(agglomProcIDs[0], agglomProcIDs[index]);
-        }
-//XXXX
 
         scalarField allSource;
 
-        const globalIndex cellOffsets(agglomeration_.cellOffsets());
+        globalIndex cellOffsets;
+        if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
+        {
+            cellOffsets.offsets() =  agglomeration_.cellOffsets(coarsestLevel);
+        }
+
         cellOffsets.gather
         (
             coarseComm,
@@ -511,16 +484,17 @@ void Foam::GAMGSolver::solveCoarsestLevel
         scalarField allCorrField;
         solverPerformance coarseSolverPerf;
 
-        if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
+        label solveComm = agglomeration_.procCommunicator(coarsestLevel);
+        label oldWarn = UPstream::warnComm;
+        UPstream::warnComm = solveComm;
+
+        if (Pstream::myProcNo(solveComm) != -1)
         {
             const lduMatrix& allMatrix = allMatrixPtr_();
 
             allCorrField.setSize(allSource.size(), 0);
 
             {
-                label solveComm = allMatrix.mesh().comm();
-                label oldWarn = UPstream::warnComm;
-                UPstream::warnComm = solveComm;
                 Pout<< "** Master:Solving on comm:" << solveComm
                     << " with procs:" << UPstream::procID(solveComm) << endl;
 
@@ -558,10 +532,10 @@ void Foam::GAMGSolver::solveCoarsestLevel
                         allSource
                     );
                 }
-                UPstream::warnComm = oldWarn;
             }
         }
 
+        UPstream::warnComm = oldWarn;
         Pout<< "done master solve." << endl;
 
         // Scatter to all processors