From fee11f0fa56b98d1bee9e52d9084452309830e2e Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 22 Feb 2013 17:26:12 +0000
Subject: [PATCH] ENH: laplacianFoam: back to solve on single processor

---
 .../laplacianFoam.C                           | 1186 ++++-------------
 1 file changed, 278 insertions(+), 908 deletions(-)

diff --git a/applications/test/laplacianFoam-communicators/laplacianFoam.C b/applications/test/laplacianFoam-communicators/laplacianFoam.C
index 3f7bf4cf6d..8a54cb256c 100644
--- a/applications/test/laplacianFoam-communicators/laplacianFoam.C
+++ b/applications/test/laplacianFoam-communicators/laplacianFoam.C
@@ -34,6 +34,8 @@ Description
 #include "globalIndex.H"
 #include "lduPrimitiveMesh.H"
 #include "processorGAMGInterface.H"
+#include "GAMGInterfaceField.H"
+#include "processorLduInterfaceField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -126,7 +128,7 @@ void print(const string& msg, const lduMesh& mesh)
         }
     }
 
-    Pout<< "Patches:" << nl;
+    Pout<< "    Patches:" << nl;
     forAll(interfaces, i)
     {
         if (interfaces.set(i))
@@ -135,7 +137,7 @@ void print(const string& msg, const lduMesh& mesh)
             {
                 const processorLduInterface& pldui =
                     refCast<const processorLduInterface>(interfaces[i]);
-                Pout<< "    " << i
+                Pout<< "        " << i
                     << " me:" << pldui.myProcNo()
                     << " nbr:" << pldui.neighbProcNo()
                     << " comm:" << pldui.comm()
@@ -144,7 +146,7 @@ void print(const string& msg, const lduMesh& mesh)
             }
 
             {
-                Pout<< "    " << i << " addressing:" << nl;
+                Pout<< "        " << i << " addressing:" << nl;
                 const labelUList& faceCells = interfaces[i].faceCells();
                 forAll(faceCells, i)
                 {
@@ -156,7 +158,6 @@ void print(const string& msg, const lduMesh& mesh)
 }
 
 
-
 template<class ProcPatch>
 lduSchedule nonBlockingSchedule
 (
@@ -203,503 +204,6 @@ lduSchedule nonBlockingSchedule
 }
 
 
-
-// Say combining procs 1,13,9,4.
-// - cells get added in order 1,4,9,13
-// - internal faces get added in order of processor and then order of
-//   neighbouring processor:
-//     - internal faces of 1
-//     - faces between 1 and 4 (keep orientation)
-//     - faces between 1 and 9
-//     - faces between 1 and 13
-//     - interfaces between 1 and other processors
-//     and then
-//     - internal faces of 4
-//     - faces between 4 and 1 (skip; already added)
-//     - faces between 4 and 9
-//     - faces between 4 and 13
-//     - interfaces between 4 and other processors
-//     etc.
-// - this still does not guarantee that upper-triangular order is preserved:
-// A cell has two processor faces. These now become internal faces and if
-// the cell numbering on the other side is not in the same order as the internal
-// faces we have upper-triangular.
-autoPtr<lduPrimitiveMesh> combineMeshes
-(
-    const label newComm,
-    const labelList& procIDs,
-    const PtrList<lduMesh>& procMeshes,
-
-    labelList& cellOffsets,     // per mesh the starting cell
-    labelList& faceOffsets,     // per mesh the starting face
-    labelList& interfaceOffsets // per mesh,per interface the starting face
-)
-{
-    // Sanity check.
-    for (label i = 1; i < procIDs.size(); i++)
-    {
-        if (procIDs[i] <= procIDs[i-1])
-        {
-            FatalErrorIn
-            (
-                "combineMeshes(const labelList&, const PtrList<lduMesh>&)"
-            )   << "Processor " << procIDs[i] << " at index " << i
-                << " should be higher numbered than its predecessor "
-                << procIDs[i-1]
-                << exit(FatalError);
-        }
-    }
-
-
-    label currentComm = procMeshes[0].comm();
-
-
-    // Cells get added in order.
-    cellOffsets.setSize(procMeshes.size()+1);
-    cellOffsets[0] = 0;
-    forAll(procMeshes, procMeshI)
-    {
-        const lduMesh& mesh = procMeshes[procMeshI];
-        cellOffsets[procMeshI+1] =
-            cellOffsets[procMeshI]
-          + mesh.lduAddr().size();
-    }
-
-    // Faces initially get added in order, sorted later
-    labelList internalFaceOffsets(procMeshes.size()+1);
-    internalFaceOffsets[0] = 0;
-    forAll(procMeshes, procMeshI)
-    {
-        const lduMesh& mesh = procMeshes[procMeshI];
-        internalFaceOffsets[procMeshI+1] =
-            internalFaceOffsets[procMeshI]
-          + mesh.lduAddr().lowerAddr().size();
-    }
-
-    // Count how faces get added. Interfaces inbetween get merged.
-
-    // Merged interfaces: map from two processors back to
-    // - procMeshes
-    // - interface in procMesh
-    EdgeMap<labelPairList> mergedMap
-    (
-        2*procMeshes[0].interfaces().size()
-    );
-
-
-    // Unmerged interfaces: map from two processors back to
-    // - procMeshes
-    // - interface in procMesh
-    EdgeMap<labelPairList> unmergedMap
-    (
-        2*procMeshes[0].interfaces().size()
-    );
-    // Per interface the size
-    //DynamicList<label> interfaceSize(2*procMeshes[0].interfaces().size());
-
-    label nBoundaryFaces = 0;
-    label nInterfaces = 0;
-    labelList nCoupledFaces(procMeshes.size(), 0);
-
-    forAll(procMeshes, procMeshI)
-    {
-        const lduInterfacePtrsList interfaces =
-            procMeshes[procMeshI].interfaces();
-
-        // Get sorted order of processors
-        forAll(interfaces, intI)
-        {
-            if (interfaces.set(intI))
-            {
-                if (isA<processorGAMGInterface>(interfaces[intI]))
-                {
-                    const processorGAMGInterface& pldui =
-                        refCast<const processorGAMGInterface>
-                        (
-                            interfaces[intI]
-                        );
-                    if (pldui.myProcNo() != procIDs[procMeshI])
-                    {
-                        FatalErrorIn("combineMeshes()")
-                            << "proc:" << procIDs[procMeshI]
-                            << " myProcNo:" << pldui.myProcNo()
-                            << abort(FatalError);
-                    }
-
-
-                    const edge procEdge
-                    (
-                        min(pldui.myProcNo(), pldui.neighbProcNo()),
-                        max(pldui.myProcNo(), pldui.neighbProcNo())
-                    );
-
-
-                    label index = findIndex(procIDs, pldui.neighbProcNo());
-
-                    if (index == -1)
-                    {
-                        // Still external interface
-                        Pout<< "external interface: myProcNo:"
-                            << pldui.myProcNo()
-                            << " nbr:" << pldui.neighbProcNo()
-                            << " size:" << pldui.faceCells().size()
-                            << endl;
-
-                        nBoundaryFaces += pldui.faceCells().size();
-
-                        EdgeMap<labelPairList>::iterator iter =
-                            unmergedMap.find(procEdge);
-
-                        if (iter != unmergedMap.end())
-                        {
-                            iter().append(labelPair(procMeshI, intI));
-                        }
-                        else
-                        {
-                            unmergedMap.insert
-                            (
-                                procEdge,
-                                labelPairList(1, labelPair(procMeshI, intI))
-                            );
-                        }
-
-                        nInterfaces++;
-                    }
-                    else //if (pldui.myProcNo() < pldui.neighbProcNo())
-                    {
-                        // Merged interface
-                        Pout<< "merged interface: myProcNo:" << pldui.myProcNo()
-                            << " nbr:" << pldui.neighbProcNo()
-                            << " size:" << pldui.faceCells().size()
-                            << endl;
-                        if (pldui.myProcNo() < pldui.neighbProcNo())
-                        {
-                            nCoupledFaces[procMeshI] +=
-                                pldui.faceCells().size();
-                        }
-
-                        EdgeMap<labelPairList>::iterator iter =
-                            mergedMap.find(procEdge);
-
-                        if (iter != mergedMap.end())
-                        {
-                            iter().append(labelPair(procMeshI, intI));
-                        }
-                        else
-                        {
-                            mergedMap.insert
-                            (
-                                procEdge,
-                                labelPairList(1, labelPair(procMeshI, intI))
-                            );
-                        }
-                    }
-                }
-                else
-                {
-                    // Still external (non proc) interface
-                    FatalErrorIn("combineMeshes()")
-                        << "At mesh from processor " << procIDs[procMeshI]
-                        << " have interface " << intI
-                        << " of unhandled type " << interfaces[intI].type()
-                        << exit(FatalError);
-
-                    nBoundaryFaces += interfaces[intI].faceCells().size();
-                    nInterfaces++;
-                }
-            }
-        }
-    }
-
-
-
-    {
-        Pout<< "Remaining interfaces:" << endl;
-        forAllConstIter(EdgeMap<labelPairList>, unmergedMap, iter)
-        {
-            Pout<< "    procEdge:" << iter.key() << endl;
-            const labelPairList& elems = iter();
-            forAll(elems, i)
-            {
-                label procMeshI = elems[i][0];
-                label interfaceI = elems[i][1];
-                const lduInterfacePtrsList interfaces =
-                    procMeshes[procMeshI].interfaces();
-                const processorGAMGInterface& pldui =
-                    refCast<const processorGAMGInterface>
-                    (
-                        interfaces[interfaceI]
-                    );
-
-                Pout<< "        proc:" << procIDs[procMeshI]
-                    << " interfaceI:" << interfaceI
-                    << " between:" << pldui.myProcNo()
-                    << " and:" << pldui.neighbProcNo()
-                    << endl;
-            }
-            Pout<< endl;
-        }
-    }
-    {
-        Pout<< "Merged interfaces:" << endl;
-        forAllConstIter(EdgeMap<labelPairList>, mergedMap, iter)
-        {
-            Pout<< "    procEdge:" << iter.key() << endl;
-            const labelPairList& elems = iter();
-            forAll(elems, i)
-            {
-                label procMeshI = elems[i][0];
-                label interfaceI = elems[i][1];
-                const lduInterfacePtrsList interfaces =
-                    procMeshes[procMeshI].interfaces();
-                const processorGAMGInterface& pldui =
-                    refCast<const processorGAMGInterface>
-                    (
-                        interfaces[interfaceI]
-                    );
-
-                Pout<< "        proc:" << procIDs[procMeshI]
-                    << " interfaceI:" << interfaceI
-                    << " between:" << pldui.myProcNo()
-                    << " and:" << pldui.neighbProcNo()
-                    << endl;
-            }
-            Pout<< endl;
-        }
-    }
-
-
-    // Adapt faceOffsets for internal interfaces
-    faceOffsets.setSize(procMeshes.size()+1);
-    faceOffsets[0] = 0;
-    forAll(procMeshes, procMeshI)
-    {
-        faceOffsets[procMeshI+1] =
-            faceOffsets[procMeshI]
-          + procMeshes[procMeshI].lduAddr().lowerAddr().size() //internal
-          + nCoupledFaces[procMeshI];   // resolved coupled faces
-    }
-
-
-    // Combine upper and lower
-    labelList allLower(faceOffsets.last(), -1);
-    labelList allUpper(allLower.size(), -1);
-
-
-    // Old internal faces and resolved coupled interfaces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    forAll(procMeshes, procMeshI)
-    {
-        const labelUList& l = procMeshes[procMeshI].lduAddr().lowerAddr();
-        const labelUList& u = procMeshes[procMeshI].lduAddr().upperAddr();
-
-        // Add internal faces
-        label allFaceI = faceOffsets[procMeshI];
-
-        forAll(l, faceI)
-        {
-            allLower[allFaceI] = cellOffsets[procMeshI]+l[faceI];
-            allUpper[allFaceI] = cellOffsets[procMeshI]+u[faceI];
-            allFaceI++;
-        }
-
-        // Add converted interfaces
-        const lduInterfacePtrsList interfaces =
-            procMeshes[procMeshI].interfaces();
-
-        forAll(interfaces, intI)
-        {
-            if (interfaces.set(intI))
-            {
-                if (isA<processorGAMGInterface>(interfaces[intI]))
-                {
-                    const processorGAMGInterface& pldui =
-                        refCast<const processorGAMGInterface>
-                        (
-                            interfaces[intI]
-                        );
-
-                    // Look up corresponding interfaces
-                    label myP = pldui.myProcNo();
-                    label nbrP = pldui.neighbProcNo();
-
-                    if (myP < nbrP)
-                    {
-                        EdgeMap<labelPairList>::const_iterator fnd =
-                            mergedMap.find(edge(myP, nbrP));
-
-                        if (fnd != mergedMap.end())
-                        {
-                            const labelPairList& elems = fnd();
-
-                            // Find nbrP in elems
-                            label nbrProcMeshI = -1;
-                            label nbrIntI = -1;
-                            if (procIDs[elems[0][0]] == nbrP)
-                            {
-                                nbrProcMeshI = elems[0][0];
-                                nbrIntI = elems[0][1];
-                            }
-                            else
-                            {
-                                nbrProcMeshI = elems[1][0];
-                                nbrIntI = elems[1][1];
-                            }
-
-                            if
-                            (
-                                elems.size() != 2
-                             || procIDs[nbrProcMeshI] != nbrP
-                            )
-                            {
-                                FatalErrorIn("combineMeshes()")
-                                    << "elems:" << elems << abort(FatalError);
-                            }
-
-
-                            const lduInterfacePtrsList nbrInterfaces =
-                                procMeshes[nbrProcMeshI].interfaces();
-
-                            const processorGAMGInterface& nbrPldui =
-                                refCast<const processorGAMGInterface>
-                                (
-                                    nbrInterfaces[nbrIntI]
-                                );
-                            const labelUList& faceCells = pldui.faceCells();
-                            const labelUList& nbrFaceCells =
-                                nbrPldui.faceCells();
-
-                            forAll(faceCells, pfI)
-                            {
-                                allLower[allFaceI] =
-                                    cellOffsets[procMeshI]+faceCells[pfI];
-                                allUpper[allFaceI] =
-                                    cellOffsets[nbrProcMeshI]+nbrFaceCells[pfI];
-                                allFaceI++;
-                            }
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-
-    // Kept interfaces
-    // ~~~~~~~~~~~~~~~
-
-    lduInterfacePtrsList coarseInterfaces(nInterfaces);
-    label allInterfaceI = 0;
-
-    forAllConstIter(EdgeMap<labelPairList>, unmergedMap, iter)
-    {
-        Pout<< "procEdge:" << iter.key() << endl;
-        const labelPairList& elems = iter();
-
-        // Count
-        label n = 0;
-
-        forAll(elems, i)
-        {
-            label procMeshI = elems[i][0];
-            label interfaceI = elems[i][1];
-            const lduInterfacePtrsList interfaces =
-                procMeshes[procMeshI].interfaces();
-            n += interfaces[interfaceI].faceCells().size();
-        }
-
-        labelField allFaceCells(n);
-        labelField allFaceRestrictAddressing(n);
-        n = 0;
-
-        forAll(elems, i)
-        {
-            label procMeshI = elems[i][0];
-            label interfaceI = elems[i][1];
-            const lduInterfacePtrsList interfaces =
-                procMeshes[procMeshI].interfaces();
-
-            const labelUList& l = interfaces[interfaceI].faceCells();
-            forAll(l, faceI)
-            {
-                allFaceCells[n] = cellOffsets[procMeshI]+l[faceI];
-                allFaceRestrictAddressing[n] = n;
-                n++;
-            }
-        }
-
-
-        // Find out local and remote processor in new communicator
-        label minProcI = UPstream::procNo
-        (
-            newComm,
-            currentComm,
-            iter.key()[0]
-        );
-        label maxProcI = UPstream::procNo
-        (
-            newComm,
-            currentComm,
-            iter.key()[1]
-        );
-
-        coarseInterfaces.set
-        (
-            allInterfaceI,
-            new processorGAMGInterface
-            (
-                allInterfaceI,
-                coarseInterfaces,
-                allFaceCells,
-                allFaceRestrictAddressing,
-                newComm,
-                minProcI,
-                maxProcI,
-                tensorField(),          // forwardT
-                Pstream::msgType()      // tag
-            )
-        );
-        allInterfaceI++;
-    }
-
-
-    // Extract faceCells from coarseInterfaces
-    labelListList coarseInterfaceAddr(coarseInterfaces.size());
-    forAll(coarseInterfaces, coarseIntI)
-    {
-        if (coarseInterfaces.set(coarseIntI))
-        {
-            coarseInterfaceAddr[coarseIntI] =
-                coarseInterfaces[coarseIntI].faceCells();
-        }
-    }
-
-
-    Pout<< "allLower:" << allLower << endl;
-    Pout<< "allUpper:" << allUpper << endl;
-    checkUpperTriangular(cellOffsets.last(), allLower, allUpper);
-
-
-    autoPtr<lduPrimitiveMesh> meshPtr
-    (
-        new lduPrimitiveMesh
-        (
-            cellOffsets.last(),         //nCells
-            allLower,                   //lower
-            allUpper,                   //upper
-            coarseInterfaceAddr,        //faceCells
-            coarseInterfaces,           //interfaces
-            nonBlockingSchedule<processorGAMGInterface>(coarseInterfaces),
-            newComm,                    //communicator
-            true                        //reuse
-        )
-    );
-    return meshPtr;
-}
-
-
-
 void sendReceive
 (
     const label comm,
@@ -929,6 +433,150 @@ void setCommunicator(fvMesh& mesh, const label newComm)
 }
 
 
+namespace Foam
+{
+    typedef UPtrList<const GAMGInterfaceField> GAMGInterfaceFieldPtrsList;
+}
+
+
+// Gather matrices from processors procIDs[1..] on procIDs[0]
+void gatherMatrices
+(
+    const labelList& procIDs,
+    const PtrList<lduMesh>& procMeshes,
+
+    const lduMatrix& mat,
+    const FieldField<Field, scalar>& interfaceBouCoeffs,
+    const FieldField<Field, scalar>& interfaceIntCoeffs,
+    const lduInterfaceFieldPtrsList& interfaces,
+
+    PtrList<lduMatrix>& otherMats,
+    PtrList<FieldField<Field, scalar> >& otherBouCoeffs,
+    PtrList<FieldField<Field, scalar> >& otherIntCoeffs,
+    PtrList<GAMGInterfaceFieldPtrsList>& otherInterfaces
+)
+{
+    const label meshComm = mat.mesh().comm();
+
+    //lduInterfacePtrsList interfaces(mesh.interfaces());
+
+    if (Pstream::myProcNo(meshComm) == procIDs[0])
+    {
+        // Master.
+        otherMats.setSize(procIDs.size()-1);
+        otherBouCoeffs.setSize(procIDs.size()-1);
+        otherIntCoeffs.setSize(procIDs.size()-1);
+        otherInterfaces.setSize(procIDs.size()-1);
+
+        for (label i = 1; i < procIDs.size(); i++)
+        {
+            const lduMesh& procMesh = procMeshes[i];
+            const lduInterfacePtrsList procInterfaces = procMesh.interfaces();
+
+            IPstream fromSlave
+            (
+                Pstream::scheduled,
+                procIDs[i],
+                0,          // bufSize
+                Pstream::msgType(),
+                meshComm
+            );
+
+            otherMats.set(i-1, new lduMatrix(procMesh, fromSlave));
+
+            // Receive number of/valid interfaces
+            boolList validTransforms(fromSlave);
+            List<int> validRanks(fromSlave);
+
+            // Size coefficients
+            otherBouCoeffs.set
+            (
+                i-1,
+                new FieldField<Field, scalar>(validTransforms.size())
+            );
+            otherIntCoeffs.set
+            (
+                i-1,
+                new FieldField<Field, scalar>(validTransforms.size())
+            );
+            otherInterfaces.set
+            (
+                i-1,
+                new GAMGInterfaceFieldPtrsList(validTransforms.size())
+            );
+
+            forAll(validTransforms, intI)
+            {
+                if (validTransforms[intI])
+                {
+                    const processorGAMGInterface& interface =
+                        refCast<const processorGAMGInterface>
+                        (
+                            procInterfaces[intI]
+                        );
+
+
+                    otherBouCoeffs[i-1].set(intI, new scalarField(fromSlave));
+                    otherIntCoeffs[i-1].set(intI, new scalarField(fromSlave));
+                    otherInterfaces[i-1].set
+                    (
+                        intI,
+                        GAMGInterfaceField::New
+                        (
+                            interface,  //procInterfaces[intI],
+                            validTransforms[intI],
+                            validRanks[intI]
+                        ).ptr()
+                    );
+                }
+            }
+        }
+    }
+    else
+    {
+        // Send to master
+        OPstream toMaster
+        (
+            Pstream::scheduled,
+            procIDs[0],
+            0,
+            Pstream::msgType(),
+            meshComm
+        );
+
+
+        // Count valid interfaces
+        boolList validTransforms(interfaceBouCoeffs.size(), false);
+        List<int> validRanks(interfaceBouCoeffs.size(), -1);
+        forAll(interfaces, intI)
+        {
+            if (interfaces.set(intI))
+            {
+                const processorLduInterfaceField& interface =
+                    refCast<const processorLduInterfaceField>
+                    (
+                        interfaces[intI]
+                    );
+
+                validTransforms[intI] = interface.doTransform();
+                validRanks[intI] = interface.rank();
+            }
+        }
+
+        toMaster << mat << validTransforms << validRanks;
+        forAll(validTransforms, intI)
+        {
+            if (validTransforms[intI])
+            {
+                toMaster
+                    << interfaceBouCoeffs[intI]
+                    << interfaceIntCoeffs[intI];
+            }
+        }
+    }
+}
+
+
 int main(int argc, char *argv[])
 {
     #include "setRootCase.H"
@@ -965,405 +613,150 @@ int main(int argc, char *argv[])
     Pout<< "procIDs newComm:" << UPstream::procID(newComm) << endl;
 
 
-// On ALL processors: get the interfaces:
-lduInterfacePtrsList interfaces(mesh.interfaces());
-const lduAddressing& addressing = mesh.lduAddr();
-
-if (Pstream::myProcNo(newComm) != -1)
-{
-    print("InitialMesh", mesh);
-
-    labelList procIDs(2);
-    procIDs[0] = 0;
-    procIDs[1] = 1;
-
-    DynamicList<label> compactMap(interfaces.size());
-    forAll(interfaces, intI)
-    {
-        if (interfaces.set(intI))
-        {
-            compactMap.append(intI);
-        }
-    }
-
-    // Collect meshes from procs 0,1 (in newComm) on 1.
-    if (Pstream::myProcNo(newComm) == procIDs[0])
-    {
-        PtrList<lduMesh> procMeshes(procIDs.size());
-
-        // Master mesh (copy for now. TBD)
-        {
-            // Convert to processorGAMGInterfaces in new communicator
-            lduInterfacePtrsList compactInterfaces(compactMap.size());
-            labelListList compactPatchAddr(compactMap.size());
-
-            forAll(compactMap, compactI)
-            {
-                label intI = compactMap[compactI];
-
-                const processorLduInterface& pldui =
-                    refCast<const processorLduInterface>(interfaces[intI]);
-
-                compactPatchAddr[compactI] = interfaces[intI].faceCells();
-                labelList faceRestrictAddressing
-                (
-                    identity(interfaces[intI].faceCells().size())
-                );
-
-
-                Pout<< "for interface:" << intI
-                    << " comm:" << pldui.comm()
-                    << " myProcNo:" << pldui.myProcNo()
-                    << " in newComm:" << UPstream::procNo
-                        (
-                            newComm,
-                            pldui.comm(),
-                            pldui.myProcNo()
-                        )
-                    << endl;
-                Pout<< "for interface:" << intI
-                    << " comm:" << pldui.comm()
-                    << " neighbProcNo:" << pldui.neighbProcNo()
-                    << " in newComm:" << UPstream::procNo
-                        (
-                            newComm,
-                            pldui.comm(),
-                            pldui.neighbProcNo()
-                        )
-                    << endl;
-
-
-                compactInterfaces.set
-                (
-                    compactI,
-                    new processorGAMGInterface
-                    (
-                        compactI,
-                        compactInterfaces,
-                        compactPatchAddr[compactI],
-                        faceRestrictAddressing,
-
-                        newComm,                        //pldui.comm(),
-                        UPstream::procNo
-                        (
-                            newComm,
-                            pldui.comm(),
-                            pldui.myProcNo()
-                        ),
-                        UPstream::procNo
-                        (
-                            newComm,
-                            pldui.comm(),
-                            pldui.neighbProcNo()
-                        ),
-                        pldui.forwardT(),
-                        pldui.tag()
-                    )
-                );
-            }
-
-
-
-            procMeshes.set
-            (
-                0,
-                new lduPrimitiveMesh
-                (
-                    addressing.size(),
-                    addressing.lowerAddr(),
-                    addressing.upperAddr(),
-                    compactPatchAddr,
-                    compactInterfaces,
-                    nonBlockingSchedule<processorGAMGInterface>
-                    (
-                        compactInterfaces
-                    ),  //patch schedule,
-                    newComm
-                )
-            );
-        }
-
-        // Slave meshes
-        for (label i = 1; i < procIDs.size(); i++)
-        {
-            //Pout<< "on master :"
-            //    << " receiving from slave " << procIDs[i]
-            //    << endl;
-
-            IPstream fromSlave
-            (
-                Pstream::scheduled,
-                procIDs[i],
-                0,          // bufSize
-                Pstream::msgType(),
-                newComm
-            );
-
-            label nCells = readLabel(fromSlave);
-            labelList lowerAddr(fromSlave);
-            labelList upperAddr(fromSlave);
-            labelListList patchAddr(fromSlave);
-            labelList comm(fromSlave);
-            labelList myProcNo(fromSlave);
-            labelList neighbProcNo(fromSlave);
-            List<tensorField> forwardT(fromSlave);
-            labelList tag(fromSlave);
-
-            // Convert to processorGAMGInterfaces
-            lduInterfacePtrsList compactInterfaces(patchAddr.size());
-            forAll(patchAddr, compactI)
-            {
-                labelList faceRestrictAddressing
-                (
-                    identity(patchAddr[compactI].size())
-                );
-
-                compactInterfaces.set
-                (
-                    compactI,
-                    new processorGAMGInterface
-                    (
-                        compactI,
-                        compactInterfaces,
-                        patchAddr[compactI],
-                        faceRestrictAddressing,
-
-                        comm[compactI],
-                        UPstream::procNo
-                        (
-                            newComm,
-                            comm[compactI],
-                            myProcNo[compactI]
-                        ),
-                        UPstream::procNo
-                        (
-                            newComm,
-                            comm[compactI],
-                            neighbProcNo[compactI]
-                        ),
-
-                        forwardT[compactI],
-                        tag[compactI]
-                    )
-                );
-            }
-
-
-            procMeshes.set
-            (
-                i,
-                new lduPrimitiveMesh
-                (
-                    nCells,
-                    lowerAddr,
-                    upperAddr,
-                    patchAddr,
-                    compactInterfaces,
-                    nonBlockingSchedule<processorGAMGInterface>
-                    (
-                        compactInterfaces
-                    ),
-                    newComm
-                )
-            );
-        }
-
-
-        // Print a bit
-        forAll(procMeshes, i)
-        {
-            const lduMesh& pMesh = procMeshes[i];
-            print("procMesh" + Foam::name(i), pMesh);
-
-            const lduAddressing& addr = pMesh.lduAddr();
-            checkUpperTriangular
-            (
-                addr.size(),
-                addr.lowerAddr(),
-                addr.upperAddr()
-            );
-        }
-
-
-        // Combine
-        labelList cellOffsets;
-        labelList faceOffsets;
-        labelList interfaceOffsets;
-        autoPtr<lduPrimitiveMesh> allMeshPtr = combineMeshes
-        (
-            newComm,
-            procIDs,
-            procMeshes,
-
-            cellOffsets,     // per mesh the starting cell
-            faceOffsets,     // per mesh the starting face
-            interfaceOffsets // per mesh,per interface the starting face
-        );
-
-        print("ALLMESH", allMeshPtr());
-    }
-    else if (findIndex(procIDs, Pstream::myProcNo(newComm)) != -1)
-    {
-        // Send to master
-
-        //- Extract info from processorGAMGInterface
-        labelListList patchAddr(compactMap.size());
-        //labelListList faceRestrictAddr(compactMap.size());
-        labelList comm(compactMap.size());
-        labelList myProcNo(compactMap.size());
-        labelList neighbProcNo(compactMap.size());
-        List<tensorField> forwardT(compactMap.size());
-        labelList tag(compactMap.size());
-
-        forAll(compactMap, compactI)
-        {
-            label intI = compactMap[compactI];
-
-            const processorLduInterface& pldui =
-                refCast<const processorLduInterface>(interfaces[intI]);
-
-            patchAddr[compactI] = interfaces[intI].faceCells();
-            //faceRestrictAddr[compactI] = pldui.faceRestrictAddresssing();
-            comm[compactI] = pldui.comm();
-            myProcNo[compactI] = pldui.myProcNo();
-            neighbProcNo[compactI] =  pldui.neighbProcNo();
-            forwardT[compactI] =  pldui.forwardT();
-            tag[compactI] =  pldui.tag();
-        }
-
-        //Pout<< "on slave sending to " << Pstream::masterNo() << endl;
-        //Pout<< "on slave nCells:" << addressing.size() << endl;
-        //Pout<< "on slave lowerAdd:" << addressing.lowerAddr().size() << endl;
-        //Pout<< "on slave upperAddr:" << addressing.upperAddr().size() << endl;
-        //Pout<< "on slave patchAddr:" << patchAddr.size() << endl;
-        //Pout<< "on slave comm:" << comm << endl;
-        //Pout<< "on slave myProcNo:" << myProcNo << endl;
-        //Pout<< "on slave neighbProcNo:" << neighbProcNo << endl;
-        //Pout<< "on slave forwardT:" << forwardT << endl;
-        //Pout<< "on slave tag:" << tag << endl;
-
-        OPstream toMaster
-        (
-            Pstream::scheduled,
-            procIDs[0],
-            0,
-            Pstream::msgType(),
-            newComm
-        );
-        toMaster
-            << addressing.size()
-            << addressing.lowerAddr()
-            << addressing.upperAddr()
-            << patchAddr
-            //<< faceRestrictAddr
-            << comm
-            << myProcNo
-            << neighbProcNo
-            << forwardT
-            << tag;
-    }
-}
-return 0;
-
-
-
+//// On ALL processors: get the interfaces:
+//lduInterfacePtrsList interfaces(mesh.interfaces());
+//PtrList<lduMesh> procMeshes;
+//
 //if (Pstream::myProcNo(newComm) != -1)
 //{
-//    scalarField diagonal(11+Pstream::myProcNo(), 11+Pstream::myProcNo());
-//    scalarField upper(22+Pstream::myProcNo(), 22+Pstream::myProcNo());
-//    scalarField lower(upper.size(), 22+Pstream::myProcNo());
+//    print("InitialMesh", mesh);
 //
-//    // Collect the local sizes
-//    const globalIndex cellOffsets
-//    (
-//        diagonal.size(),
-//        Pstream::msgType(),
-//        newComm,
-//        true
-//    );
+//    labelList procIDs(3);
+//    procIDs[0] = 0;
+//    procIDs[1] = 1;
+//    procIDs[2] = 2;
 //
-//    const globalIndex faceOffsets
-//    (
-//        upper.size(),
-//        Pstream::msgType(),
-//        newComm,
-//        true
-//    );
+////XXXXXX
+//    // Collect meshes from procs 0,1 (in newComm) on 1.
+//    lduPrimitiveMesh::gather(mesh, procIDs, procMeshes);
 //
-//
-//    scalarField allDiagonal;
-//    scalarField allUpper;
-//    scalarField allLower;
-//
-//    collect
-//    (
-//        newComm,
-//        cellOffsets,
-//        faceOffsets,
-//
-//        diagonal,
-//        upper,
-//        lower,
-//
-//        allDiagonal,
-//        allUpper,
-//        allLower
-//    );
-//
-//    Pout<< "diagonal:" << diagonal << endl;
-//    Pout<< "allDiagonal:" << allDiagonal << endl;
-//
-//    Pout<< "upper:" << upper << endl;
-//    Pout<< "allUpper:" << allUpper << endl;
-//
-//    Pout<< "lower:" << lower << endl;
-//    Pout<< "allLower:" << allLower << endl;
-//
-//}
-
-
-
-//if (Pstream::myProcNo(newComm) != -1)
-//{
-//    // Interface coefficients. The destination number of interfaces can differ
-//    // We have pre-sorted the overall match so we can just slot them in.
-//    FieldField<Field, scalar> interfaceBouCoeffs;
-//    //FieldField<Field, scalar> interfaceIntCoeffs;
+//    if (Pstream::myProcNo(newComm) == procIDs[0])
 //    {
-//        label nInt = 3+Pstream::myProcNo();
-//        interfaceBouCoeffs.setSize(nInt);
-//        forAll(interfaceBouCoeffs, intI)
+//        // Print a bit
+//        forAll(procMeshes, i)
 //        {
-//            interfaceBouCoeffs.set(intI, new scalarField(intI+1));
-//            scalarField& fld = interfaceBouCoeffs[intI];
-//            fld = Pstream::myProcNo();
+//            const lduMesh& pMesh = procMeshes[i];
+//            print("procMesh" + Foam::name(i), pMesh);
+//
+//            const lduAddressing& addr = pMesh.lduAddr();
+//            checkUpperTriangular
+//            (
+//                addr.size(),
+//                addr.lowerAddr(),
+//                addr.upperAddr()
+//            );
+//        }
+//
+//
+//        // Combine
+//        labelList cellOffsets;
+//        labelListList faceMap;
+//        labelListList boundaryMap;
+//        labelListListList boundaryFaceMap;
+//        //autoPtr<lduPrimitiveMesh> allMeshPtr = combineMeshes
+//        //(
+//        //    newComm,
+//        //    procIDs,
+//        //    procMeshes,
+//        //
+//        //    cellOffsets,        // per mesh the starting cell
+//        //    faceMap,            // per mesh the starting face
+//        //    boundaryMap,        // per mesh,per interface the starting face
+//        //    boundaryFaceMap
+//        //);
+//        //const lduPrimitiveMesh& allMesh = allMeshPtr();
+//        const lduPrimitiveMesh allMesh
+//        (
+//            newComm,
+//            procIDs,
+//            procMeshes,
+//
+//            cellOffsets,
+//            faceMap,
+//            boundaryMap,
+//            boundaryFaceMap
+//        );
+//
+//
+//        print("ALLMESH", allMesh);
+//
+//        forAll(procMeshes, procMeshI)
+//        {
+//            const lduMesh& pMesh = procMeshes[procMeshI];
+//            //const lduAddressing& pAddressing = pMesh.lduAddr();
+//
+//            Pout<< "procMesh" << procMeshI << endl
+//                << "    cells start at:" << cellOffsets[procMeshI] << endl
+//                << "    faces to to:" << faceMap[procMeshI] << endl;
+//
+//            lduInterfacePtrsList interfaces = pMesh.interfaces();
+//            forAll(interfaces, intI)
+//            {
+//                Pout<< "    patch:" << intI
+//                    << " becomes patch:" << boundaryMap[procMeshI][intI]
+//                    << endl;
+//
+//                Pout<< "    patch:" << intI
+//                    << " faces become faces:"
+//                    << boundaryFaceMap[procMeshI][intI]
+//                    << endl;
+//            }
 //        }
 //    }
 //
-//    const globalIndex interfaceBouOffsets
-//    (
-//        interfaceBouCoeffs.size(),
-//        Pstream::msgType(),
-//        newComm,
-//        true
-//    );
 //
-//    FieldField<Field, scalar> allInterfaceBouCoeffs;
-//    sendReceive
+//    // Construct test data
+//    // ~~~~~~~~~~~~~~~~~~~
+//
+//    GAMGInterfaceFieldPtrsList interfaces(interfaces.size());
+//    FieldField<Field, scalar> interfaceBouCoeffs(interfaces.size());
+//    FieldField<Field, scalar> interfaceIntCoeffs(interfaces.size());
+//
+//    forAll(interfaces, intI)
+//    {
+//        if (interfaces.set(intI))
+//        {
+//            label size = interfaces[intI].size();
+//
+//            interfaces.set
+//            (
+//                intI,
+//                GAMGInterfaceField::New
+//                (
+//                    mesh.interfaces()[intI],
+//                    interfaces[intI]
+//                )
+//            );
+//            interfaceBouCoeffs.set(intI, new scalarField(size, 111));
+//            interfaceIntCoeffs.set(intI, new scalarField(size, 222));
+//        }
+//    }
+//
+//
+//    PtrList<lduMatrix> otherMats;
+//    PtrList<FieldField<Field, scalar> > otherBouCoeffs;
+//    PtrList<FieldField<Field, scalar> > otherIntCoeffs;
+//    PtrList<GAMGInterfaceFieldPtrsList> otherInterfaces;
+//    gatherMatrices
 //    (
-//        newComm,
-//        Pstream::msgType(),
-//        interfaceBouOffsets,
+//        procIDs,
+//        procMeshes,
+//
+//        mat,
 //        interfaceBouCoeffs,
+//        interfaceIntCoeffs,
+//        interfaces,
 //
-//        allInterfaceBouCoeffs
+//        otherMats,
+//        otherBouCoeffs,
+//        otherIntCoeffs,
+//        otherInterfaces
 //    );
-//
-//    Pout<< "interfaceBouCoeffs   :" << interfaceBouCoeffs << endl;
-//    Pout<< "allInterfaceBouCoeffs:" << allInterfaceBouCoeffs << endl;
-//
+////XXXXXX
 //}
 
+
+
     {
         Pout<< "World:" << UPstream::worldComm
             << " procID:" << 2
@@ -1372,30 +765,6 @@ return 0;
             << " rank2:" << UPstream::procNo(newComm, UPstream::worldComm, 2)
             << endl;
     }
-    //
-    //
-    //{
-    //    setCommunicator(mesh, newComm);
-    //    Pout<< "Mesh comm:" << mesh.comm() << endl;
-    //    forAll(pbm, patchI)
-    //    {
-    //        if (isA<processorPolyPatch>(pbm[patchI]))
-    //        {
-    //            const processorPolyPatch& ppp =
-    //            refCast
-    //            <
-    //                const processorPolyPatch
-    //            >(pbm[patchI]);
-    //
-    //
-    //            Pout<< "    " << ppp.name()
-    //                << " myRank:" << ppp.myProcNo()
-    //                << " nbrRank:" << ppp.neighbProcNo()
-    //                << endl;
-    //        }
-    //    }
-    //    setCommunicator(mesh, UPstream::worldComm);
-    //}
 
 
     while (simple.loop())
@@ -1406,7 +775,8 @@ return 0;
         {
             fvScalarMatrix Teqn
             (
-                fvm::ddt(T) - fvm::laplacian(DT, T)
+                //fvm::ddt(T) - fvm::laplacian(DT, T)
+                fvm::laplacian(DT, T)
             );