diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 2559fe7344fcecb3a9505145734a2f4f021a82bf..2c01641321e1085592aaff8666f44afa72ea0205 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -320,6 +320,14 @@ dummyAgglomeration = $(GAMGAgglomerations)/dummyAgglomeration
 $(dummyAgglomeration)/dummyAgglomeration.C
 
 
+GAMGProcAgglomerations = $(GAMG)/GAMGProcAgglomerations
+
+GAMGProcAgglomeration = $(GAMGProcAgglomerations)/GAMGProcAgglomeration
+$(GAMGProcAgglomeration)/GAMGProcAgglomeration.C
+masterCoarsest = $(GAMGProcAgglomerations)/masterCoarsest
+$(masterCoarsest)/masterCoarsest.C
+
+
 meshes/lduMesh/lduMesh.C
 meshes/lduMesh/lduPrimitiveMesh.C
 
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 6d806b7789bc1ae6b75db02b16d43fecfe57f253..63405e58cc0a432eaabff833dda5468163af52ea 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C
@@ -248,12 +248,6 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
         }
     }
 
-    Pout<< "Total faces:" << faceRestrictAddr.size()
-        << " nDissapear:" << nDissapear
-        << " unflipped:" << faceRestrictAddr.size()-nFlipped
-        << " flipped:" << nFlipped
-        << endl;
-
 
 
     // Clear the temporary storage for the coarse cell data
@@ -385,73 +379,6 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
 }
 
 
-void Foam::GAMGAgglomeration::gatherMeshes
-(
-    const label meshComm,
-    const labelList& procAgglomMap,
-    const labelList& procIDs,
-    const label allMeshComm,
-
-    const lduMesh& myMesh,
-
-    autoPtr<lduPrimitiveMesh>& allMesh,
-    labelList& procCellOffsets,
-    labelListList& procFaceMap,
-    labelListList& procBoundaryMap,
-    labelListListList& procBoundaryFaceMap
-)
-{
-    label oldWarn = UPstream::warnComm;
-    UPstream::warnComm = meshComm;
-
-    //Pout<< "GAMGAgglomeration : gathering myMesh (level="
-    //    << levelIndex
-    //    << ") using communicator " << meshComm << endl;
-
-    PtrList<lduMesh> otherMeshes;
-    lduPrimitiveMesh::gather(meshComm, myMesh, procIDs, otherMeshes);
-
-    if (Pstream::myProcNo(meshComm) == procIDs[0])
-    {
-        // Combine all addressing
-        // ~~~~~~~~~~~~~~~~~~~~~~
-
-        //Pout<< "Own Mesh " << myMesh.lduAddr().size() << endl;
-        //forAll(otherMeshes, i)
-        //{
-        //    Pout<< "    otherMesh " << i << " "
-        //        << otherMeshes[i].lduAddr().size()
-        //        << endl;
-        //}
-
-        labelList procFaceOffsets;
-
-        allMesh.reset
-        (
-            new lduPrimitiveMesh
-            (
-                allMeshComm,
-                procAgglomMap,
-
-                procIDs,
-                myMesh,
-                otherMeshes,
-
-                procCellOffsets,
-                procFaceOffsets,
-                procFaceMap,
-                procBoundaryMap,
-                procBoundaryFaceMap
-            )
-        );
-
-        //Pout<< "** Agglomerated Mesh " << allMesh().info();
-    }
-
-    UPstream::warnComm = oldWarn;
-}
-
-
 void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
 (
     const label meshComm,
@@ -462,49 +389,8 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
     const label levelIndex
 )
 {
-    Pout<< "** GAMGAgglomeration:procAgglomerateLduAddressing **" << endl;
-    Pout<< "level:" << levelIndex
-        << " havefinemesh:" << meshLevels_.set(levelIndex-1) << endl;
-
     const lduMesh& myMesh = meshLevels_[levelIndex-1];
 
-    Pout<< "my mesh                     :" << myMesh.info();
-    Pout<< "nCoarseCells                :" << nCells_[levelIndex] << endl;
-    Pout<< "restrictAddressing_         :"
-        << restrictAddressing_[levelIndex].size()
-        << " max:" << max(restrictAddressing_[levelIndex]) << endl;
-    Pout<< "faceRestrictAddressing_     :"
-        << faceRestrictAddressing_[levelIndex].size()
-        << " max:" << max(faceRestrictAddressing_[levelIndex]) << endl;
-    Pout<< "patchFaceRestrictAddressing_:" << endl;
-    forAll(patchFaceRestrictAddressing_[levelIndex], patchI)
-    {
-        const labelList& map =
-            patchFaceRestrictAddressing_[levelIndex][patchI];
-
-        if (map.size())
-        {
-            Pout<< "    patch:" << patchI
-                << " size:" << map.size()
-                << " max:" << max(map)
-                << endl;
-        }
-    }
-    Pout<< "interfaceLevels_:" << endl;
-    forAll(interfaceLevels_[levelIndex], patchI)
-    {
-        if (interfaceLevels_[levelIndex].set(patchI))
-        {
-            Pout<< "    patch:" << patchI
-                << " interface:" << interfaceLevels_[levelIndex][patchI].type()
-                << " size:"
-                << interfaceLevels_[levelIndex][patchI].faceCells().size()
-                << endl;
-        }
-    }
-
-
-
 
     label oldWarn = UPstream::warnComm;
     UPstream::warnComm = meshComm;
@@ -522,21 +408,36 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
     procBoundaryFaceMap_.set(levelIndex, new labelListListList(0));
 
     autoPtr<lduPrimitiveMesh> allMesh;
-    gatherMeshes
-    (
-        meshComm,
-        procAgglomMap,
-        procIDs,
-        allMeshComm,
+    {
+        // Collect meshes
+        PtrList<lduMesh> otherMeshes;
+        lduPrimitiveMesh::gather(meshComm, myMesh, procIDs, otherMeshes);
 
-        myMesh,
+        if (Pstream::myProcNo(meshComm) == procIDs[0])
+        {
+            // Combine all addressing
 
-        allMesh,
-        procCellOffsets_[levelIndex],
-        procFaceMap_[levelIndex],
-        procBoundaryMap_[levelIndex],
-        procBoundaryFaceMap_[levelIndex]
-    );
+            labelList procFaceOffsets;
+            allMesh.reset
+            (
+                new lduPrimitiveMesh
+                (
+                    allMeshComm,
+                    procAgglomMap,
+
+                    procIDs,
+                    myMesh,
+                    otherMeshes,
+
+                    procCellOffsets_[levelIndex],
+                    procFaceOffsets,
+                    procFaceMap_[levelIndex],
+                    procBoundaryMap_[levelIndex],
+                    procBoundaryFaceMap_[levelIndex]
+                )
+            );
+        }
+    }
 
 
     if (Pstream::myProcNo(meshComm) == procIDs[0])
@@ -566,7 +467,6 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
     // Combine restrict addressing
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-
     procAgglomerateRestrictAddressing
     (
         meshComm,
@@ -578,19 +478,6 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
     {
         clearLevel(levelIndex);
     }
-    else
-    {
-        //const lduPrimitiveMesh& levelMesh = meshLevels_[levelIndex-1];
-        //Pout<< "** DONE GAMGAgglomeration:procAgglomerateLduAddressing **"
-        //    << endl;
-        //Pout<< "my mesh:" << levelMesh.info();
-        //Pout<< "nCoarseCells:" << nCells_[levelIndex] << endl;
-        //Pout<< "restrictAddressing_:"
-        //    << restrictAddressing_[levelIndex].size()
-        //    << " max:" << max(restrictAddressing_[levelIndex]) << endl;
-        //Pout<< "** DONE GAMGAgglomeration:procAgglomerateLduAddressing **"
-        //    << endl;
-    }
 
     UPstream::warnComm = oldWarn;
 }
@@ -603,16 +490,6 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
     const label levelIndex
 )
 {
-    //Pout<< "** GAMGAgglomeration:procAgglomerateRestrictAddressing **"
-    //    << endl;
-    //Pout<< "level:" << levelIndex << endl;
-    //Pout<< "procIDs:" << procIDs << endl;
-    //Pout<< "myProcNo:" << UPstream::myProcNo(comm) << endl;
-    //Pout<< "procCellOffsets_:" << procCellOffsets_[levelIndex] << endl;
-
-    //const lduMesh& levelMesh = meshLevel(levelIndex);
-    //Pout<< "fine:" << restrictAddressing_[levelIndex].size() << endl;
-
     // Collect number of cells
     labelList nFineCells;
     gatherList
@@ -622,9 +499,6 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
         restrictAddressing_[levelIndex].size(),
         nFineCells
     );
-    //const labelList& offsets = procCellOffsets_[levelIndex];
-    //Pout<< "offsets:" << offsets << endl;
-    //Pout<< "nFineCells:" << nFineCells << endl;
 
     labelList offsets(nFineCells.size()+1);
     {
@@ -634,10 +508,8 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
             offsets[i+1] = offsets[i] + nFineCells[i];
         }
     }
-    //Pout<< "offsets:" << offsets << endl;
 
     // Combine and renumber nCoarseCells
-    //Pout<< "Starting nCells_ ..." << endl;
     labelList nCoarseCells;
     gatherList
     (
@@ -648,7 +520,6 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
     );
 
     // (cell)restrictAddressing
-    //Pout<< "Starting restrictAddressing_ ..." << endl;
     const globalIndex cellOffsetter(offsets);
 
     labelList procRestrictAddressing;
@@ -671,15 +542,9 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
                 coarseCellOffsets[i+1] = coarseCellOffsets[i]+nCoarseCells[i];
             }
         }
-        //label nOldCoarseCells = nCells_[levelIndex];
 
         nCells_[levelIndex] = coarseCellOffsets.last();
 
-        //Pout<< "Finished nCoarseCells_ ..."
-        //    << " was:" << nOldCoarseCells
-        //    << " now:" << nCells_[levelIndex] << endl;
-
-
         // Renumber consecutively
         for (label procI = 1; procI < procIDs.size(); procI++)
         {
@@ -695,17 +560,8 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
             }
         }
 
-        //Pout<< "coarseCellOffsets:" << coarseCellOffsets << endl;
         restrictAddressing_[levelIndex].transfer(procRestrictAddressing);
-        //Pout<< "restrictAddressing_:"
-        //    << restrictAddressing_[levelIndex].size()
-        //    << " min:" << min(restrictAddressing_[levelIndex])
-        //    << " max:" << max(restrictAddressing_[levelIndex])
-        //    << endl;
-        //Pout<< "Finished restrictAddressing_ ..." << endl;
     }
-    //Pout<< "** DONE GAMGAgglomeration:procAgglomerateRestrictAddressing **"
-    //    << endl;
 }
 
 
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 bdf1a315f2f05e34f8f966faadb8ae9f7dc1189f..8832fb81a7b5bf9ac04451750a67f120d0c5283d 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
@@ -29,6 +29,7 @@ License
 #include "Time.H"
 #include "dlLibraryTable.H"
 #include "GAMGInterface.H"
+#include "GAMGProcAgglomeration.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -57,10 +58,8 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
 
     // Have procCommunicator_ always, even if not procAgglomerating
     procCommunicator_.setSize(nCreatedLevels + 1);
-    if (processorAgglomerate_)
+    if (processorAgglomerate())
     {
-        Pout<< "GAMGAgglomeration : compacting from " << procAgglomMap_.size()
-            << " to " << nCreatedLevels << " levels" << endl;
         procAgglomMap_.setSize(nCreatedLevels);
         agglomProcIDs_.setSize(nCreatedLevels);
         procCellOffsets_.setSize(nCreatedLevels);
@@ -68,308 +67,7 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
         procBoundaryMap_.setSize(nCreatedLevels);
         procBoundaryFaceMap_.setSize(nCreatedLevels);
 
-
-        if (debug)
-        {
-            Pout<< nl << "Starting mesh overview" << endl;
-            for (label levelI = 0; levelI <= size(); levelI++)
-            {
-                if (hasMeshLevel(levelI))
-                {
-                    const lduMesh& fineMesh = meshLevel(levelI);
-                    const lduInterfacePtrsList& interfaces =
-                        interfaceLevel(levelI);
-
-                    Pout<< "Level " << levelI << " fine mesh:"<< nl
-                        << "    nCells:"
-                        << fineMesh.lduAddr().size() << nl
-                        << "    nFaces:"
-                        << fineMesh.lduAddr().lowerAddr().size() << nl
-                        << "    nInterfaces:" << interfaces.size()
-                        << endl;
-
-                    forAll(interfaces, i)
-                    {
-                        if (interfaces.set(i))
-                        {
-                            Pout<< "        " << i
-                                << "\tsize:" << interfaces[i].faceCells().size()
-                                << endl;
-                        }
-                    }
-                    Pout<< endl;
-                }
-                else
-                {
-                    Pout<< "Level " << levelI << " has no fine mesh:" << nl
-                        << endl;
-                }
-
-                if
-                (
-                    levelI < restrictAddressing_.size()
-                 && restrictAddressing_.set(levelI)
-                )
-                {
-                    const labelList& cellRestrict = restrictAddressing(levelI);
-                    const labelList& faceRestrict =
-                        faceRestrictAddressing(levelI);
-
-                    Pout<< "Level " << levelI << " agglomeration:" << nl
-                        << "    nCoarseCells:" << nCells(levelI) << nl
-                        << "    nCoarseFaces:" << nFaces(levelI) << nl
-                        << "    cellRestriction:"
-                        << " size:" << cellRestrict.size()
-                        << " max:" << max(cellRestrict)
-                        << nl
-                        << "    faceRestriction:"
-                        << " size:" << faceRestrict.size()
-                        << " max:" << max(faceRestrict)
-                        << nl;
-
-
-                    const labelListList& patchFaceRestrict =
-                        patchFaceRestrictAddressing(levelI);
-                    forAll(patchFaceRestrict, i)
-                    {
-                        if (patchFaceRestrict[i].size())
-                        {
-                            const labelList& faceRestrict =
-                                patchFaceRestrict[i];
-                            Pout<< "        " << i
-                                << " size:" << faceRestrict.size()
-                                << " max:" << max(faceRestrict)
-                                << nl;
-                        }
-                    }
-                    Pout<< endl;
-                }
-                Pout<< endl;
-            }
-            Pout<< endl;
-        }
-
-
-
-
-//XXXXXX
-        // As a test: agglomerate finelevel 2, coarselevel 3
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        label fineLevelIndex = 2;
-
-        // Get the coarse mesh
-        const lduMesh& levelMesh = meshLevels_[fineLevelIndex];
-        label levelComm = levelMesh.comm();
-
-        // Processor restriction map: per processor the coarse processor
-        labelList procAgglomMap(UPstream::nProcs(levelComm));
-        {
-            label half = (procAgglomMap.size()+1)/2;
-            for (label i = 0; i < half; i++)
-            {
-                procAgglomMap[i] = 0;
-            }
-            for (label i = half; i < procAgglomMap.size(); i++)
-            {
-                procAgglomMap[i] = 1;
-            }
-        }
-
-        // Master processor
-        labelList masterProcs;
-        // Local processors that agglomerate. agglomProcIDs[0] is in
-        // masterProc.
-        List<int> agglomProcIDs;
-        calculateRegionMaster
-        (
-            levelComm,
-            procAgglomMap,
-            masterProcs,
-            agglomProcIDs
-        );
-
-        Pout<< "procAgglomMap:" << procAgglomMap << endl;
-        Pout<< "agglomProcIDs:" << agglomProcIDs << endl;
-        Pout<< "masterProcs:" << masterProcs << endl;
-
-
-        // Allocate a communicator for the processor-agglomerated matrix
-        label procAgglomComm = UPstream::allocateCommunicator
-        (
-            levelComm,
-            masterProcs
-        );
-
-
-        // Above should really be determined by a procesor-agglomeration
-        // engine. Use procesor agglomeration maps to do the actual
-        // collecting.
-
-        if (Pstream::myProcNo(levelComm) != -1)
-        {
-            // Collect meshes and restrictAddressing onto master
-            // Overwrites the fine mesh (meshLevels_[index-1]) and addressing
-            // from fine mesh to coarse mesh (restrictAddressing_[index]).
-            procAgglomerateLduAddressing
-            (
-                levelComm,
-                procAgglomMap,
-                agglomProcIDs,
-                procAgglomComm,
-
-                fineLevelIndex               //fine level index
-            );
-
-            // Combine restrict addressing only onto master
-            for
-            (
-                label levelI = fineLevelIndex+1;
-                levelI < meshLevels_.size();
-                levelI++
-            )
-            {
-                Pout<< "Starting procAgglomerateRestrictAddressing level:"
-                    << levelI << endl;
-
-                procAgglomerateRestrictAddressing
-                (
-                    levelComm,
-                    agglomProcIDs,
-                    levelI
-                );
-                Pout<< "Finished procAgglomerateRestrictAddressing level:"
-                    << levelI << endl;
-            }
-
-            if (Pstream::myProcNo(levelComm) == agglomProcIDs[0])
-            {
-                // On master. Recreate coarse meshes from restrict addressing
-                for
-                (
-                    label levelI = fineLevelIndex;
-                    levelI < meshLevels_.size();
-                    levelI++
-                )
-                {
-                    Pout<< "Starting agglomerateLduAddressing level:" << levelI
-                        << endl;
-                    agglomerateLduAddressing(levelI);
-                    Pout<< "Finished agglomerateLduAddressing level:" << levelI
-                        << endl;
-                }
-            }
-            else
-            {
-                // Agglomerated away. Clear mesh storage.
-                for
-                (
-                    label levelI = fineLevelIndex+1;
-                    levelI <= size();
-                    levelI++
-                )
-                {
-                    clearLevel(levelI);
-                }
-            }
-        }
-    }
-
-    // Print a bit
-    if (debug)
-    {
-        Pout<< nl << "Mesh overview" << endl;
-        for (label levelI = 0; levelI <= size(); levelI++)
-        {
-            if (hasMeshLevel(levelI))
-            {
-                const lduMesh& fineMesh = meshLevel(levelI);
-                const lduInterfacePtrsList& interfaces =
-                    interfaceLevel(levelI);
-
-                Pout<< "Level " << levelI << " fine mesh:"<< nl
-                    << "    nCells:"
-                    << fineMesh.lduAddr().size() << nl
-                    << "    nFaces:"
-                    << fineMesh.lduAddr().lowerAddr().size() << nl
-                    << "    nInterfaces:" << interfaces.size()
-                    << endl;
-
-                forAll(interfaces, i)
-                {
-                    if (interfaces.set(i))
-                    {
-                        Pout<< "        " << i
-                            << "\tsize:" << interfaces[i].faceCells().size()
-                            << endl;
-                    }
-                }
-
-                Pout<< fineMesh.info() << endl;
-
-                Pout<< endl;
-            }
-            else
-            {
-                Pout<< "Level " << levelI << " has no fine mesh:" << nl
-                    << endl;
-            }
-
-            if
-            (
-                levelI < restrictAddressing_.size()
-             && restrictAddressing_.set(levelI)
-            )
-            {
-                const labelList& cellRestrict = restrictAddressing(levelI);
-                const labelList& faceRestrict =
-                    faceRestrictAddressing(levelI);
-
-                Pout<< "Level " << levelI << " agglomeration:" << nl
-                    << "    nCoarseCells:" << nCells(levelI) << nl
-                    << "    nCoarseFaces:" << nFaces(levelI) << nl
-                    << "    cellRestriction:"
-                    << " size:" << cellRestrict.size()
-                    << " max:" << max(cellRestrict)
-                    << nl
-                    << "    faceRestriction:"
-                    << " size:" << faceRestrict.size()
-                    << " max:" << max(faceRestrict)
-                    << nl;
-
-                const labelListList& patchFaceRestrict =
-                    patchFaceRestrictAddressing(levelI);
-                forAll(patchFaceRestrict, i)
-                {
-                    if (patchFaceRestrict[i].size())
-                    {
-                        const labelList& faceRestrict =
-                            patchFaceRestrict[i];
-                        Pout<< "        " << i
-                            << " size:" << faceRestrict.size()
-                            << " max:" << max(faceRestrict)
-                            << nl;
-                    }
-                }
-            }
-            if
-            (
-                levelI < procCellOffsets_.size()
-             && procCellOffsets_.set(levelI)
-            )
-            {
-                Pout<< "    procCellOffsets:" << procCellOffsets_[levelI]
-                    << nl
-                    << "    procAgglomMap:" << procAgglomMap_[levelI]
-                    << nl
-                    << "    procIDs:" << agglomProcIDs_[levelI]
-                    << nl
-                    << "    comm:" << procCommunicator_[levelI]
-                    << endl;
-            }
-
-            Pout<< endl;
-        }
-        Pout<< endl;
+        procAgglomeratorPtr_().agglomerate();
     }
 }
 
@@ -402,9 +100,16 @@ Foam::GAMGAgglomeration::GAMGAgglomeration
     (
         readLabel(controlDict.lookup("nCellsInCoarsestLevel"))
     ),
-    processorAgglomerate_
+    procAgglomeratorPtr_
     (
-        controlDict.lookupOrDefault<bool>("processorAgglomerate", false)
+        controlDict.found("processorAgglomerator")
+      ? GAMGProcAgglomeration::New
+        (
+            controlDict.lookup("processorAgglomerator"),
+            *this,
+            controlDict
+        )
+      : autoPtr<GAMGProcAgglomeration>(NULL)
     ),
 
     nCells_(maxLevels_),
@@ -420,10 +125,8 @@ Foam::GAMGAgglomeration::GAMGAgglomeration
     interfaceLevels_(maxLevels_ + 1)
 {
     procCommunicator_.setSize(maxLevels_ + 1, -1);
-    if (processorAgglomerate_)
+    if (processorAgglomerate())
     {
-        Pout<< "GAMGAgglomeration : sizing to " << maxLevels_
-            << " levels" << endl;
         procAgglomMap_.setSize(maxLevels_);
         agglomProcIDs_.setSize(maxLevels_);
         procCellOffsets_.setSize(maxLevels_);
@@ -596,8 +299,6 @@ void Foam::GAMGAgglomeration::clearLevel(const label i)
 {
     if (hasMeshLevel(i))
     {
-        Pout<< "Clearing out level " << i << endl;
-
         meshLevels_.set(i - 1, NULL);
 
         if (i < nCells_.size())
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 bcf4b8cec901a1f51a94271d81f73b88db250936..faaf62bbc1ed080c5a864fb2e66440c2b37b2638 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H
@@ -53,6 +53,7 @@ namespace Foam
 class lduMesh;
 class lduMatrix;
 class mapDistribute;
+class GAMGProcAgglomeration;
 
 /*---------------------------------------------------------------------------*\
                     Class GAMGAgglomeration Declaration
@@ -64,16 +65,6 @@ class GAMGAgglomeration
 {
 protected:
 
-    // Protected data types
-
-        //- Enumeration how to handle the mapping of flipped faces
-        enum faceMappingType
-        {
-            IGNORESIGN, // do not apply sign on flipped faces
-            APPLYSIGN,  // apply sign
-            DONTMAP     // do not sign mapped faces
-        };
-
     // Protected data
 
         //- Max number of levels
@@ -82,8 +73,7 @@ protected:
         //- Number of cells in coarsest level
         const label nCellsInCoarsestLevel_;
 
-        //- Whether to agglomerate across processors
-        const bool processorAgglomerate_;
+        autoPtr<GAMGProcAgglomeration> procAgglomeratorPtr_;
 
 
         //- The number of cells in each level
@@ -129,7 +119,7 @@ protected:
         PtrList<lduInterfacePtrsList> interfaceLevels_;
 
 
-        //- Processor agglomeration
+        // Processor agglomeration
 
             //- Per level, per processor the processor it agglomerates into
             mutable PtrList<labelList> procAgglomMap_;
@@ -154,6 +144,8 @@ protected:
             mutable PtrList<labelListListList> procBoundaryFaceMap_;
 
 
+    // Protected Member Functions
+
         //- Assemble coarse mesh addressing
         void agglomerateLduAddressing(const label fineLevelIndex);
 
@@ -163,23 +155,6 @@ protected:
         //- Check the need for further agglomeration
         bool continueAgglomerating(const label nCoarseCells) const;
 
-        //- Collect and combine processor meshes onto procIDs[0]
-        static void gatherMeshes
-        (
-            const label comm,
-            const labelList& procAgglomMap,
-            const labelList& procIDs,
-            const label allMeshComm,
-
-            const lduMesh& myMesh,
-
-            autoPtr<lduPrimitiveMesh>& allMesh,
-            labelList& procCellOffsets,
-            labelListList& procFaceMap,
-            labelListList& procBoundaryMap,
-            labelListListList& procBoundaryFaceMap
-        );
-
         //- Gather value from all procIDs onto procIDs[0]
         template<class Type>
         static void gatherList
@@ -195,6 +170,36 @@ protected:
         void clearLevel(const label leveli);
 
 
+        // Processor agglomeration
+
+            //- Collect and combine processor meshes into allMesh.
+            //  - allMeshComm   : communicator for combined mesh.
+            //  - procAgglomMap : per processor the new agglomerated processor
+            //                    (rank in allMeshComm!). Global information.
+            //  - procIDs       : local information: same for all in
+            //                    agglomerated processor.
+            void procAgglomerateLduAddressing
+            (
+                const label comm,
+                const labelList& procAgglomMap,
+                const labelList& procIDs,
+                const label allMeshComm,
+                const label levelIndex
+            );
+
+            //- Collect and combine basic restriction addressing:
+            //      nCells_
+            //      restrictAddressing_
+            void procAgglomerateRestrictAddressing
+            (
+                const label comm,
+                const labelList& procIDs,
+                const label levelIndex
+            );
+
+
+private:
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -206,6 +211,9 @@ protected:
 
 public:
 
+    //- Declare friendship with GAMGProcAgglomeration
+    friend class GAMGProcAgglomeration;
+
     //- Runtime type information
     TypeName("GAMGAgglomeration");
 
@@ -388,12 +396,6 @@ public:
         // has been agglomerated). This is just for convenience and consistency
         // with GAMGSolver notation.
 
-            //- Whether to agglomerate across processors
-            bool processorAgglomerate() const
-            {
-                return processorAgglomerate_;
-            }
-
             //- Given fine to coarse processor map determine:
             //  - for each coarse processor a master (minimum of the fine
             //    processors)
@@ -407,64 +409,11 @@ public:
                 List<int>& agglomProcIDs
             );
 
-            template<class Type>
-            static void mapField
-            (
-                const faceMappingType mapType,
-                const labelList& map,
-                const UList<Type>& fld,
-                UList<Type>& mappedFld
-            );
-
-            //- Helper: collect and order faces. Requires allFaceFld,
-            //  allPatchFld to have been sized already.
-            template<class Type>
-            static void gatherFaceField
-            (
-                const label comm,
-                const labelList& procIDs,
-
-                const faceMappingType mapType,
-
-                // Field on internal faces
-                const labelListList& procFaceMap,
-                const UList<Type>& faceFld,
-
-                // Field on patches faces
-                const labelListList& procPatchMap,
-                const labelListListList& procPatchFaceMap,
-                const UList<List<Type> >& patchFld,
-
-                List<Type>& allFaceFld,
-                List<List<Type> >& allPatchFld,
-                const int tag = Pstream::msgType()
-            );
-
-
-            //- Collect and combine processor meshes into allMesh.
-            //  - allMeshComm   : communicator for combined mesh.
-            //  - procAgglomMap : per processor the new agglomerated processor
-            //                    (rank in allMeshComm!). Global information.
-            //  - procIDs       : local information: same for all in
-            //                    agglomerated processor.
-            void procAgglomerateLduAddressing
-            (
-                const label comm,
-                const labelList& procAgglomMap,
-                const labelList& procIDs,
-                const label allMeshComm,
-                const label levelIndex
-            );
-
-            //- Collect and combine basic restriction addressing:
-            //      nCells_
-            //      restrictAddressing_
-            void procAgglomerateRestrictAddressing
-            (
-                const label comm,
-                const labelList& procIDs,
-                const label levelIndex
-            );
+            //- Whether to agglomerate across processors
+            bool processorAgglomerate() const
+            {
+                return procAgglomeratorPtr_.valid();
+            }
 
             //- Mapping from processor to agglomerated processor (global, all
             //  processors have the same information). Note that level is
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/GAMGProcAgglomeration/GAMGProcAgglomeration.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/GAMGProcAgglomeration/GAMGProcAgglomeration.C
new file mode 100644
index 0000000000000000000000000000000000000000..b8a1e59a74522138ccb6cff6815c0c81f66b03c3
--- /dev/null
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/GAMGProcAgglomeration/GAMGProcAgglomeration.C
@@ -0,0 +1,274 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 "GAMGProcAgglomeration.H"
+#include "GAMGAgglomeration.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(GAMGProcAgglomeration, 0);
+    defineRunTimeSelectionTable(GAMGProcAgglomeration, GAMGAgglomeration);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::GAMGProcAgglomeration::printStats
+(
+    Ostream& os,
+    GAMGAgglomeration& agglom
+) const
+{
+    for (label levelI = 0; levelI <= agglom.size(); levelI++)
+    {
+        if (agglom.hasMeshLevel(levelI))
+        {
+            const lduMesh& fineMesh = agglom.meshLevel(levelI);
+            const lduInterfacePtrsList& interfaces =
+                agglom.interfaceLevel(levelI);
+
+            os  << "Level " << levelI << " fine mesh:"<< nl
+                << "    nCells:"
+                << fineMesh.lduAddr().size() << nl
+                << "    nFaces:"
+                << fineMesh.lduAddr().lowerAddr().size() << nl
+                << "    nInterfaces:" << interfaces.size()
+                << endl;
+
+            forAll(interfaces, i)
+            {
+                if (interfaces.set(i))
+                {
+                    os  << "        " << i
+                        << "\tsize:" << interfaces[i].faceCells().size()
+                        << endl;
+                }
+            }
+
+            os  << fineMesh.info() << endl;
+
+            os  << endl;
+        }
+        else
+        {
+            os  << "Level " << levelI << " has no fine mesh:" << nl
+                << endl;
+        }
+
+        if
+        (
+            levelI < agglom.restrictAddressing_.size()
+         && agglom.restrictAddressing_.set(levelI)
+        )
+        {
+            const labelList& cellRestrict =
+                agglom.restrictAddressing(levelI);
+            const labelList& faceRestrict =
+                agglom.faceRestrictAddressing(levelI);
+
+            os  << "Level " << levelI << " agglomeration:" << nl
+                << "    nCoarseCells:" << agglom.nCells(levelI) << nl
+                << "    nCoarseFaces:" << agglom.nFaces(levelI) << nl
+                << "    cellRestriction:"
+                << " size:" << cellRestrict.size()
+                << " max:" << max(cellRestrict)
+                << nl
+                << "    faceRestriction:"
+                << " size:" << faceRestrict.size()
+                << " max:" << max(faceRestrict)
+                << nl;
+
+            const labelListList& patchFaceRestrict =
+                agglom.patchFaceRestrictAddressing(levelI);
+            forAll(patchFaceRestrict, i)
+            {
+                if (patchFaceRestrict[i].size())
+                {
+                    const labelList& faceRestrict =
+                        patchFaceRestrict[i];
+                    os  << "        " << i
+                        << " size:" << faceRestrict.size()
+                        << " max:" << max(faceRestrict)
+                        << nl;
+                }
+            }
+        }
+        if
+        (
+            levelI < agglom.procCellOffsets_.size()
+         && agglom.procCellOffsets_.set(levelI)
+        )
+        {
+            os  << "    procCellOffsets:" << agglom.procCellOffsets_[levelI]
+                << nl
+                << "    procAgglomMap:" << agglom.procAgglomMap_[levelI]
+                << nl
+                << "    procIDs:" << agglom.agglomProcIDs_[levelI]
+                << nl
+                << "    comm:" << agglom.procCommunicator_[levelI]
+                << endl;
+        }
+
+        os  << endl;
+    }
+    os  << endl;
+}
+
+
+bool Foam::GAMGProcAgglomeration::agglomerate
+(
+    const label fineLevelIndex,
+    const labelList& procAgglomMap,
+    const labelList& masterProcs,
+    const List<int>& agglomProcIDs,
+    const label procAgglomComm
+)
+{
+    const lduMesh& levelMesh = agglom_.meshLevels_[fineLevelIndex];
+    label levelComm = levelMesh.comm();
+
+    if (Pstream::myProcNo(levelComm) != -1)
+    {
+        // Collect meshes and restrictAddressing onto master
+        // Overwrites the fine mesh (meshLevels_[index-1]) and addressing
+        // from fine mesh to coarse mesh (restrictAddressing_[index]).
+        agglom_.procAgglomerateLduAddressing
+        (
+            levelComm,
+            procAgglomMap,
+            agglomProcIDs,
+            procAgglomComm,
+
+            fineLevelIndex               //fine level index
+        );
+
+        // Combine restrict addressing only onto master
+        for
+        (
+            label levelI = fineLevelIndex+1;
+            levelI < agglom_.meshLevels_.size();
+            levelI++
+        )
+        {
+            agglom_.procAgglomerateRestrictAddressing
+            (
+                levelComm,
+                agglomProcIDs,
+                levelI
+            );
+        }
+
+        if (Pstream::myProcNo(levelComm) == agglomProcIDs[0])
+        {
+            // On master. Recreate coarse meshes from restrict addressing
+            for
+            (
+                label levelI = fineLevelIndex;
+                levelI < agglom_.meshLevels_.size();
+                levelI++
+            )
+            {
+                agglom_.agglomerateLduAddressing(levelI);
+            }
+        }
+        else
+        {
+            // Agglomerated away. Clear mesh storage.
+            for
+            (
+                label levelI = fineLevelIndex+1;
+                levelI <= agglom_.size();
+                levelI++
+            )
+            {
+                agglom_.clearLevel(levelI);
+            }
+        }
+    }
+
+    // Should check!
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::GAMGProcAgglomeration::GAMGProcAgglomeration
+(
+    GAMGAgglomeration& agglom,
+    const dictionary& controlDict
+)
+:
+    agglom_(agglom)
+{}
+
+
+Foam::autoPtr<Foam::GAMGProcAgglomeration> Foam::GAMGProcAgglomeration::New
+(
+    const word& type,
+    GAMGAgglomeration& agglom,
+    const dictionary& controlDict
+)
+{
+    if (debug)
+    {
+        Info<< "GAMGProcAgglomeration::New(const word&, GAMGAgglomeration&"
+               ", const dictionary&) : "
+               "constructing GAMGProcAgglomeration"
+            << endl;
+    }
+
+    GAMGAgglomerationConstructorTable::iterator cstrIter =
+        GAMGAgglomerationConstructorTablePtr_->find(type);
+
+    if (cstrIter == GAMGAgglomerationConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "GAMGProcAgglomeration::New(const word&, GAMGAgglomeration&"
+            ", const dictionary&) "
+        )   << "Unknown GAMGProcAgglomeration type "
+            << type << " for GAMGAgglomeration " << agglom.type() << nl << nl
+            << "Valid GAMGProcAgglomeration types are :" << endl
+            << GAMGAgglomerationConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<GAMGProcAgglomeration>(cstrIter()(agglom, controlDict));
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::GAMGProcAgglomeration::~GAMGProcAgglomeration()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/GAMGProcAgglomeration/GAMGProcAgglomeration.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/GAMGProcAgglomeration/GAMGProcAgglomeration.H
new file mode 100644
index 0000000000000000000000000000000000000000..cf27f44bc0fbe2b5394fb92009e5024e142f0359
--- /dev/null
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/GAMGProcAgglomeration/GAMGProcAgglomeration.H
@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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::GAMGProcAgglomeration
+
+Description
+    Processor agglomeration of GAMGAgglomerations.
+
+SourceFiles
+    GAMGProcAgglomeration.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef GAMGProcAgglomeration_H
+#define GAMGProcAgglomeration_H
+
+#include "runTimeSelectionTables.H"
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class GAMGAgglomeration;
+
+/*---------------------------------------------------------------------------*\
+                    Class GAMGProcAgglomeration Declaration
+\*---------------------------------------------------------------------------*/
+
+class GAMGProcAgglomeration
+{
+
+protected:
+
+    // Protected data
+
+        //- Reference to agglomeration
+        GAMGAgglomeration& agglom_;
+
+    // Protected Member Functions
+
+        //- Debug: write agglomeration info
+        void printStats(Ostream& os, GAMGAgglomeration& agglom) const;
+
+        //- Agglomerate a level. Return true if anything has changed
+        bool agglomerate
+        (
+            const label fineLevelIndex,
+            const labelList& procAgglomMap,
+            const labelList& masterProcs,
+            const List<int>& agglomProcIDs,
+            const label procAgglomComm
+        );
+
+private:
+
+    // Private data
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        GAMGProcAgglomeration(const GAMGProcAgglomeration&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const GAMGProcAgglomeration&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("GAMGProcAgglomeration");
+
+
+    // Declare run-time constructor selection tables
+
+        //- Runtime selection table for pure geometric agglomerators
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            GAMGProcAgglomeration,
+            GAMGAgglomeration,
+            (
+                GAMGAgglomeration& agglom,
+                const dictionary& controlDict
+            ),
+            (
+                agglom,
+                controlDict
+            )
+        );
+
+
+    // Constructors
+
+        //- Construct given agglomerator and controls
+        GAMGProcAgglomeration
+        (
+            GAMGAgglomeration& agglom,
+            const dictionary& controlDict
+        );
+
+
+    // Selectors
+
+        //- Return the selected agglomerator
+        static autoPtr<GAMGProcAgglomeration> New
+        (
+            const word& type,
+            GAMGAgglomeration& agglom,
+            const dictionary& controlDict
+        );
+
+
+    //- Destructor
+    virtual ~GAMGProcAgglomeration();
+
+
+    // Member Functions
+
+       //- Modify agglomeration. Return true if modified
+        virtual bool agglomerate() = 0;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/masterCoarsest/masterCoarsest.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/masterCoarsest/masterCoarsest.C
new file mode 100644
index 0000000000000000000000000000000000000000..6a88ef4112a9993411bfbd89c4e459982f8caf28
--- /dev/null
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/masterCoarsest/masterCoarsest.C
@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 "masterCoarsest.H"
+#include "addToRunTimeSelectionTable.H"
+#include "GAMGAgglomeration.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(masterCoarsest, 0);
+
+    addToRunTimeSelectionTable
+    (
+        GAMGProcAgglomeration,
+        masterCoarsest,
+        GAMGAgglomeration
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::masterCoarsest::masterCoarsest
+(
+    GAMGAgglomeration& agglom,
+    const dictionary& controlDict
+)
+:
+    GAMGProcAgglomeration(agglom, controlDict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::masterCoarsest::agglomerate()
+{
+    if (debug)
+    {
+        Pout<< nl << "Starting mesh overview" << endl;
+        printStats(Pout, agglom_);
+    }
+
+    if (agglom_.size() >= 1)
+    {
+        // Agglomerate one but last level (since also agglomerating
+        // restrictAddressing)
+        label fineLevelIndex = agglom_.size()-1;
+
+        // Get the coarse mesh
+        const lduMesh& levelMesh = agglom_.meshLevel(fineLevelIndex-1);
+        label levelComm = levelMesh.comm();
+
+        // Processor restriction map: per processor the coarse processor
+        labelList procAgglomMap(UPstream::nProcs(levelComm));
+        //{
+        //    label half = (procAgglomMap.size()+1)/2;
+        //    for (label i = 0; i < half; i++)
+        //    {
+        //        procAgglomMap[i] = 0;
+        //    }
+        //    for (label i = half; i < procAgglomMap.size(); i++)
+        //    {
+        //        procAgglomMap[i] = 1;
+        //    }
+        //}
+        procAgglomMap = 0;
+
+        // Master processor
+        labelList masterProcs;
+        // Local processors that agglomerate. agglomProcIDs[0] is in
+        // masterProc.
+        List<int> agglomProcIDs;
+        GAMGAgglomeration::calculateRegionMaster
+        (
+            levelComm,
+            procAgglomMap,
+            masterProcs,
+            agglomProcIDs
+        );
+
+        // Allocate a communicator for the processor-agglomerated matrix
+        label procAgglomComm = UPstream::allocateCommunicator
+        (
+            levelComm,
+            masterProcs
+        );
+
+
+        // Above should really be determined by a procesor-agglomeration
+        // engine. Use procesor agglomeration maps to do the actual
+        // collecting.
+
+        if (Pstream::myProcNo(levelComm) != -1)
+        {
+            GAMGProcAgglomeration::agglomerate
+            (
+                fineLevelIndex,
+                procAgglomMap,
+                masterProcs,
+                agglomProcIDs,
+                procAgglomComm
+            );
+        }
+    }
+
+    // Print a bit
+    if (debug)
+    {
+        Pout<< nl << "Agglomerated mesh overview" << endl;
+        printStats(Pout, agglom_);
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/masterCoarsest/masterCoarsest.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/masterCoarsest/masterCoarsest.H
new file mode 100644
index 0000000000000000000000000000000000000000..1b7828612cadba6cdb646a21b91256129c0c94b4
--- /dev/null
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/masterCoarsest/masterCoarsest.H
@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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::masterCoarsest
+
+Description
+    Processor agglomeration of GAMGAgglomerations.
+
+SourceFiles
+    masterCoarsest.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef masterCoarsest_H
+#define masterCoarsest_H
+
+#include "GAMGProcAgglomeration.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class GAMGAgglomeration;
+
+/*---------------------------------------------------------------------------*\
+                    Class masterCoarsest Declaration
+\*---------------------------------------------------------------------------*/
+
+class masterCoarsest
+:
+    public GAMGProcAgglomeration
+{
+    // Private data
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        masterCoarsest(const masterCoarsest&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const masterCoarsest&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("masterCoarsest");
+
+
+    // Constructors
+
+        //- Construct given agglomerator and controls
+        masterCoarsest
+        (
+            GAMGAgglomeration& agglom,
+            const dictionary& controlDict
+        );
+
+
+    //- Destructor
+    virtual ~masterCoarsest()
+    {}
+
+
+    // Member Functions
+
+       //- Modify agglomeration. Return true if modified
+        virtual bool agglomerate();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//#ifdef NoRepository
+//#   include "masterCoarsestTemplates.C"
+//#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //