diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
index be850e59b0..ceb0661489 100644
--- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
+++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
@@ -29,6 +29,14 @@ License
 #include "procLduInterface.H"
 #include "cyclicLduInterface.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(LUscalarMatrix, 0);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::LUscalarMatrix::LUscalarMatrix(const scalarSquareMatrix& matrix)
@@ -134,8 +142,39 @@ Foam::LUscalarMatrix::LUscalarMatrix
         convert(ldum, interfaceCoeffs, interfaces);
     }
 
-    if (Pstream::master(comm_))
+    if (debug && Pstream::master(comm_))
     {
+        label nRows = n();
+        label nColumns = m();
+
+        Pout<< "LUscalarMatrix : size:" << nRows << endl;
+        for (label rowI = 0; rowI < nRows; rowI++)
+        {
+            const scalar* row = operator[](rowI);
+
+            Pout<< "cell:" << rowI << " diagCoeff:" << row[rowI] << endl;
+
+            Pout<< "    connects to upper cells :";
+            for (label columnI = rowI+1; columnI < nColumns; columnI++)
+            {
+                if (mag(row[columnI]) > SMALL)
+                {
+                    Pout<< ' ' << columnI << " (coeff:" << row[columnI] << ")";
+                }
+            }
+            Pout<< endl;
+            Pout<< "    connects to lower cells :";
+            for (label columnI = 0; columnI < rowI; columnI++)
+            {
+                if (mag(row[columnI]) > SMALL)
+                {
+                    Pout<< ' ' << columnI << " (coeff:" << row[columnI] << ")";
+                }
+            }
+            Pout<< endl;
+        }
+        Pout<< endl;
+
         pivotIndices_.setSize(n());
         LUDecompose(*this, pivotIndices_);
     }
diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
index b01c88224e..9c96d8dd23 100644
--- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
+++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
@@ -87,6 +87,9 @@ class LUscalarMatrix
 
 public:
 
+    // Declare name of the class and its debug switch
+    ClassName("LUscalarMatrix");
+
     // Constructors
 
         //- Construct from scalarSquareMatrix and perform LU decomposition
diff --git a/src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.C b/src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.C
index 11d27fabf8..60f5148ad7 100644
--- a/src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.C
+++ b/src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -99,8 +99,20 @@ void Foam::GAMGPreconditioner::precondition
     // Create the smoothers for all levels
     PtrList<lduMatrix::smoother> smoothers;
 
+    // Scratch fields if processor-agglomerated coarse level meshes
+    // are bigger than original. Usually not needed
+    scalarField ApsiScratch;
+    scalarField finestCorrectionScratch;
+
     // Initialise the above data structures
-    initVcycle(coarseCorrFields, coarseSources, smoothers);
+    initVcycle
+    (
+        coarseCorrFields,
+        coarseSources,
+        smoothers,
+        ApsiScratch,
+        finestCorrectionScratch
+    );
 
     for (label cycle=0; cycle<nVcycles_; cycle++)
     {
@@ -112,6 +124,14 @@ void Foam::GAMGPreconditioner::precondition
             AwA,
             finestCorrection,
             finestResidual,
+
+            (ApsiScratch.size() ? ApsiScratch : AwA),
+            (
+                finestCorrectionScratch.size()
+              ? finestCorrectionScratch
+              : finestCorrection
+            ),
+
             coarseCorrFields,
             coarseSources,
             cmpt
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 e8cecab5e4..6d806b7789 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C
@@ -270,12 +270,19 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
         interfaceLevels_[fineLevelIndex];
 
     // Create coarse-level interfaces
+    primitiveInterfaces_.set
+    (
+        fineLevelIndex + 1,
+        new PtrList<const lduInterface>(fineInterfaces.size())
+    );
+    PtrList<const lduInterface>& coarsePrimInterfaces =
+        primitiveInterfaces_[fineLevelIndex + 1];
+
     interfaceLevels_.set
     (
         fineLevelIndex + 1,
         new lduInterfacePtrsList(fineInterfaces.size())
     );
-
     lduInterfacePtrsList& coarseInterfaces =
         interfaceLevels_[fineLevelIndex + 1];
 
@@ -317,7 +324,7 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
     {
         if (fineInterfaces.set(inti))
         {
-            coarseInterfaces.set
+            coarsePrimInterfaces.set
             (
                 inti,
                 GAMGInterface::New
@@ -336,6 +343,7 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
                 ).ptr()
             );
 
+            coarseInterfaces.set(inti, &coarsePrimInterfaces[inti]);
             coarseInterfaceAddr[inti] = coarseInterfaces[inti].faceCells();
             nPatchFaces[inti] = coarseInterfaceAddr[inti].size();
             patchFineToCoarse[inti] = refCast<const GAMGInterface>
@@ -408,13 +416,13 @@ void Foam::GAMGAgglomeration::gatherMeshes
         // Combine all addressing
         // ~~~~~~~~~~~~~~~~~~~~~~
 
-        Pout<< "Own Mesh " << myMesh.lduAddr().size() << endl;
-        forAll(otherMeshes, i)
-        {
-            Pout<< "    otherMesh " << i << " "
-                << otherMeshes[i].lduAddr().size()
-                << endl;
-        }
+        //Pout<< "Own Mesh " << myMesh.lduAddr().size() << endl;
+        //forAll(otherMeshes, i)
+        //{
+        //    Pout<< "    otherMesh " << i << " "
+        //        << otherMeshes[i].lduAddr().size()
+        //        << endl;
+        //}
 
         labelList procFaceOffsets;
 
@@ -437,7 +445,7 @@ void Foam::GAMGAgglomeration::gatherMeshes
             )
         );
 
-        Pout<< "** Agglomerated Mesh " << allMesh().info();
+        //Pout<< "** Agglomerated Mesh " << allMesh().info();
     }
 
     UPstream::warnComm = oldWarn;
@@ -459,7 +467,6 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
         << " havefinemesh:" << meshLevels_.set(levelIndex-1) << endl;
 
     const lduMesh& myMesh = meshLevels_[levelIndex-1];
-    //label meshComm = myMesh.comm();
 
     Pout<< "my mesh                     :" << myMesh.info();
     Pout<< "nCoarseCells                :" << nCells_[levelIndex] << endl;
@@ -510,8 +517,6 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
     // These could only be set on the master procs but it is
     // quite convenient to also have them on the slaves
     procCellOffsets_.set(levelIndex, new labelList(0));
-    //procRestrictAddressing_.set(levelIndex, new labelField(0));
-    //procFaceRestrictAddressing_.set(levelIndex, new labelList(0));
     procFaceMap_.set(levelIndex, new labelListList(0));
     procBoundaryMap_.set(levelIndex, new labelListList(0));
     procBoundaryFaceMap_.set(levelIndex, new labelListListList(0));
@@ -575,46 +580,16 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
     }
     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<< "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;
-        //    }
-        //}
-        Pout<< "** DONE GAMGAgglomeration:procAgglomerateLduAddressing **"
-            << endl;
+        //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;
@@ -628,14 +603,15 @@ 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<< "** 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;
+    //Pout<< "fine:" << restrictAddressing_[levelIndex].size() << endl;
 
     // Collect number of cells
     labelList nFineCells;
@@ -648,7 +624,7 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
     );
     //const labelList& offsets = procCellOffsets_[levelIndex];
     //Pout<< "offsets:" << offsets << endl;
-    Pout<< "nFineCells:" << nFineCells << endl;
+    //Pout<< "nFineCells:" << nFineCells << endl;
 
     labelList offsets(nFineCells.size()+1);
     {
@@ -658,10 +634,10 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
             offsets[i+1] = offsets[i] + nFineCells[i];
         }
     }
-    Pout<< "offsets:" << offsets << endl;
+    //Pout<< "offsets:" << offsets << endl;
 
     // Combine and renumber nCoarseCells
-    Pout<< "Starting nCells_ ..." << endl;
+    //Pout<< "Starting nCells_ ..." << endl;
     labelList nCoarseCells;
     gatherList
     (
@@ -672,7 +648,7 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
     );
 
     // (cell)restrictAddressing
-    Pout<< "Starting restrictAddressing_ ..." << endl;
+    //Pout<< "Starting restrictAddressing_ ..." << endl;
     const globalIndex cellOffsetter(offsets);
 
     labelList procRestrictAddressing;
@@ -695,12 +671,13 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
                 coarseCellOffsets[i+1] = coarseCellOffsets[i]+nCoarseCells[i];
             }
         }
-        label nOldCoarseCells = nCells_[levelIndex];
+        //label nOldCoarseCells = nCells_[levelIndex];
 
         nCells_[levelIndex] = coarseCellOffsets.last();
-        Pout<< "Finished nCoarseCells_ ..."
-            << " was:" << nOldCoarseCells
-            << " now:" << nCells_[levelIndex] << endl;
+
+        //Pout<< "Finished nCoarseCells_ ..."
+        //    << " was:" << nOldCoarseCells
+        //    << " now:" << nCells_[levelIndex] << endl;
 
 
         // Renumber consecutively
@@ -718,17 +695,17 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
             }
         }
 
-        Pout<< "coarseCellOffsets:" << coarseCellOffsets << endl;
+        //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<< "restrictAddressing_:"
+        //    << restrictAddressing_[levelIndex].size()
+        //    << " min:" << min(restrictAddressing_[levelIndex])
+        //    << " max:" << max(restrictAddressing_[levelIndex])
+        //    << endl;
+        //Pout<< "Finished restrictAddressing_ ..." << endl;
     }
-    Pout<< "** DONE GAMGAgglomeration:procAgglomerateRestrictAddressing **"
-        << endl;
+    //Pout<< "** DONE GAMGAgglomeration:procAgglomerateRestrictAddressing **"
+    //    << endl;
 }
 
 
@@ -777,4 +754,55 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
 //}
 
 
+void Foam::GAMGAgglomeration::calculateRegionMaster
+(
+    const label comm,
+    const labelList& procAgglomMap,
+    labelList& masterProcs,
+    List<int>& agglomProcIDs
+)
+{
+    // 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() = min(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(comm);
+    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]);
+}
+
+
+
 // ************************************************************************* //
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 d7f490e2c4..bdf1a315f2 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
@@ -52,6 +52,7 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
     nPatchFaces_.setSize(nCreatedLevels);
     patchFaceRestrictAddressing_.setSize(nCreatedLevels);
     meshLevels_.setSize(nCreatedLevels);
+    primitiveInterfaces_.setSize(nCreatedLevels + 1);
     interfaceLevels_.setSize(nCreatedLevels + 1);
 
     // Have procCommunicator_ always, even if not procAgglomerating
@@ -62,9 +63,6 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
             << " to " << nCreatedLevels << " levels" << endl;
         procAgglomMap_.setSize(nCreatedLevels);
         agglomProcIDs_.setSize(nCreatedLevels);
-        //procMeshLevels_.setSize(nCreatedLevels);
-        //procRestrictAddressing_.setSize(nCreatedLevels);
-        //procFaceRestrictAddressing_.setSize(nCreatedLevels);
         procCellOffsets_.setSize(nCreatedLevels);
         procFaceMap_.setSize(nCreatedLevels);
         procBoundaryMap_.setSize(nCreatedLevels);
@@ -165,60 +163,34 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
 
         // 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
+        );
 
-        {
-            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
-            agglomProcIDs = findIndices(procAgglomMap, myAgglom);
-            // Make sure the master is the first element.
-            label index = findIndex
-            (
-                agglomProcIDs,
-                agglomToMaster[myAgglom]
-            );
-            Swap(agglomProcIDs[0], agglomProcIDs[index]);
-        }
         Pout<< "procAgglomMap:" << procAgglomMap << endl;
         Pout<< "agglomProcIDs:" << agglomProcIDs << endl;
+        Pout<< "masterProcs:" << masterProcs << endl;
 
 
         // Allocate a communicator for the processor-agglomerated matrix
@@ -259,10 +231,6 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
                 Pout<< "Starting procAgglomerateRestrictAddressing level:"
                     << levelI << endl;
 
-                //procAgglomMap_.set(levelI, new labelList(procAgglomMap));
-                //agglomProcIDs_.set(levelI, new labelList(agglomProcIDs));
-                //procCommunicator_[levelI] = procAgglomComm;
-
                 procAgglomerateRestrictAddressing
                 (
                     levelComm,
@@ -335,6 +303,9 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
                             << endl;
                     }
                 }
+
+                Pout<< fineMesh.info() << endl;
+
                 Pout<< endl;
             }
             else
@@ -445,6 +416,7 @@ Foam::GAMGAgglomeration::GAMGAgglomeration
     patchFaceRestrictAddressing_(maxLevels_),
 
     meshLevels_(maxLevels_),
+    primitiveInterfaces_(maxLevels_ + 1),
     interfaceLevels_(maxLevels_ + 1)
 {
     procCommunicator_.setSize(maxLevels_ + 1, -1);
@@ -454,9 +426,6 @@ Foam::GAMGAgglomeration::GAMGAgglomeration
             << " levels" << endl;
         procAgglomMap_.setSize(maxLevels_);
         agglomProcIDs_.setSize(maxLevels_);
-        //procMeshLevels_.setSize(maxLevels_);
-        //procRestrictAddressing_.setSize(maxLevels_);
-        //procFaceRestrictAddressing_.setSize(maxLevels_);
         procCellOffsets_.setSize(maxLevels_);
         procFaceMap_.setSize(maxLevels_);
         procBoundaryMap_.setSize(maxLevels_);
@@ -573,42 +542,6 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
 
 Foam::GAMGAgglomeration::~GAMGAgglomeration()
 {
-//    // Temporary store the user-defined communicators so we can delete them
-//    labelHashSet communicators(meshLevels_.size());
-//    forAll(meshLevels_, leveli)
-//    {
-//        communicators.insert(meshLevels_[leveli].comm());
-//    }
-    //forAll(procMeshLevels_, leveli)
-    //{
-    //    if (procMeshLevels_.set(leveli))
-    //    {
-    //        communicators.insert(procMeshLevels_[leveli].comm());
-    //    }
-    //}
-
-//    Pout<< "~GAMGAgglomeration() : current communicators:" << communicators
-//        << endl;
-
-    // Clear the interface storage by hand.
-    // It is a list of ptrs not a PtrList for consistency of the interface
-    for (label leveli=1; leveli<interfaceLevels_.size(); leveli++)
-    {
-        lduInterfacePtrsList& curLevel = interfaceLevels_[leveli];
-
-        forAll(curLevel, i)
-        {
-            if (curLevel.set(i))
-            {
-                delete curLevel(i);
-            }
-        }
-    }
-
-//    forAllConstIter(labelHashSet, communicators, iter)
-//    {
-//        UPstream::freeCommunicator(iter.key());
-//    }
     forAllReverse(procCommunicator_, i)
     {
         if (procCommunicator_[i] != -1)
@@ -676,7 +609,8 @@ void Foam::GAMGAgglomeration::clearLevel(const label i)
             faceFlipMap_.set(i, NULL);
             nPatchFaces_.set(i, NULL);
             patchFaceRestrictAddressing_.set(i, NULL);
-            interfaceLevels_.set(i, NULL);  // TBD: delete storage
+            primitiveInterfaces_.set(i, NULL);
+            interfaceLevels_.set(i, NULL);
         }
     }
 }
@@ -706,13 +640,6 @@ bool Foam::GAMGAgglomeration::hasProcMesh(const label leveli) const
 }
 
 
-//const Foam::lduMesh& Foam::GAMGAgglomeration::procMeshLevel
-//(const label leveli) const
-//{
-//    return procMeshLevels_[leveli];
-//}
-
-
 Foam::label Foam::GAMGAgglomeration::procCommunicator(const label leveli) const
 {
     return procCommunicator_[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 d8baab09d9..8590cff0c2 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H
@@ -62,8 +62,7 @@ class GAMGAgglomeration
 :
     public MeshObject<lduMesh, GAMGAgglomeration>
 {
-//protected:
-public:
+protected:
 
     // Protected data types
 
@@ -86,6 +85,7 @@ public:
         //- Whether to agglomerate across processors
         const bool processorAgglomerate_;
 
+
         //- The number of cells in each level
         labelList nCells_;
 
@@ -122,8 +122,10 @@ public:
         //- Hierarchy of mesh addressing
         PtrList<lduPrimitiveMesh> meshLevels_;
 
-        //- Hierarchy interfaces.
-        //  Warning: Needs to be deleted explicitly.
+        //- Hierarchy of interfaces
+        PtrList<PtrList<const lduInterface> > primitiveInterfaces_;
+
+        //- Hierarchy of interfaces in lduInterfacePtrsList form
         PtrList<lduInterfacePtrsList> interfaceLevels_;
 
 
@@ -136,13 +138,6 @@ public:
             //  the 'master' of the cluster.
             mutable PtrList<labelList> agglomProcIDs_;
 
-            //- Combined mesh for all processors
-            //mutable PtrList<lduPrimitiveMesh> procMeshLevels_;
-
-            //mutable PtrList<labelField> procRestrictAddressing_;
-
-            //mutable PtrList<labelList> procFaceRestrictAddressing_;
-
             //- Communicator for given level
             mutable labelList procCommunicator_;
 
@@ -185,17 +180,6 @@ public:
             labelListListList& procBoundaryFaceMap
         );
 
-        ////- Gather value from all procIDs onto procIDs[0]
-        //static void gatherList
-        //(
-        //    const label comm,
-        //    const labelList& procIDs,
-        //
-        //    const label myVal,
-        //    labelList& vals,
-        //    const int tag = Pstream::msgType()
-        //);
-
         //- Gather value from all procIDs onto procIDs[0]
         template<class Type>
         static void gatherList
@@ -420,7 +404,18 @@ public:
                 return processorAgglomerate_;
             }
 
-
+            //- Given fine to coarse processor map determine:
+            //  - for each coarse processor a master (minimum of the fine
+            //    processors)
+            //  - for each coarse processor the set of fine processors
+            //    (element 0 is the master processor)
+            static void calculateRegionMaster
+            (
+                const label comm,
+                const labelList& procAgglomMap,
+                labelList& masterProcs,
+                List<int>& agglomProcIDs
+            );
 
             template<class Type>
             static void mapField
@@ -488,18 +483,14 @@ public:
             const labelList& procAgglomMap(const label fineLeveli) const;
 
             //- Set of processors to agglomerate. Element 0 is the
-            //  master processor. (local, same only on those processor that
+            //  master processor. (local, same only on those processors that
             //  agglomerate)
             const labelList& agglomProcIDs(const label fineLeveli) const;
 
             //- Check that level has combined mesh
             bool hasProcMesh(const label fineLeveli) const;
 
-            //- Combined mesh
-            //const lduMesh& procMeshLevel(const label leveli) const;
-
-            //- Communicator for procMesh (stored separately so also works
-            //  even if no mesh is stored on this processor)
+            //- Communicator for current level or -1
             label procCommunicator(const label fineLeveli) const;
 
             //- Mapping from processor to procMesh cells
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C
index 18b0f90269..2fb1cbe39b 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C
@@ -78,22 +78,13 @@ Foam::GAMGSolver::GAMGSolver
     agglomeration_(GAMGAgglomeration::New(matrix_, controlDict_)),
 
     matrixLevels_(agglomeration_.size()),
-    interfaceLevels_(agglomeration_.size()),
     primitiveInterfaceLevels_(agglomeration_.size()),
+    interfaceLevels_(agglomeration_.size()),
     interfaceLevelsBouCoeffs_(agglomeration_.size()),
     interfaceLevelsIntCoeffs_(agglomeration_.size())
 {
     readControls();
 
-    //if (agglomeration_.processorAgglomerate())
-    //{
-    //    procMatrixLevels_.setSize(agglomeration_.size());
-    //    procInterfaceLevels_.setSize(agglomeration_.size());
-    //    procPrimitiveInterfaces_.setSize(agglomeration_.size());
-    //    procInterfaceLevelsBouCoeffs_.setSize(agglomeration_.size());
-    //    procInterfaceLevelsIntCoeffs_.setSize(agglomeration_.size());
-    //}
-
     if (agglomeration_.processorAgglomerate())
     {
         forAll(agglomeration_, fineLevelIndex)
@@ -105,15 +96,9 @@ Foam::GAMGSolver::GAMGSolver
 
                 if
                 (
-                    //   !agglomeration_.hasMeshLevel(fineLevelIndex+1)
-                    //||  (
-                    //        agglomeration_.nCells(fineLevelIndex)
-                    //     != agglomeration_.meshLevel(fineLevelIndex+1).
-                    //            lduAddr().size()
-                    //    )
-
-                    (fineLevelIndex+1) < agglomeration_.procBoundaryMap_.size()
-                 && agglomeration_.procBoundaryMap_.set(fineLevelIndex+1)
+                    (fineLevelIndex+1) < agglomeration_.size()
+                 //&& agglomeration_.procBoundaryMap_.set(fineLevelIndex+1)
+                 && agglomeration_.hasProcMesh(fineLevelIndex+1)
                 )
                 {
                     Pout<< "Level:" << fineLevelIndex
@@ -291,29 +276,33 @@ Foam::GAMGSolver::GAMGSolver
             Pout<< "GAMGSolver :"
                 << " coarsestLevel:" << coarsestLevel << endl;
 
-            const lduMesh& coarsestMesh = matrixLevels_[coarsestLevel].mesh();
+            if (matrixLevels_.set(coarsestLevel))
+            {
+                const lduMesh& coarsestMesh =
+                    matrixLevels_[coarsestLevel].mesh();
 
-            label coarseComm = coarsestMesh.comm();
-            label oldWarn = UPstream::warnComm;
-            UPstream::warnComm = coarseComm;
+                label coarseComm = coarsestMesh.comm();
+                label oldWarn = UPstream::warnComm;
+                UPstream::warnComm = coarseComm;
 
-            Pout<< "Solve direct on coasestmesh (level=" << coarsestLevel
-                << ") using communicator " << coarseComm << endl;
+                Pout<< "Solve direct on coasestmesh (level=" << coarsestLevel
+                    << ") using communicator " << coarseComm << endl;
 
-            coarsestLUMatrixPtr_.set
-            (
-                new LUscalarMatrix
+                coarsestLUMatrixPtr_.set
                 (
-                    matrixLevels_[coarsestLevel],
-                    interfaceLevelsBouCoeffs_[coarsestLevel],
-                    interfaceLevels_[coarsestLevel]
-                )
-            );
+                    new LUscalarMatrix
+                    (
+                        matrixLevels_[coarsestLevel],
+                        interfaceLevelsBouCoeffs_[coarsestLevel],
+                        interfaceLevels_[coarsestLevel]
+                    )
+                );
 
-            UPstream::warnComm = oldWarn;
+                UPstream::warnComm = oldWarn;
+            }
         }
-        else if (agglomeration_.processorAgglomerate())
-        {
+        //else if (agglomeration_.processorAgglomerate())
+        //{
             //// Pick a level to processor agglomerate
             //label agglomLevel = matrixLevels_.size() - 1;//1;
             //
@@ -424,7 +413,7 @@ Foam::GAMGSolver::GAMGSolver
             //
             //
             //UPstream::warnComm = oldWarn;
-        }
+        //}
     }
     else
     {
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
index 3497895241..ed49d00827 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
@@ -109,9 +109,6 @@ class GAMGSolver
         //- Direct or iteratively solve the coarsest level
         bool directSolveCoarsest_;
 
-        //- Agglomerate processors
-        bool processorAgglomerate_;
-
         //- The agglomeration
         const GAMGAgglomeration& agglomeration_;
 
@@ -119,11 +116,11 @@ class GAMGSolver
         PtrList<lduMatrix> matrixLevels_;
 
         //- Hierarchy of interfaces.
-        //  Warning: Needs to be deleted explicitly.
-        PtrList<lduInterfaceFieldPtrsList> interfaceLevels_;
-
         PtrList<PtrList<lduInterfaceField> > primitiveInterfaceLevels_;
 
+        //- Hierarchy of interfaces in lduInterfaceFieldPtrs form
+        PtrList<lduInterfaceFieldPtrsList> interfaceLevels_;
+
         //- Hierarchy of interface boundary coefficients
         PtrList<FieldField<Field, scalar> > interfaceLevelsBouCoeffs_;
 
@@ -133,15 +130,6 @@ class GAMGSolver
         //- LU decompsed coarsest matrix
         autoPtr<LUscalarMatrix> coarsestLUMatrixPtr_;
 
-        //// Processor matrix agglomeration
-
-            ////- Combined coarsest level matrix for all processors
-            //PtrList<lduMatrix> procMatrixLevels_;
-            //PtrList<lduInterfaceFieldPtrsList> procInterfaceLevels_;
-            //PtrList<PtrList<lduInterfaceField> > procPrimitiveInterfaces_;
-            //PtrList<FieldField<Field, scalar> > procInterfaceLevelsBouCoeffs_;
-            //PtrList<FieldField<Field, scalar> > procInterfaceLevelsIntCoeffs_;
-
 
     // Private Member Functions
 
@@ -184,6 +172,7 @@ class GAMGSolver
         (
             const label fineLevelIndex,
             const lduInterfacePtrsList& coarseMeshInterfaces,
+            PtrList<lduInterfaceField>& coarsePrimInterfaces,
             lduInterfaceFieldPtrsList& coarseInterfaces,
             FieldField<Field, scalar>& coarseInterfaceBouCoeffs,
             FieldField<Field, scalar>& coarseInterfaceIntCoeffs
@@ -255,7 +244,6 @@ class GAMGSolver
             scalarField& field,
             scalarField& Acf,
             const lduMatrix& A,
-            const int comm,
             const FieldField<Field, scalar>& interfaceLevelBouCoeffs,
             const lduInterfaceFieldPtrsList& interfaceLevel,
             const scalarField& source,
@@ -267,7 +255,9 @@ class GAMGSolver
         (
             PtrList<scalarField>& coarseCorrFields,
             PtrList<scalarField>& coarseSources,
-            PtrList<lduMatrix::smoother>& smoothers
+            PtrList<lduMatrix::smoother>& smoothers,
+            scalarField& scratch1,
+            scalarField& scratch2
         ) const;
 
 
@@ -280,6 +270,10 @@ class GAMGSolver
             scalarField& Apsi,
             scalarField& finestCorrection,
             scalarField& finestResidual,
+
+            scalarField& scratch1,
+            scalarField& scratch2,
+
             PtrList<scalarField>& coarseCorrFields,
             PtrList<scalarField>& coarseSources,
             const direction cmpt=0
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C
index e3754e47cc..9f40158872 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C
@@ -46,7 +46,6 @@ void Foam::GAMGSolver::agglomerateMatrix
         const label nCoarseCells = agglomeration_.nCells(fineLevelIndex);
 
 
-
         Pout<< "agglomerateMatrix : I have fine mesh:"
             << fineMatrix.mesh().lduAddr().size()
             << " at fine level:" << fineLevelIndex
@@ -82,6 +81,15 @@ void Foam::GAMGSolver::agglomerateMatrix
             interfaceLevel(fineLevelIndex);
 
         // Create coarse-level interfaces
+        primitiveInterfaceLevels_.set
+        (
+            fineLevelIndex,
+            new PtrList<lduInterfaceField>(fineInterfaces.size())
+        );
+
+        PtrList<lduInterfaceField>& coarsePrimInterfaces =
+            primitiveInterfaceLevels_[fineLevelIndex];
+
         interfaceLevels_.set
         (
             fineLevelIndex,
@@ -114,6 +122,7 @@ void Foam::GAMGSolver::agglomerateMatrix
         (
             fineLevelIndex,
             coarseMeshInterfaces,
+            coarsePrimInterfaces,
             coarseInterfaces,
             coarseInterfaceBouCoeffs,
             coarseInterfaceIntCoeffs
@@ -192,12 +201,12 @@ void Foam::GAMGSolver::agglomerateMatrix
 }
 
 
-//XXXXX
 // Agglomerate only the interface coefficients.
 void Foam::GAMGSolver::agglomerateInterfaceCoefficients
 (
     const label fineLevelIndex,
     const lduInterfacePtrsList& coarseMeshInterfaces,
+    PtrList<lduInterfaceField>& coarsePrimInterfaces,
     lduInterfaceFieldPtrsList& coarseInterfaces,
     FieldField<Field, scalar>& coarseInterfaceBouCoeffs,
     FieldField<Field, scalar>& coarseInterfaceIntCoeffs
@@ -233,7 +242,7 @@ void Foam::GAMGSolver::agglomerateInterfaceCoefficients
                     coarseMeshInterfaces[inti]
                 );
 
-            coarseInterfaces.set
+            coarsePrimInterfaces.set
             (
                 inti,
                 GAMGInterfaceField::New
@@ -242,6 +251,11 @@ void Foam::GAMGSolver::agglomerateInterfaceCoefficients
                     fineInterfaces[inti]
                 ).ptr()
             );
+            coarseInterfaces.set
+            (
+                inti,
+                &coarsePrimInterfaces[inti]
+            );
 
             const labelList& faceRestrictAddressing = patchFineToCoarse[inti];
 
@@ -364,9 +378,9 @@ void Foam::GAMGSolver::gatherMatrices
     }
     else
     {
-        Pout<< "GAMGSolver::gatherMatrices :"
-            << " sending from:" << Pstream::myProcNo(meshComm)
-            << " to master:" << procIDs[0] << endl;
+        //Pout<< "GAMGSolver::gatherMatrices :"
+        //    << " sending from:" << Pstream::myProcNo(meshComm)
+        //    << " to master:" << procIDs[0] << endl;
 
         // Send to master
 
@@ -503,7 +517,6 @@ void Foam::GAMGSolver::procAgglomerateMatrix
 
 Pout<< "Agglomerating onto mesh:" << allMesh.info() << endl;
 
-
         allMatrixPtr.reset(new lduMatrix(allMesh));
         lduMatrix& allMatrix = allMatrixPtr();
 
@@ -666,8 +679,9 @@ Pout<< "Agglomerating onto mesh:" << allMesh.info() << endl;
 
                     const labelList& map = boundaryFaceMap[procI][procIntI];
 
-Pout<< "    from proc:" << procI << " interface:" << procIntI
-<< " mapped to faces:" << map << endl;
+                    //Pout<< "    from proc:" << procI
+                    //<< " interface:" << procIntI
+                    //<< " mapped to faces:" << map << endl;
 
                     const scalarField& procBou = procBouCoeffs[procIntI];
                     const scalarField& procInt = procIntCoeffs[procIntI];
@@ -697,6 +711,7 @@ Pout<< "    from proc:" << procI << " interface:" << procIntI
                     const scalarField& procBou = procBouCoeffs[procIntI];
                     const scalarField& procInt = procIntCoeffs[procIntI];
 
+
                     forAll(map, i)
                     {
                         if (map[i] >= 0)
@@ -779,22 +794,21 @@ Pout<< "    from proc:" << procI << " interface:" << procIntI
                     << allInterfaces[intI].interface().
                         faceCells().size()
                     << endl;
-                const scalarField& bouCoeffs = allInterfaceBouCoeffs[intI];
-                const scalarField& intCoeffs = allInterfaceIntCoeffs[intI];
 
-                forAll(bouCoeffs, faceI)
-                {
-                    Pout<< "        " << faceI
-                        << "\tbou:" << bouCoeffs[faceI]
-                        << "\tint:" << intCoeffs[faceI]
-                        << endl;
-                }
+                //const scalarField& bouCoeffs = allInterfaceBouCoeffs[intI];
+                //const scalarField& intCoeffs = allInterfaceIntCoeffs[intI];
+                //forAll(bouCoeffs, faceI)
+                //{
+                //    Pout<< "        " << faceI
+                //        << "\tbou:" << bouCoeffs[faceI]
+                //        << "\tint:" << intCoeffs[faceI]
+                //        << endl;
+                //}
             }
         }
     }
     UPstream::warnComm = oldWarn;
 }
-//XXXX
 
 
 void Foam::GAMGSolver::procAgglomerateMatrix
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C
index d7f73c1266..b81455ff4b 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C
@@ -33,7 +33,6 @@ void Foam::GAMGSolver::scale
     scalarField& field,
     scalarField& Acf,
     const lduMatrix& A,
-    const int comm,
     const FieldField<Field, scalar>& interfaceLevelBouCoeffs,
     const lduInterfaceFieldPtrsList& interfaceLevel,
     const scalarField& source,
@@ -59,12 +58,7 @@ void Foam::GAMGSolver::scale
     }
 
     vector2D scalingVector(scalingFactorNum, scalingFactorDenom);
-    //A.mesh().reduce
-    //(
-    //    scalingVector,
-    //    sumOp<vector2D>()
-    //);
-    Foam::reduce(scalingVector, sumOp<vector2D>(), Pstream::msgType(), comm);
+    A.mesh().reduce(scalingVector, sumOp<vector2D>());
 
     scalar sf = scalingVector.x()/stabilise(scalingVector.y(), VSMALL);
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
index 61022c353d..ae2cfa21ad 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
@@ -27,7 +27,6 @@ License
 #include "ICCG.H"
 #include "BICCG.H"
 #include "SubField.H"
-#include "globalIndex.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -81,8 +80,20 @@ Foam::solverPerformance Foam::GAMGSolver::solve
         // Create the smoothers for all levels
         PtrList<lduMatrix::smoother> smoothers;
 
+        // Scratch fields if processor-agglomerated coarse level meshes
+        // are bigger than original. Usually not needed
+        scalarField scratch1;
+        scalarField scratch2;
+
         // Initialise the above data structures
-        initVcycle(coarseCorrFields, coarseSources, smoothers);
+        initVcycle
+        (
+            coarseCorrFields,
+            coarseSources,
+            smoothers,
+            scratch1,
+            scratch2
+        );
 
         do
         {
@@ -94,11 +105,22 @@ Foam::solverPerformance Foam::GAMGSolver::solve
                 Apsi,
                 finestCorrection,
                 finestResidual,
+
+                (scratch1.size() ? scratch1 : Apsi),
+                (scratch2.size() ? scratch2 : finestCorrection),
+
                 coarseCorrFields,
                 coarseSources,
                 cmpt
             );
 
+
+            //Pout<< "finestCorrection:" << finestCorrection << endl;
+            //Pout<< "finestResidual:" << finestResidual << endl;
+            //Pout<< "psi:" << psi << endl;
+            //Pout<< "Apsi:" << Apsi << endl;
+
+
             // Calculate finest level residual field
             matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
             finestResidual = source;
@@ -133,6 +155,10 @@ void Foam::GAMGSolver::Vcycle
     scalarField& Apsi,
     scalarField& finestCorrection,
     scalarField& finestResidual,
+
+    scalarField& scratch1,
+    scalarField& scratch2,
+
     PtrList<scalarField>& coarseCorrFields,
     PtrList<scalarField>& coarseSources,
     const direction cmpt
@@ -154,8 +180,6 @@ void Foam::GAMGSolver::Vcycle
     // Residual restriction (going to coarser levels)
     for (label leveli = 0; leveli < coarsestLevel; leveli++)
     {
-        Pout<< "Restriction for level:" << leveli << endl;
-
         if (coarseSources.set(leveli + 1))
         {
             // If the optional pre-smoothing sweeps are selected
@@ -178,16 +202,15 @@ void Foam::GAMGSolver::Vcycle
 
                 scalarField::subField ACf
                 (
-                    Apsi,
+                    scratch1,
                     coarseCorrFields[leveli].size()
                 );
+                //scalarField ACf(coarseCorrFields[leveli].size(), VGREAT);
 
                 // Scale coarse-grid correction field
                 // but not on the coarsest level because it evaluates to 1
                 if (scaleCorrection_ && leveli < coarsestLevel - 1)
                 {
-                    int comm = matrixLevels_[leveli].mesh().comm();
-
                     scale
                     (
                         coarseCorrFields[leveli],
@@ -195,8 +218,9 @@ void Foam::GAMGSolver::Vcycle
                         (
                             ACf.operator const scalarField&()
                         ),
+                        //ACf,
+
                         matrixLevels_[leveli],
-                        comm,
                         interfaceLevelsBouCoeffs_[leveli],
                         interfaceLevels_[leveli],
                         coarseSources[leveli],
@@ -211,6 +235,7 @@ void Foam::GAMGSolver::Vcycle
                     (
                         ACf.operator const scalarField&()
                     ),
+                    //ACf,
                     coarseCorrFields[leveli],
                     interfaceLevelsBouCoeffs_[leveli],
                     interfaceLevels_[leveli],
@@ -240,7 +265,6 @@ void Foam::GAMGSolver::Vcycle
     // Solve Coarsest level with either an iterative or direct solver
     if (coarseCorrFields.set(coarsestLevel))
     {
-        Pout<< "Coarsest solve for level:" << coarsestLevel << endl;
         solveCoarsestLevel
         (
             coarseCorrFields[coarsestLevel],
@@ -260,28 +284,27 @@ void Foam::GAMGSolver::Vcycle
 
     for (label leveli = coarsestLevel - 1; leveli >= 0; leveli--)
     {
-        Pout<< "Smoothing and prolongation for level:" << leveli << endl;
-
         if (coarseCorrFields.set(leveli))
         {
-            Pout<< "Prolonging from " << leveli + 1 << " up to "
-                << leveli << endl;
-            //// Create a field for the pre-smoothed correction field
-            //// as a sub-field of the finestCorrection which is not
-            //// currently being used
-            //scalarField::subField preSmoothedCoarseCorrField
+            // Create a field for the pre-smoothed correction field
+            // as a sub-field of the finestCorrection which is not
+            // currently being used
+            scalarField::subField preSmoothedCoarseCorrField
+            (
+                scratch2,
+                coarseCorrFields[leveli].size()
+            );
+            //scalarField preSmoothedCoarseCorrField
             //(
-            //    finestCorrection,
-            //    coarseCorrFields[leveli].size()
+            //    coarseCorrFields[leveli].size(),
+            //    VGREAT
             //);
-            scalarField preSmoothedCoarseCorrField;
 
             // Only store the preSmoothedCoarseCorrField if pre-smoothing is
             // used
             if (nPreSweeps_)
             {
-                //preSmoothedCoarseCorrField.assign(coarseCorrFields[leveli]);
-                preSmoothedCoarseCorrField = coarseCorrFields[leveli];
+                preSmoothedCoarseCorrField.assign(coarseCorrFields[leveli]);
             }
 
             agglomeration_.prolongField
@@ -296,22 +319,24 @@ void Foam::GAMGSolver::Vcycle
                 true
             );
 
-            Pout<< "Doing stuff at level " << leveli << endl;
 
-            //// Create A.psi for this coarse level as a sub-field of Apsi
-            //scalarField::subField ACf
+            // Create A.psi for this coarse level as a sub-field of Apsi
+            scalarField::subField ACf
+            (
+               scratch1,
+                coarseCorrFields[leveli].size()
+            );
+            scalarField& ACfRef =
+                const_cast<scalarField&>(ACf.operator const scalarField&());
+            //scalarField ACfRef
             //(
-            //    Apsi,
-            //    coarseCorrFields[leveli].size()
+            //    coarseCorrFields[leveli].size(),
+            //    VGREAT
             //);
-            //scalarField& ACfRef =
-            //    const_cast<scalarField&>(ACf.operator const scalarField&());
-            scalarField ACfRef(coarseCorrFields[leveli].size());
 
 
             if (interpolateCorrection_)
             {
-                Pout<< "doing interpolate." << endl;
                 interpolate
                 (
                     coarseCorrFields[leveli],
@@ -322,35 +347,22 @@ void Foam::GAMGSolver::Vcycle
                     coarseSources[leveli],
                     cmpt
                 );
-                Pout<< "done interpolate." << endl;
             }
 
             // Scale coarse-grid correction field
             // but not on the coarsest level because it evaluates to 1
             if (scaleCorrection_ && leveli < coarsestLevel - 1)
             {
-                //int comm =
-                //(
-                //    matrixLevels_.set(leveli+1)
-                //  ? matrixLevels_[leveli+1].mesh().comm()
-                //  : matrixLevels_[leveli].mesh().comm()
-                //);
-                int comm = matrixLevels_[leveli].mesh().comm();
-
-
-                Pout<< "doing scale with comm:" << comm << endl;
                 scale
                 (
                     coarseCorrFields[leveli],
                     ACfRef,
                     matrixLevels_[leveli],
-                    comm,
                     interfaceLevelsBouCoeffs_[leveli],
                     interfaceLevels_[leveli],
                     coarseSources[leveli],
                     cmpt
                 );
-                Pout<< "done scale with comm:" << comm << endl;
             }
 
             // Only add the preSmoothedCoarseCorrField if pre-smoothing is
@@ -360,7 +372,6 @@ void Foam::GAMGSolver::Vcycle
                 coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
             }
 
-            Pout<< "doing smooth." << endl;
             smoothers[leveli + 1].smooth
             (
                 coarseCorrFields[leveli],
@@ -372,24 +383,20 @@ void Foam::GAMGSolver::Vcycle
                     maxPostSweeps_
                 )
             );
-            Pout<< "done smooth." << endl;
         }
     }
 
     // Prolong the finest level correction
-    Pout<< "Doing Prolong to finest level" << endl;
     agglomeration_.prolongField
     (
         finestCorrection,
         coarseCorrFields[0],
         0,
-        true            //false               // no proc agglomeration for now
+        true
     );
-    Pout<< "Done Prolong to finest level" << endl;
 
     if (interpolateCorrection_)
     {
-        Pout<< "doing interpolate on finest level" << endl;
         interpolate
         (
             finestCorrection,
@@ -400,26 +407,21 @@ void Foam::GAMGSolver::Vcycle
             finestResidual,
             cmpt
         );
-        Pout<< "done interpolate on finest level" << endl;
     }
 
     if (scaleCorrection_)
     {
         // Scale the finest level correction
-        int comm = matrix_.mesh().comm();
-        Pout<< "doing scale on finest level with comm:" << comm << endl;
         scale
         (
             finestCorrection,
             Apsi,
             matrix_,
-            comm,
             interfaceBouCoeffs_,
             interfaces_,
             finestResidual,
             cmpt
         );
-        Pout<< "done scale on finest level with comm:" << comm << endl;
     }
 
     forAll(psi, i)
@@ -427,7 +429,6 @@ void Foam::GAMGSolver::Vcycle
         psi[i] += finestCorrection[i];
     }
 
-    Pout<< "Doing smooth on finest level" << endl;
     smoothers[0].smooth
     (
         psi,
@@ -435,7 +436,6 @@ void Foam::GAMGSolver::Vcycle
         cmpt,
         nFinestSweeps_
     );
-    Pout<< "Done smooth on finest level" << endl;
 }
 
 
@@ -443,9 +443,13 @@ void Foam::GAMGSolver::initVcycle
 (
     PtrList<scalarField>& coarseCorrFields,
     PtrList<scalarField>& coarseSources,
-    PtrList<lduMatrix::smoother>& smoothers
+    PtrList<lduMatrix::smoother>& smoothers,
+    scalarField& scratch1,
+    scalarField& scratch2
 ) const
 {
+    label maxSize = matrix_.diag().size();
+
     coarseCorrFields.setSize(matrixLevels_.size());
     coarseSources.setSize(matrixLevels_.size());
     smoothers.setSize(matrixLevels_.size() + 1);
@@ -471,15 +475,7 @@ void Foam::GAMGSolver::initVcycle
         {
             label nCoarseCells = agglomeration_.nCells(leveli);
 
-            Pout<< "initVCucle level:" << leveli << " nCoarseCells:"
-                << nCoarseCells << endl;
-
             coarseSources.set(leveli, new scalarField(nCoarseCells));
-
-            //if (!matrixLevels_.set(leveli))
-            //{
-            //    coarseCorrFields.set(leveli, new scalarField(nCoarseCells));
-            //}
         }
 
         if (matrixLevels_.set(leveli))
@@ -487,12 +483,10 @@ void Foam::GAMGSolver::initVcycle
             const lduMatrix& mat = matrixLevels_[leveli];
 
             label nCoarseCells = mat.diag().size();
-            Pout<< "initVCucle level:" << leveli << " matrix size:"
-                << nCoarseCells << endl;
+
+            maxSize = max(maxSize, nCoarseCells);
 
             coarseCorrFields.set(leveli, new scalarField(nCoarseCells));
-            //coarseCorrFields.set(leveli, new scalarField(nCoarseCells));
-            //coarseSources.set(leveli, new scalarField(nCoarseCells));
 
             smoothers.set
             (
@@ -509,6 +503,13 @@ void Foam::GAMGSolver::initVcycle
             );
         }
     }
+
+    if (maxSize > matrix_.diag().size())
+    {
+        // Allocate some scratch storage
+        scratch1.setSize(maxSize);
+        scratch2.setSize(maxSize);
+    }
 }
 
 
@@ -520,10 +521,6 @@ void Foam::GAMGSolver::solveCoarsestLevel
 {
     const label coarsestLevel = matrixLevels_.size() - 1;
 
-    Pout<< "solveCoarsestLevel :"
-        << " coarsestLevel:" << coarsestLevel << endl;
-
-
     label coarseComm = matrixLevels_[coarsestLevel].mesh().comm();
     label oldWarn = UPstream::warnComm;
     UPstream::warnComm = coarseComm;
@@ -682,9 +679,6 @@ void Foam::GAMGSolver::solveCoarsestLevel
         {
             coarseSolverPerf.print(Info(coarseComm));
         }
-
-        Pout<< "GC: coarsestSource   :" << coarsestSource << endl;
-        Pout<< "GC: coarsestCorrField:" << coarsestCorrField << endl;
     }
 
     UPstream::warnComm = oldWarn;
diff --git a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C
index 72cd0a30e4..f33504a8b0 100644
--- a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C
+++ b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C
@@ -67,15 +67,13 @@ void Foam::lduPrimitiveMesh::checkUpperTriangular
     {
         if (u[faceI] < l[faceI])
         {
-            //FatalErrorIn
-            WarningIn
+            FatalErrorIn
             (
                 "checkUpperTriangular"
                 "(const label, const labelUList&, const labelUList&)"
             )   << "Reversed face. Problem at face " << faceI
                 << " l:" << l[faceI] << " u:" << u[faceI]
-                //<< abort(FatalError);
-                << endl;
+                << abort(FatalError);
         }
         if (l[faceI] < 0 || u[faceI] < 0 || u[faceI] >= size)
         {
@@ -93,8 +91,7 @@ void Foam::lduPrimitiveMesh::checkUpperTriangular
     {
         if (l[faceI-1] > l[faceI])
         {
-            //FatalErrorIn
-            WarningIn
+            FatalErrorIn
             (
                 "checkUpperTriangular"
                 "(const label, const labelUList&, const labelUList&)"
@@ -102,16 +99,14 @@ void Foam::lduPrimitiveMesh::checkUpperTriangular
                 << " Problem at face " << faceI
                 << " l:" << l[faceI] << " u:" << u[faceI]
                 << " previous l:" << l[faceI-1]
-                //<< abort(FatalError);
-                << endl;
+                << abort(FatalError);
         }
         else if (l[faceI-1] == l[faceI])
         {
             // Same cell.
             if (u[faceI-1] > u[faceI])
             {
-                //FatalErrorIn
-                WarningIn
+                FatalErrorIn
                 (
                     "checkUpperTriangular"
                     "(const label, const labelUList&, const labelUList&)"
@@ -119,8 +114,7 @@ void Foam::lduPrimitiveMesh::checkUpperTriangular
                     << " Problem at face " << faceI
                     << " l:" << l[faceI] << " u:" << u[faceI]
                     << " previous u:" << u[faceI-1]
-                    //<< abort(FatalError);
-                    << endl;
+                    << abort(FatalError);
             }
         }
     }
@@ -139,6 +133,79 @@ Foam::label Foam::lduPrimitiveMesh::size(const PtrList<lduMesh>& meshes)
 }
 
 
+Foam::labelList Foam::lduPrimitiveMesh::upperTriOrder
+(
+    const label nCells,
+    const labelUList& lower,
+    const labelUList& upper
+)
+{
+    labelList nNbrs(nCells, 0);
+
+    // Count number of upper neighbours
+    forAll(lower, faceI)
+    {
+        if (upper[faceI] < lower[faceI])
+        {
+            FatalErrorIn("lduPrimitiveMesh::upperTriOrder(..)")
+                << "Problem at face:" << faceI
+                << " lower:" << lower[faceI]
+                << " upper:" << upper[faceI]
+                << exit(FatalError);
+        }
+        nNbrs[lower[faceI]]++;
+    }
+
+    // Construct cell-upper cell addressing
+    labelList offsets(nCells+1);
+    offsets[0] = 0;
+    forAll(nNbrs, cellI)
+    {
+        offsets[cellI+1] = offsets[cellI]+nNbrs[cellI];
+    }
+
+    nNbrs = offsets;
+
+    labelList cellToFaces(offsets.last());
+    forAll(upper, faceI)
+    {
+        label cellI = lower[faceI];
+        cellToFaces[nNbrs[cellI]++] = faceI;
+    }
+
+    // Sort
+
+    labelList oldToNew(lower.size());
+
+    labelList order;
+    labelList nbr;
+
+    label newFaceI = 0;
+
+    for (label cellI = 0; cellI < nCells; cellI++)
+    {
+        label startOfCell = offsets[cellI];
+        label nNbr = offsets[cellI+1] - startOfCell;
+
+        nbr.setSize(nNbr);
+        order.setSize(nNbr);
+        forAll(order, i)
+        {
+            nbr[i] = upper[cellToFaces[offsets[cellI]+i]];
+        }
+        sortedOrder(nbr, order);
+
+        forAll(order, i)
+        {
+            label index = order[i];
+            oldToNew[cellToFaces[startOfCell + index]] = newFaceI++;
+        }
+    }
+
+    return oldToNew;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::lduPrimitiveMesh::lduPrimitiveMesh
@@ -159,36 +226,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     interfaces_(interfaces),
     patchSchedule_(ps),
     comm_(comm)
-{
-    //Pout<< "lduPrimitiveMesh :"
-    //    << " nCells:" << nCells
-    //    << " l:" << lowerAddr_.size()
-    //    << " u:" << upperAddr_.size()
-    //    << " pa:" << patchAddr_.size()
-    //    << " interfaces:" << interfaces_.size()
-    //    << " comm:" << comm_
-    //    << endl;
-    //forAll(interfaces_, i)
-    //{
-    //    if (interfaces_.set(i))
-    //    {
-    //        if (isA<processorLduInterface>(interfaces_[i]))
-    //        {
-    //            const processorLduInterface& pi = refCast
-    //            <
-    //                const processorLduInterface
-    //            >(interfaces_[i]);
-    //
-    //            Pout<< "    patch:" << i
-    //                << " size:" << patchAddr_[i].size()
-    //                << " myProcNo:" << pi.myProcNo()
-    //                << " neighbProcNo:" << pi.neighbProcNo()
-    //                << " comm:" << pi.comm()
-    //                << endl;
-    //        }
-    //    }
-    //}
-}
+{}
 
 
 Foam::lduPrimitiveMesh::lduPrimitiveMesh
@@ -210,95 +248,40 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     interfaces_(interfaces, reUse),
     patchSchedule_(ps),
     comm_(comm)
+{}
+
+
+Foam::lduPrimitiveMesh::lduPrimitiveMesh
+(
+    const label nCells,
+    labelList& l,
+    labelList& u,
+    labelListList& pa,
+    const Xfer<PtrList<const lduInterface> >& primitiveInterfaces,
+    const lduSchedule& ps,
+    const label comm
+)
+:
+    lduAddressing(nCells),
+    lowerAddr_(l, true),
+    upperAddr_(u, true),
+    patchAddr_(pa, true),
+    primitiveInterfaces_(primitiveInterfaces),
+    patchSchedule_(ps),
+    comm_(comm)
 {
-    //Pout<< "lduPrimitiveMesh :"
-    //    << " nCells:" << nCells
-    //    << " l:" << lowerAddr_.size()
-    //    << " u:" << upperAddr_.size()
-    //    << " pa:" << patchAddr_.size()
-    //    << " interfaces:" << interfaces_.size()
-    //    << " comm:" << comm_
-    //    << endl;
-    //forAll(interfaces_, i)
-    //{
-    //    if (interfaces_.set(i))
-    //    {
-    //        if (isA<processorLduInterface>(interfaces_[i]))
-    //        {
-    //            const processorLduInterface& pi = refCast
-    //            <
-    //                const processorLduInterface
-    //            >(interfaces_[i]);
-    //
-    //            Pout<< "    patch:" << i
-    //                << " size:" << patchAddr_[i].size()
-    //                << " myProcNo:" << pi.myProcNo()
-    //                << " neighbProcNo:" << pi.neighbProcNo()
-    //                << " comm:" << pi.comm()
-    //                << endl;
-    //        }
-    //    }
-    //}
+    // Create interfaces
+    interfaces_.setSize(primitiveInterfaces_.size());
+    forAll(primitiveInterfaces_, i)
+    {
+        if (primitiveInterfaces_.set(i))
+        {
+            interfaces_.set(i, &primitiveInterfaces_[i]);
+        }
+    }
 }
 
 
-//Foam::lduPrimitiveMesh::lduPrimitiveMesh
-//(
-//    const label nCells,
-//    labelList& l,
-//    labelList& u,
-//    labelListList& pa,
-//    PtrList<const lduInterface>& primitiveInterfaces,
-//    const lduSchedule& ps,
-//    const label comm,
-//    bool reUse
-//)
-//:
-//    lduAddressing(nCells),
-//    lowerAddr_(l, reUse),
-//    upperAddr_(u, reUse),
-//    patchAddr_(pa, reUse),
-//    primitiveInterfaces_(primitiveInterfaces, reUse),
-//    patchSchedule_(ps),
-//    comm_(comm)
-//{
-//    interfaces_.setSize(primitiveInterfaces_.size());
-//    forAll(primitiveInterfaces, intI)
-//    {
-//        interfaces_.set(intI, &primitiveInterfaces_[intI]);
-//    }
-//
-//    //Pout<< "lduPrimitiveMesh :"
-//    //    << " nCells:" << nCells
-//    //    << " l:" << lowerAddr_.size()
-//    //    << " u:" << upperAddr_.size()
-//    //    << " pa:" << patchAddr_.size()
-//    //    << " interfaces:" << interfaces_.size()
-//    //    << " comm:" << comm_
-//    //    << endl;
-//    //forAll(interfaces_, i)
-//    //{
-//    //    if (interfaces_.set(i))
-//    //    {
-//    //        if (isA<processorLduInterface>(interfaces_[i]))
-//    //        {
-//    //            const processorLduInterface& pi = refCast
-//    //            <
-//    //                const processorLduInterface
-//    //            >(interfaces_[i]);
-//    //
-//    //            Pout<< "    patch:" << i
-//    //                << " size:" << patchAddr_[i].size()
-//    //                << " myProcNo:" << pi.myProcNo()
-//    //                << " neighbProcNo:" << pi.neighbProcNo()
-//    //                << " comm:" << pi.comm()
-//    //                << endl;
-//    //        }
-//    //    }
-//    //}
-//}
-
-
 Foam::lduPrimitiveMesh::lduPrimitiveMesh
 (
     const label comm,
@@ -343,10 +326,14 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     const label nMeshes = otherMeshes.size()+1;
 
     const label myAgglom = procAgglomMap[UPstream::myProcNo(currentComm)];
-    Pout<< "I am " << UPstream::myProcNo(currentComm)
-        << " agglomerating into " << myAgglom
-        << " as are " << findIndices(procAgglomMap, myAgglom)
-        << endl;
+
+    if (lduPrimitiveMesh::debug)
+    {
+        Pout<< "I am " << UPstream::myProcNo(currentComm)
+            << " agglomerating into " << myAgglom
+            << " as are " << findIndices(procAgglomMap, myAgglom)
+            << endl;
+    }
 
 
     forAll(procIDs, i)
@@ -448,11 +435,14 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                     else if (agglom0 == myAgglom && agglom1 == myAgglom)
                     {
                         // Merged interface
-                        Pout<< "merged interface: myProcNo:"
-                            << pldui.myProcNo()
-                            << " nbr:" << pldui.neighbProcNo()
-                            << " size:" << ldui.faceCells().size()
-                            << endl;
+                        if (debug)
+                        {
+                            Pout<< "merged interface: myProcNo:"
+                                << pldui.myProcNo()
+                                << " nbr:" << pldui.neighbProcNo()
+                                << " size:" << ldui.faceCells().size()
+                                << endl;
+                        }
 
                         label nbrProcMeshI = findIndex
                         (
@@ -484,11 +474,14 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                     }
                     else
                     {
-                        Pout<< "external interface: myProcNo:"
-                            << pldui.myProcNo()
-                            << " nbr:" << pldui.neighbProcNo()
-                            << " size:" << ldui.faceCells().size()
-                            << endl;
+                        if (debug)
+                        {
+                            Pout<< "external interface: myProcNo:"
+                                << pldui.myProcNo()
+                                << " nbr:" << pldui.neighbProcNo()
+                                << " size:" << ldui.faceCells().size()
+                                << endl;
+                        }
 
                         EdgeMap<labelPairList>::iterator iter =
                             unmergedMap.find(procEdge);
@@ -524,7 +517,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
 
 
 
-    //if (debug)
+    if (debug)
     {
         Pout<< "Remaining interfaces:" << endl;
         forAllConstIter(EdgeMap<labelPairList>, unmergedMap, iter)
@@ -553,7 +546,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
             Pout<< endl;
         }
     }
-    //if (debug)
+    if (debug)
     {
         Pout<< "Merged interfaces:" << endl;
         forAllConstIter(EdgeMap<labelPairList>, mergedMap, iter)
@@ -632,6 +625,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
             allFaceI++;
         }
 
+
         // Add merged interfaces
         const lduInterfacePtrsList interfaces = procMesh.interfaces();
 
@@ -670,7 +664,26 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                             label nbrIntI = -1;
                             forAll(elems, i)
                             {
-                                if (elems[i][0] == nbrProcMeshI)
+                                label procI = elems[i][0];
+                                label interfaceI = elems[i][1];
+                                const lduInterfacePtrsList interfaces =
+                                    mesh
+                                    (
+                                        myMesh,
+                                        otherMeshes,
+                                        procI
+                                    ).interfaces();
+                                const processorLduInterface& pldui =
+                                    refCast<const processorLduInterface>
+                                    (
+                                        interfaces[interfaceI]
+                                    );
+
+                                if
+                                (
+                                    elems[i][0] == nbrProcMeshI
+                                 && pldui.neighbProcNo() == procMeshI
+                                )
                                 {
                                     nbrIntI = elems[i][1];
                                     break;
@@ -700,6 +713,17 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                             const labelUList& nbrFaceCells =
                                 nbrInterfaces[nbrIntI].faceCells();
 
+                            if (faceCells.size() != nbrFaceCells.size())
+                            {
+                                FatalErrorIn
+                                (
+                                    "lduPrimitiveMesh::lduPrimitiveMesh(..)"
+                                )   << "faceCells:" << faceCells
+                                    << " nbrFaceCells:" << nbrFaceCells
+                                    << abort(FatalError);
+                            }
+
+
                             labelList& bfMap =
                                 boundaryFaceMap[procMeshI][intI];
                             labelList& nbrBfMap =
@@ -726,6 +750,67 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     }
 
 
+    // Sort upper-tri order
+    {
+        labelList oldToNew
+        (
+            upperTriOrder
+            (
+                cellOffsets.last(), //nCells
+                lowerAddr_,
+                upperAddr_
+            )
+        );
+
+        forAll(faceMap, procMeshI)
+        {
+            labelList& map = faceMap[procMeshI];
+            forAll(map, i)
+            {
+                if (map[i] >= 0)
+                {
+                    map[i] = oldToNew[map[i]];
+                }
+                else
+                {
+                    label allFaceI = -map[i]-1;
+                    map[i] = -oldToNew[allFaceI]-1;
+                }
+            }
+        }
+
+
+        inplaceReorder(oldToNew, lowerAddr_);
+        inplaceReorder(oldToNew, upperAddr_);
+
+        forAll(boundaryFaceMap, procI)
+        {
+            const labelList& bMap = boundaryMap[procI];
+            forAll(bMap, intI)
+            {
+                if (bMap[intI] == -1)
+                {
+                    // Merged interface
+                    labelList& bfMap = boundaryFaceMap[procI][intI];
+
+                    forAll(bfMap, i)
+                    {
+                        if (bfMap[i] >= 0)
+                        {
+                            bfMap[i] = oldToNew[bfMap[i]];
+                        }
+                        else
+                        {
+                            label allFaceI = -bfMap[i]-1;
+                            bfMap[i] = (-oldToNew[allFaceI]-1);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+
     // Kept interfaces
     // ~~~~~~~~~~~~~~~
 
@@ -736,7 +821,6 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
 
     forAllConstIter(EdgeMap<labelPairList>, unmergedMap, iter)
     {
-        Pout<< "agglom procEdge:" << iter.key() << endl;
         const labelPairList& elems = iter();
 
         // Sort processors in increasing order
@@ -758,10 +842,10 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                 procMeshI
             ).interfaces();
 
-            Pout<< "    adding interface:" << " procMeshI:" << procMeshI
-                << " interface:" << interfaceI
-                << " size:" << interfaces[interfaceI].faceCells().size()
-                << endl;
+            //Pout<< "    adding interface:" << " procMeshI:" << procMeshI
+            //    << " interface:" << interfaceI
+            //    << " size:" << interfaces[interfaceI].faceCells().size()
+            //    << endl;
 
             n += interfaces[interfaceI].faceCells().size();
         }
@@ -850,13 +934,16 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
         );
         interfaces_.set(allInterfaceI, &primitiveInterfaces_[allInterfaceI]);
 
-        Pout<< "Created " << interfaces_[allInterfaceI].type()
-            << " interface at " << allInterfaceI
-            << " comm:" << comm_
-            << " myProcNo:" << myAgglom
-            << " neighbProcNo:" << neighbProcNo
-            << " nFaces:" << allFaceCells.size()
-            << endl;
+        if (debug)
+        {
+            Pout<< "Created " << interfaces_[allInterfaceI].type()
+                << " interface at " << allInterfaceI
+                << " comm:" << comm_
+                << " myProcNo:" << myAgglom
+                << " neighbProcNo:" << neighbProcNo
+                << " nFaces:" << allFaceCells.size()
+                << endl;
+        }
 
 
         allInterfaceI++;
@@ -992,12 +1079,6 @@ void Foam::lduPrimitiveMesh::gather
             comm
         );
 
-        //Pout<< "sent nCells:" << addressing.size()
-        //    << "sent lowerAddr:" << addressing.lowerAddr().size()
-        //    << "sent upperAddr:" << addressing.upperAddr()
-        //    << "sent validInterface:" << validInterface
-        //    << endl;
-
         toMaster
             << addressing.size()
             << addressing.lowerAddr()
diff --git a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H
index 1bd01e9fa6..f78dad68d8 100644
--- a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H
+++ b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H
@@ -81,6 +81,13 @@ class lduPrimitiveMesh
         //- Get size of all meshes
         static label size(const PtrList<lduMesh>&);
 
+        static labelList upperTriOrder
+        (
+            const label nCells,
+            const labelUList& lower,
+            const labelUList& upper
+        );
+
         //- Check if in upper-triangular ordering
         static void checkUpperTriangular
         (
@@ -129,6 +136,18 @@ public:
             bool reUse
         );
 
+        //- Construct from components and re-use storage.
+        lduPrimitiveMesh
+        (
+            const label nCells,
+            labelList& l,
+            labelList& u,
+            labelListList& pa,
+            const Xfer<PtrList<const lduInterface> >& primitiveInterfaces,
+            const lduSchedule& ps,
+            const label comm
+        );
+
         //- Construct by combining multiple meshes. The meshes come from
         //  processors procIDs:
         //  procIDs[0] : local processor (myMesh)