ENH: GAMGAgglomeration: moved lduMesh agglomeration into GAMGAgglomeration

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C

@@ -255,4 +255,101 @@ const Foam::lduInterfacePtrsList& Foam::GAMGAgglomeration::interfaceLevel
 }
 
 
+void Foam::GAMGAgglomeration::gatherMeshes
+(
+    const labelList& procIDs,
+    const label allMeshComm
+) const
+{
+    const lduMesh& coarsestMesh = meshLevels_.last();
+
+    label coarseComm = coarsestMesh.comm();
+
+    label oldWarn = UPstream::warnComm;
+    UPstream::warnComm = coarseComm;
+
+    Pout<< "GAMGAgglomeration : gathering coarsestmesh (level="
+        << meshLevels_.size()-1
+        << ") using communicator " << coarseComm << endl;
+
+    lduPrimitiveMesh::gather(coarsestMesh, procIDs, otherMeshes_);
+
+    Pout<< "** Own Mesh " << coarsestMesh.info() << endl;
+    forAll(otherMeshes_, i)
+    {
+        Pout<< "** otherMesh " << i << " "
+            << otherMeshes_[i].info()
+            << endl;
+    }
+    Pout<< endl;
+
+    if (Pstream::myProcNo(coarseComm) == procIDs[0])
+    {
+        // Agglomerate all addressing
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        allMeshPtr_.reset
+        (
+            new lduPrimitiveMesh
+            (
+                allMeshComm,
+                procIDs,
+                coarsestMesh,
+                otherMeshes_,
+
+                cellOffsets_,
+                faceMap_,
+                boundaryMap_,
+                boundaryFaceMap_
+            )
+        );
+    }
+
+    UPstream::warnComm = oldWarn;
+}
+
+
+const Foam::lduPrimitiveMesh& Foam::GAMGAgglomeration::allMesh() const
+{
+    if (!allMeshPtr_.valid())
+    {
+        FatalErrorIn("GAMGAgglomeration::allMesh() const")
+            << "Processor-agglomerated mesh not constructed. Did you call"
+            << " gatherMeshes?"
+            << exit(FatalError);
+    }
+    return allMeshPtr_();
+}
+
+
+const Foam::labelList& Foam::GAMGAgglomeration::cellOffsets() const
+{
+    return cellOffsets_;
+}
+
+
+const Foam::labelListList& Foam::GAMGAgglomeration::faceMap() const
+{
+    return faceMap_;
+}
+
+
+const Foam::labelListList& Foam::GAMGAgglomeration::boundaryMap() const
+{
+    return boundaryMap_;
+}
+
+
+const Foam::labelListListList& Foam::GAMGAgglomeration::boundaryFaceMap() const
+{
+    return boundaryFaceMap_;
+}
+
+
+const Foam::PtrList<Foam::lduMesh>& Foam::GAMGAgglomeration::otherMeshes() const
+{
+    return otherMeshes_;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H

@@ -99,6 +99,30 @@ protected:
         //  Warning: Needs to be deleted explicitly.
         PtrList<lduInterfacePtrsList> interfaceLevels_;
 
+
+        //- Processor agglomeration
+
+            //- Combined coarsest meshes for all processors
+            mutable autoPtr<lduPrimitiveMesh> allMeshPtr_;
+
+            //- Per remote processor the local equivalent. TBD: hack. Can be
+            //  deleted once lduMatrices get received and merged into single
+            //  one without constructing local ones
+            mutable PtrList<lduMesh> otherMeshes_;
+
+            //- Mapping from processor to coarsestAllMesh cells
+            mutable labelList cellOffsets_;
+
+            //- Mapping from processor to coarsestAllMesh face
+            mutable labelListList faceMap_;
+
+            //- Mapping from processor to coarsestAllMesh boundary
+            mutable labelListList boundaryMap_;
+
+            //- Mapping from processor to coarsestAllMesh boundary face
+            mutable labelListListList boundaryFaceMap_;
+
+
         //- Assemble coarse mesh addressing
         void agglomerateLduAddressing(const label fineLevelIndex);
 
@@ -250,6 +274,38 @@ public:
                 const Field<Type>& cf,
                 const label coarseLevelIndex
             ) const;
+
+
+        // Procesor agglomeration
+
+            //- Collect and combine processor meshes into allMesh. Pass in
+            //  communicator for combined mesh.
+            void gatherMeshes
+            (
+                const labelList& procIDs,
+                const label allMeshComm
+            ) const;
+
+            //- Combined mesh
+            const lduPrimitiveMesh& allMesh() const;
+
+            //- Mapping from processor to coarsestAllMesh cells
+            const labelList& cellOffsets() const;
+
+            //- Mapping from processor to coarsestAllMesh face
+            const labelListList& faceMap() const;
+
+            //- Mapping from processor to coarsestAllMesh boundary
+            const labelListList& boundaryMap() const;
+
+            //- Mapping from processor to coarsestAllMesh boundary face
+            const labelListListList& boundaryFaceMap() const;
+
+            //- Per remote processor the local equivalent. TBD: hack. Can be
+            //  deleted once lduMatrices get received and merged into single
+            //  one without constructing local ones
+            const PtrList<lduMesh>& otherMeshes() const;
+
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C

@@ -45,10 +45,12 @@ namespace Foam
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
+// Gather matrices.
+// Note: matrices get constructed with dummy mesh
 void Foam::GAMGSolver::gatherMatrices
 (
     const labelList& procIDs,
-    const PtrList<lduMesh>& procMeshes,
+    const lduMesh& dummyMesh,
     const label meshComm,
 
     const lduMatrix& mat,
@@ -59,16 +61,14 @@ void Foam::GAMGSolver::gatherMatrices
     PtrList<lduMatrix>& otherMats,
     PtrList<FieldField<Field, scalar> >& otherBouCoeffs,
     PtrList<FieldField<Field, scalar> >& otherIntCoeffs,
-    PtrList<lduInterfaceFieldPtrsList>& otherInterfaces
+    List<boolList>& otherTransforms,
+    List<List<int> >& otherRanks
 ) const
 {
     Pout<< "GAMGSolver::gatherMatrices :"
         << " collecting matrices on procs:" << procIDs
         << " using comm:" << meshComm << endl;
 
-
-    //lduInterfacePtrsList interfaces(mesh.interfaces());
-
     if (Pstream::myProcNo(meshComm) == procIDs[0])
     {
         // Master.
@@ -76,85 +76,59 @@ void Foam::GAMGSolver::gatherMatrices
         Pout<< "GAMGSolver::gatherMatrices :"
             << " master:" << Pstream::myProcNo(meshComm) << endl;
 
-        for (label i=0; i < procMeshes.size(); i++)
-        {
-            if (meshComm != procMeshes[i].comm())
-            {
-                FatalErrorIn("GAMGSolver::gatherMatrices()")
-                    << "Inconsistent communicator :"
-                    << "master processor:" << procIDs[0]
-                    << " comm:" << meshComm
-                    << " slave processor:" << procIDs[i]
-                    << " comm:" << procMeshes[i].comm()
-                    << abort(FatalError);
-            }
-        }
-
 
         otherMats.setSize(procIDs.size()-1);
         otherBouCoeffs.setSize(procIDs.size()-1);
         otherIntCoeffs.setSize(procIDs.size()-1);
-        otherInterfaces.setSize(procIDs.size()-1);
+        otherTransforms.setSize(procIDs.size()-1);
+        otherRanks.setSize(procIDs.size()-1);
 
-        for (label i = 1; i < procIDs.size(); i++)
+        for (label procI = 1; procI < procIDs.size(); procI++)
         {
-            const lduMesh& procMesh = procMeshes[i];
+            label otherI = procI-1;
-            const lduInterfacePtrsList procInterfaces = procMesh.interfaces();
 
             IPstream fromSlave
             (
                 Pstream::scheduled,
-                procIDs[i],
+                procIDs[procI],
                 0,          // bufSize
                 Pstream::msgType(),
                 meshComm
             );
 
-            otherMats.set(i-1, new lduMatrix(procMesh, fromSlave));
+            otherMats.set(otherI, new lduMatrix(dummyMesh, fromSlave));
 
             // Receive number of/valid interfaces
-            boolList validTransforms(fromSlave);
+            boolList& procTransforms = otherTransforms[otherI];
-            List<int> validRanks(fromSlave);
+            List<int>& procRanks = otherRanks[otherI];
+
+            fromSlave >> procTransforms;
+            fromSlave >> procRanks;
 
             // Size coefficients
             otherBouCoeffs.set
             (
-                i-1,
+                otherI,
-                new FieldField<Field, scalar>(validTransforms.size())
+                new FieldField<Field, scalar>(procRanks.size())
             );
             otherIntCoeffs.set
             (
-                i-1,
+                otherI,
-                new FieldField<Field, scalar>(validTransforms.size())
+                new FieldField<Field, scalar>(procRanks.size())
             );
-            otherInterfaces.set
+            forAll(procRanks, intI)
-            (
-                i-1,
-                new lduInterfaceFieldPtrsList(validTransforms.size())
-            );
-
-            forAll(validTransforms, intI)
             {
-                if (validTransforms[intI])
+                if (procRanks[intI] != -1)
                 {
-                    const processorGAMGInterface& interface =
+                    otherBouCoeffs[otherI].set
-                        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
+                        new scalarField(fromSlave)
-                        (
+                    );
-                            interface,  //procInterfaces[intI],
+                    otherIntCoeffs[otherI].set
-                            validTransforms[intI],
+                    (
-                            validRanks[intI]
+                        intI,
-                        ).ptr()
+                        new scalarField(fromSlave)
                     );
                 }
             }
@@ -169,8 +143,8 @@ void Foam::GAMGSolver::gatherMatrices
         // Send to master
 
         // Count valid interfaces
-        boolList validTransforms(interfaceBouCoeffs.size(), false);
+        boolList procTransforms(interfaceBouCoeffs.size(), false);
-        List<int> validRanks(interfaceBouCoeffs.size(), -1);
+        List<int> procRanks(interfaceBouCoeffs.size(), -1);
         forAll(interfaces, intI)
         {
             if (interfaces.set(intI))
@@ -181,8 +155,8 @@ void Foam::GAMGSolver::gatherMatrices
                         interfaces[intI]
                     );
 
-                validTransforms[intI] = interface.doTransform();
+                procTransforms[intI] = interface.doTransform();
-                validRanks[intI] = interface.rank();
+                procRanks[intI] = interface.rank();
             }
         }
 
@@ -198,10 +172,10 @@ void Foam::GAMGSolver::gatherMatrices
             meshComm
         );
 
-        toMaster << mat << validTransforms << validRanks;
+        toMaster << mat << procTransforms << procRanks;
-        forAll(validTransforms, intI)
+        forAll(procRanks, intI)
         {
-            if (validTransforms[intI])
+            if (procRanks[intI] != -1)
             {
                 toMaster
                     << interfaceBouCoeffs[intI]
@@ -310,60 +284,8 @@ Foam::GAMGSolver::GAMGSolver
 
             // All processors to master
             const List<int>& procIDs = UPstream::procID(coarseComm);
-
-            // Gather all meshes onto procIDs[0]
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            // Note: move into GAMGAgglomeration
-
             Pout<< "procIDs:" << procIDs << endl;
 
-            PtrList<lduMesh> procMeshes;
-            lduPrimitiveMesh::gather(coarsestMesh, procIDs, procMeshes);
-
-            forAll(procMeshes, procI)
-            {
-                Pout<< "** procMesh " << procI << " "
-                    << procMeshes[procI].info()
-                    << endl;
-            }
-            Pout<< endl;
-
-
-            // Gather all matrix coefficients onto procIDs[0]
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-            PtrList<lduMatrix> otherMats;
-            PtrList<FieldField<Field, scalar> > otherBouCoeffs;
-            PtrList<FieldField<Field, scalar> > otherIntCoeffs;
-            PtrList<lduInterfaceFieldPtrsList> otherInterfaces;
-            gatherMatrices
-            (
-                procIDs,
-                procMeshes,
-                coarseComm,
-
-                coarsestMatrix,
-                coarsestBouCoeffs,
-                coarsestIntCoeffs,
-                coarsestInterfaces,
-
-                otherMats,
-                otherBouCoeffs,
-                otherIntCoeffs,
-                otherInterfaces
-            );
-
-            Pout<< "Own matrix:" << coarsestMatrix.info() << endl;
-
-            forAll(otherMats, i)
-            {
-                Pout<< "** otherMats " << i << " "
-                    << otherMats[i].info()
-                    << endl;
-            }
-            Pout<< endl;
-
-
             // Allocate a communicator for single processor
             label masterComm = UPstream::allocateCommunicator
             (
@@ -375,31 +297,56 @@ Foam::GAMGSolver::GAMGSolver
                 << " out of " << procIDs << endl;
 
 
+            // Gather all meshes onto procIDs[0]
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            agglomeration_.gatherMeshes(procIDs, masterComm);
+
+
+            // Gather all matrix coefficients onto procIDs[0]
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            PtrList<lduMatrix> otherMats;
+            PtrList<FieldField<Field, scalar> > otherBouCoeffs;
+            PtrList<FieldField<Field, scalar> > otherIntCoeffs;
+            List<boolList> otherTransforms;
+            List<List<int> > otherRanks;
+            gatherMatrices
+            (
+                procIDs,
+                coarsestMesh,
+                coarseComm,
+
+                coarsestMatrix,
+                coarsestBouCoeffs,
+                coarsestIntCoeffs,
+                coarsestInterfaces,
+
+                otherMats,
+                otherBouCoeffs,
+                otherIntCoeffs,
+                otherTransforms,
+                otherRanks
+            );
+
+
             if (Pstream::myProcNo(coarseComm) == procIDs[0])
             {
-                // Agglomerate all addressing
-                // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-                allMeshPtr_.reset
-                (
-                    new lduPrimitiveMesh
-                    (
-                        masterComm,
-                        procIDs,
-                        procMeshes,
-
-                        cellOffsets_,
-                        faceMap_,
-                        boundaryMap_,
-                        boundaryFaceMap_
-                    )
-                );
-                const lduMesh& allMesh = allMeshPtr_();
-
-
                 // Agglomerate all matrix
                 // ~~~~~~~~~~~~~~~~~~~~~~
 
+                Pout<< "Own matrix:" << coarsestMatrix.info() << endl;
+
+                forAll(otherMats, i)
+                {
+                    Pout<< "** otherMats " << i << " "
+                        << otherMats[i].info()
+                        << endl;
+                }
+                Pout<< endl;
+
+
+                const lduMesh& allMesh = agglomeration_.allMesh();
+
                 allMatrixPtr_.reset(new lduMatrix(allMesh));
                 lduMatrix& allMatrix = allMatrixPtr_();
 
@@ -420,7 +367,7 @@ Foam::GAMGSolver::GAMGSolver
                         (
                             allDiag,
                             otherMats[i].diag().size(),
-                            cellOffsets_[i+1]
+                            agglomeration_.cellOffsets()[i+1]
                         ).assign
                         (
                             otherMats[i].diag()
@@ -430,23 +377,35 @@ Foam::GAMGSolver::GAMGSolver
                 if (coarsestMatrix.hasLower())
                 {
                     scalarField& allLower = allMatrix.lower();
-                    UIndirectList<scalar>(allLower, faceMap_[0]) =
+                    UIndirectList<scalar>
-                        coarsestMatrix.lower();
+                    (
+                        allLower,
+                        agglomeration_.faceMap()[0]
+                    ) = coarsestMatrix.lower();
                     forAll(otherMats, i)
                     {
-                        UIndirectList<scalar>(allLower, faceMap_[i+1]) =
+                        UIndirectList<scalar>
-                            otherMats[i].lower();
+                        (
+                            allLower,
+                            agglomeration_.faceMap()[i+1]
+                        ) = otherMats[i].lower();
                     }
                 }
                 if (coarsestMatrix.hasUpper())
                 {
                     scalarField& allUpper = allMatrix.upper();
-                    UIndirectList<scalar>(allUpper, faceMap_[0]) =
+                    UIndirectList<scalar>
-                        coarsestMatrix.upper();
+                    (
+                        allUpper,
+                        agglomeration_.faceMap()[0]
+                    ) = coarsestMatrix.upper();
                     forAll(otherMats, i)
                     {
-                        UIndirectList<scalar>(allUpper, faceMap_[i+1]) =
+                        UIndirectList<scalar>
-                            otherMats[i].upper();
+                        (
+                            allUpper,
+                            agglomeration_.faceMap()[i+1]
+                        ) = otherMats[i].upper();
                     }
                 }
 
@@ -470,7 +429,7 @@ Foam::GAMGSolver::GAMGSolver
                 }
 
                 labelList nBounFaces(allMeshInterfaces.size());
-                forAll(boundaryMap_, procI)
+                forAll(agglomeration_.boundaryMap(), procI)
                 {
                     const FieldField<Field, scalar>& procBouCoeffs
                     (
@@ -484,15 +443,8 @@ Foam::GAMGSolver::GAMGSolver
                       ? coarsestIntCoeffs
                       : otherIntCoeffs[procI-1]
                     );
-                    const lduInterfaceFieldPtrsList procInterfaces
-                    (
-                        (procI == 0)
-                      ? coarsestInterfaces
-                      : otherInterfaces[procI-1]
-                    );
 
-
+                    const labelList& bMap = agglomeration_.boundaryMap()[procI];
-                    const labelList& bMap = boundaryMap_[procI];
                     forAll(bMap, procIntI)
                     {
                         label allIntI = bMap[procIntI];
@@ -506,11 +458,27 @@ Foam::GAMGSolver::GAMGSolver
                             {
                                 // Construct lduInterfaceField
 
-                                const processorGAMGInterfaceField& procInt =
+                                bool doTransform = false;
-                                    refCast<const processorGAMGInterfaceField>
+                                int rank = -1;
-                                    (
+                                if (procI == 0)
-                                        procInterfaces[procIntI]
+                                {
-                                    );
+                                    const processorGAMGInterfaceField& procInt =
+                                        refCast
+                                        <
+                                            const processorGAMGInterfaceField
+                                        >
+                                        (
+                                            coarsestInterfaces[procIntI]
+                                        );
+                                    doTransform = procInt.doTransform();
+                                    rank = procInt.rank();
+                                }
+                                else
+                                {
+                                    doTransform =
+                                        otherTransforms[procI-1][procIntI];
+                                    rank = otherRanks[procI-1][procIntI];
+                                }
 
                                 allInterfaces_.set
                                 (
@@ -521,8 +489,8 @@ Foam::GAMGSolver::GAMGSolver
                                         (
                                             allMeshInterfaces[allIntI]
                                         ),
-                                        procInt.doTransform(),
+                                        doTransform,
-                                        procInt.rank()
+                                        rank
                                     ).ptr()
                                 );
                             }
@@ -535,8 +503,8 @@ Foam::GAMGSolver::GAMGSolver
                             scalarField& allInt =
                                 allInterfaceIntCoeffs_[allIntI];
 
-                            const labelList& map =
+                            const labelList& map = agglomeration_.
-                                boundaryFaceMap_[procI][procIntI];
+                                boundaryFaceMap()[procI][procIntI];
                             const scalarField& procBou =
                                 procBouCoeffs[procIntI];
                             const scalarField& procInt =
@@ -558,27 +526,15 @@ Foam::GAMGSolver::GAMGSolver
                         {
                             // Boundary has become internal face
 
-                            // Problem: we don't know whether face is flipped
+                            const labelList& map = agglomeration_.
-                            // or not. Either we interrogate the original
+                                boundaryFaceMap()[procI][procIntI];
-                            // mesh (requires procMeshes which we don't want)
-                            // or we need to store this information somehow
-                            // in the mapping (boundaryMap_, boundaryFaceMap_)
-
-
-                            //Hack. See above.
-                            //const labelList& fCells =
-                            //    procInterfaces[procI].faceCells();
-                            const labelList& fCells =
-                                procMeshes[procI].lduAddr().patchAddr(procIntI);
-
-                            const labelList& map =
-                                boundaryFaceMap_[procI][procIntI];
 
                             const labelList& allU =
                                 allMesh.lduAddr().upperAddr();
                             const labelList& allL =
                                 allMesh.lduAddr().lowerAddr();
-                            const label off = cellOffsets_[procI];
+                            const label off = agglomeration_.
+                                cellOffsets()[procI];
 
                             const scalarField& procBou =
                                 procBouCoeffs[procIntI];
@@ -594,12 +550,12 @@ Foam::GAMGSolver::GAMGSolver
                                     if (coarsestMatrix.hasUpper())
                                     {
                                         allMatrix.upper()[allFaceI] =
-                                            procBou[i];
+                                            -procBou[i];
                                     }
                                     if (coarsestMatrix.hasLower())
                                     {
                                         allMatrix.lower()[allFaceI] =
-                                            procInt[i];
+                                            -procInt[i];
                                     }
                                 }
                                 else
@@ -609,12 +565,12 @@ Foam::GAMGSolver::GAMGSolver
                                     if (coarsestMatrix.hasUpper())
                                     {
                                         allMatrix.upper()[allFaceI] =
-                                            procInt[i];
+                                            -procInt[i];
                                     }
                                     if (coarsestMatrix.hasLower())
                                     {
                                         allMatrix.lower()[allFaceI] =
-                                            procBou[i];
+                                            -procBou[i];
                                     }
                                 }
 
@@ -626,6 +582,15 @@ Foam::GAMGSolver::GAMGSolver
                                   ? map[i]
                                   : -map[i]-1
                                 );
+
+                                const labelList& fCells =
+                                    lduPrimitiveMesh::mesh
+                                    (
+                                        coarsestMesh,
+                                        agglomeration_.otherMeshes(),
+                                        procI
+                                    ).lduAddr().patchAddr(procIntI);
+
                                 label allCellI = off + fCells[i];
 
                                 if
@@ -649,6 +614,8 @@ Foam::GAMGSolver::GAMGSolver
                         }
                     }
                 }
+
+                Pout<< "** Assembled allMatrix:" << allMatrix.info() << endl;
             }
 
             UPstream::warnComm = oldWarn;

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H

@@ -131,20 +131,6 @@ class GAMGSolver
         //- LU decompsed coarsest matrix
         autoPtr<LUscalarMatrix> coarsestLUMatrixPtr_;
 
-
-
-        // Processor addressing agglomeration - to be moved into
-        // GAMGAgglomeration
-
-            //- Combined coarsest meshes for all processors
-            autoPtr<lduPrimitiveMesh> allMeshPtr_;
-
-            //- Mapping from processor to coarsestAllMesh
-            labelList cellOffsets_;
-            labelListList faceMap_;
-            labelListList boundaryMap_;
-            labelListListList boundaryFaceMap_;
-
         // Processor matrix agglomeration
 
             //- Combined coarsest level matrix for all processors
@@ -187,7 +173,8 @@ class GAMGSolver
         void gatherMatrices
         (
             const labelList& procIDs,
-            const PtrList<lduMesh>& procMeshes,
+            const lduMesh& myMesh,
+            const PtrList<lduMesh>& otherMeshes,
             const label meshComm,
 
             const lduMatrix& mat,
@@ -200,14 +187,22 @@ class GAMGSolver
             PtrList<FieldField<Field, scalar> >& otherIntCoeffs,
             PtrList<lduInterfaceFieldPtrsList>& otherInterfaces
         ) const;
-
+        void gatherMatrices
-        //- Gather field in allMatrix form
-        void gatherField
         (
-            const label meshComm,
             const labelList& procIDs,
-            const scalarField& coarsestSource,
+            const lduMesh& dummyMesh,
-            scalarField& allSource
+            const label meshComm,
+
+            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,
+            List<boolList>& otherTransforms,
+            List<List<int> >& otherRanks
         ) const;
 
         //-  Interpolate the correction after injected prolongation

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C

@@ -27,71 +27,10 @@ License
 #include "ICCG.H"
 #include "BICCG.H"
 #include "SubField.H"
-
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::GAMGSolver::gatherField
-(
-    const label meshComm,
-    const labelList& procIDs,
-    const scalarField& coarsestSource,
-    scalarField& allSource
-) const
-{
-    if (Pstream::myProcNo(meshComm) == procIDs[0])
-    {
-        // Master.
-        const lduMesh& allMesh = allMeshPtr_();
-
-        allSource.setSize(allMesh.lduAddr().size());
-
-        SubList<scalar>
-        (
-            allSource,
-            coarsestSource.size(),
-            0
-        ).assign(coarsestSource);
-
-        for (label i = 1; i < procIDs.size(); i++)
-        {
-            IPstream fromSlave
-            (
-                Pstream::scheduled,
-                procIDs[i],
-                0,          // bufSize
-                Pstream::msgType(),
-                meshComm
-            );
-
-            SubList<scalar> slaveSlots
-            (
-                allSource,
-                cellOffsets_[i+1]-cellOffsets_[i],
-                cellOffsets_[i]
-            );
-
-            UList<scalar>& l = slaveSlots;
-
-            fromSlave >> l;
-        }
-    }
-    else
-    {
-        // Send to master
-        OPstream toMaster
-        (
-            Pstream::scheduled,
-            procIDs[0],
-            0,
-            Pstream::msgType(),
-            meshComm
-        );
-        toMaster << coarsestSource;
-    }
-}
-
-
 Foam::solverPerformance Foam::GAMGSolver::solve
 (
     scalarField& psi,
@@ -519,7 +458,9 @@ void Foam::GAMGSolver::solveCoarsestLevel
         const List<int>& procIDs = UPstream::procID(coarseComm);
 
         scalarField allSource;
-        gatherField
+
+        const globalIndex cellOffsets(agglomeration_.cellOffsets());
+        cellOffsets.gather
         (
             coarseComm,
             procIDs,
@@ -578,55 +519,25 @@ void Foam::GAMGSolver::solveCoarsestLevel
                 }
                 UPstream::warnComm = oldWarn;
             }
-
-
-            // Distribute correction
-            coarsestCorrField = SubList<scalar>
-            (
-                allCorrField,
-                coarsestSource.size(),
-                0
-            );
-            for (label i=1; i < procIDs.size(); i++)
-            {
-                Pout<< "** master sending to " << procIDs[i] << endl;
-
-                OPstream toSlave
-                (
-                    Pstream::scheduled,
-                    procIDs[i],
-                    0,          // bufSize
-                    Pstream::msgType(),
-                    coarseComm
-                );
-
-                toSlave << SubList<scalar>
-                (
-                    allCorrField,
-                    cellOffsets_[i+1]-cellOffsets_[i],
-                    cellOffsets_[i]
-                );
-            }
-        }
-        else
-        {
-            Pout<< "** slave receiving from " << procIDs[0] << endl;
-            IPstream fromMaster
-            (
-                Pstream::scheduled,
-                procIDs[0],
-                0,          // bufSize
-                Pstream::msgType(),
-                coarseComm
-            );
-            fromMaster >> coarsestCorrField;
         }
 
+        // Scatter to all processors
+        coarsestCorrField.setSize(coarsestSource.size());
+        cellOffsets.scatter
+        (
+            coarseComm,
+            procIDs,
+            allCorrField,
+            coarsestCorrField
+        );
 
         if (debug >= 2)
         {
             coarseSolverPerf.print(Info(coarseComm));
         }
+
+        Pout<< "procAgglom: coarsestSource   :" << coarsestSource << endl;
+        Pout<< "procAgglom: coarsestCorrField:" << coarsestCorrField << endl;
     }
     else
     {
@@ -672,6 +583,9 @@ void Foam::GAMGSolver::solveCoarsestLevel
         {
             coarseSolverPerf.print(Info(coarseComm));
         }
+
+        Pout<< "GC: coarsestSource   :" << coarsestSource << endl;
+        Pout<< "GC: coarsestCorrField:" << coarsestCorrField << endl;
     }
 
     UPstream::warnComm = oldWarn;

src/OpenFOAM/meshes/lduMesh/lduMesh.C

@@ -91,6 +91,36 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const InfoProxy<lduMesh>& ip)
         }
     }
 
+
+    // Print actual contents
+    {
+        const labelList& l = addr.lowerAddr();
+        const labelList& u = addr.upperAddr();
+        forAll(l, faceI)
+        {
+            Pout<< "        face:" << faceI << " l:" << l[faceI]
+                << " u:" << u[faceI] << endl;
+        }
+        forAll(interfaces, i)
+        {
+            if (interfaces.set(i))
+            {
+                const labelUList& faceCells = addr.patchAddr(i);
+                if (faceCells.size())
+                {
+                    Pout<< "    patch:" << i
+                        << " type:" << interfaces[i].type() << endl;
+
+                    forAll(faceCells, i)
+                    {
+                        Pout<< "        " << i << " own:" << faceCells[i]
+                            << endl;
+                    }
+                }
+            }
+        }
+    }
+
     os.check("Ostream& operator<<(Ostream&, const lduMesh&");
 
     return os;

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C

@@ -217,22 +217,12 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
 }
 
 
-//Foam::lduPrimitiveMesh::lduPrimitiveMesh(Istream& is)
-//:
-//    lduAddressing(readLabel(is)),
-//    lowerAddr_(is),
-//    upperAddr_(is),
-//    patchAddr_(is),
-//    interfaces_(is),
-//    patchSchedule_(is),
-//    comm_(readLabel(is))
-//{}
-
 Foam::lduPrimitiveMesh::lduPrimitiveMesh
 (
     const label comm,
     const labelList& procIDs,
-    const PtrList<lduMesh>& procMeshes,
+    const lduMesh& myMesh,
+    const PtrList<lduMesh>& otherMeshes,
 
     labelList& cellOffsets,
     labelListList& faceMap,
@@ -240,7 +230,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     labelListListList& boundaryFaceMap
 )
 :
-    lduAddressing(size(procMeshes)),
+    lduAddressing(myMesh.lduAddr().size() + size(otherMeshes)),
     lowerAddr_(0),
     upperAddr_(0),
     patchAddr_(0),
@@ -261,40 +251,49 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                 << procIDs[i-1]
                 << exit(FatalError);
         }
-        if (procMeshes[i].comm() != procMeshes[i-1].comm())
+    }
+
+    const label currentComm = myMesh.comm();
+
+    forAll(otherMeshes, i)
+    {
+        if (otherMeshes[i].comm() != currentComm)
         {
             FatalErrorIn
             (
                 "lduPrimitiveMesh::lduPrimitiveMesh(..)"
-            )   << "Communicator " << procMeshes[i].comm() << " at index " << i
+            )   << "Communicator " << otherMeshes[i].comm()
+                << " at index " << i
                 << " differs from that of predecessor "
-                << procMeshes[i-1].comm()
+                << currentComm
                 << exit(FatalError);
         }
     }
 
-    const label currentComm = procMeshes[0].comm();
+    const label nMeshes = otherMeshes.size()+1;
 
     // Cells get added in order.
-    cellOffsets.setSize(procMeshes.size()+1);
+    cellOffsets.setSize(nMeshes+1);
     cellOffsets[0] = 0;
-    forAll(procMeshes, procMeshI)
+    for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
     {
-        const lduMesh& mesh = procMeshes[procMeshI];
+        const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
+
         cellOffsets[procMeshI+1] =
             cellOffsets[procMeshI]
-          + mesh.lduAddr().size();
+          + procMesh.lduAddr().size();
     }
 
     // Faces initially get added in order, sorted later
-    labelList internalFaceOffsets(procMeshes.size()+1);
+    labelList internalFaceOffsets(nMeshes+1);
     internalFaceOffsets[0] = 0;
-    forAll(procMeshes, procMeshI)
+    for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
     {
-        const lduMesh& mesh = procMeshes[procMeshI];
+        const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
+
         internalFaceOffsets[procMeshI+1] =
             internalFaceOffsets[procMeshI]
-          + mesh.lduAddr().lowerAddr().size();
+          + procMesh.lduAddr().lowerAddr().size();
     }
 
     // Count how faces get added. Interfaces inbetween get merged.
@@ -304,7 +303,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     // - interface in procMesh
     EdgeMap<labelPairList> mergedMap
     (
-        2*procMeshes[0].interfaces().size()
+        2*myMesh.interfaces().size()
     );
 
     // Unmerged interfaces: map from two processors back to
@@ -312,18 +311,18 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     // - interface in procMesh
     EdgeMap<labelPairList> unmergedMap(mergedMap.size());
 
-    boundaryMap.setSize(procMeshes.size());
+    boundaryMap.setSize(nMeshes);
-    boundaryFaceMap.setSize(procMeshes.size());
+    boundaryFaceMap.setSize(nMeshes);
 
 
     label nBoundaryFaces = 0;
     label nInterfaces = 0;
-    labelList nCoupledFaces(procMeshes.size(), 0);
+    labelList nCoupledFaces(nMeshes, 0);
 
-    forAll(procMeshes, procMeshI)
+    for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
     {
         const lduInterfacePtrsList interfaces =
-            procMeshes[procMeshI].interfaces();
+            mesh(myMesh, otherMeshes, procMeshI).interfaces();
 
         // Inialise all boundaries as merged
         boundaryMap[procMeshI].setSize(interfaces.size(), -1);
@@ -449,7 +448,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                 label procMeshI = elems[i][0];
                 label interfaceI = elems[i][1];
                 const lduInterfacePtrsList interfaces =
-                    procMeshes[procMeshI].interfaces();
+                    mesh(myMesh, otherMeshes, procMeshI).interfaces();
                 const processorLduInterface& pldui =
                     refCast<const processorLduInterface>
                     (
@@ -476,7 +475,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                 label procMeshI = elems[i][0];
                 label interfaceI = elems[i][1];
                 const lduInterfacePtrsList interfaces =
-                    procMeshes[procMeshI].interfaces();
+                    mesh(myMesh, otherMeshes, procMeshI).interfaces();
                 const processorLduInterface& pldui =
                     refCast<const processorLduInterface>
                     (
@@ -495,12 +494,13 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
 
 
     // Adapt faceOffsets for internal interfaces
-    labelList faceOffsets(procMeshes.size()+1);
+    labelList faceOffsets(nMeshes+1);
     faceOffsets[0] = 0;
-    faceMap.setSize(procMeshes.size());
+    faceMap.setSize(nMeshes);
-    forAll(procMeshes, procMeshI)
+    for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
     {
-        label nInternal = procMeshes[procMeshI].lduAddr().lowerAddr().size();
+        const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
+        label nInternal = procMesh.lduAddr().lowerAddr().size();
 
         faceOffsets[procMeshI+1] =
             faceOffsets[procMeshI]
@@ -524,10 +524,12 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     // Old internal faces and resolved coupled interfaces
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    forAll(procMeshes, procMeshI)
+    for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
     {
-        const labelUList& l = procMeshes[procMeshI].lduAddr().lowerAddr();
+        const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
-        const labelUList& u = procMeshes[procMeshI].lduAddr().upperAddr();
+
+        const labelUList& l = procMesh.lduAddr().lowerAddr();
+        const labelUList& u = procMesh.lduAddr().upperAddr();
 
         // Add internal faces
         label allFaceI = faceOffsets[procMeshI];
@@ -540,8 +542,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
         }
 
         // Add merged interfaces
-        const lduInterfacePtrsList interfaces =
+        const lduInterfacePtrsList interfaces = procMesh.interfaces();
-            procMeshes[procMeshI].interfaces();
 
         forAll(interfaces, intI)
         {
@@ -596,7 +597,13 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
 
 
                             const lduInterfacePtrsList nbrInterfaces =
-                                procMeshes[nbrProcMeshI].interfaces();
+                                mesh
+                                (
+                                    myMesh,
+                                    otherMeshes,
+                                    nbrProcMeshI
+                                ).interfaces();
+
 
                             const labelUList& faceCells =
                                 interfaces[intI].faceCells();
@@ -648,7 +655,13 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
             label procMeshI = elems[i][0];
             label interfaceI = elems[i][1];
             const lduInterfacePtrsList interfaces =
-                procMeshes[procMeshI].interfaces();
+                mesh
+                (
+                    myMesh,
+                    otherMeshes,
+                    procMeshI
+                ).interfaces();
+
             n += interfaces[interfaceI].faceCells().size();
         }
 
@@ -661,7 +674,12 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
             label procMeshI = elems[i][0];
             label interfaceI = elems[i][1];
             const lduInterfacePtrsList interfaces =
-                procMeshes[procMeshI].interfaces();
+                mesh
+                (
+                    myMesh,
+                    otherMeshes,
+                    procMeshI
+                ).interfaces();
 
             boundaryMap[procMeshI][interfaceI] = allInterfaceI;
             labelList& bfMap = boundaryFaceMap[procMeshI][interfaceI];
@@ -755,87 +773,39 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+const Foam::lduMesh& Foam::lduPrimitiveMesh::mesh
+(
+    const lduMesh& myMesh,
+    const PtrList<lduMesh>& otherMeshes,
+    const label meshI
+)
+{
+    return (meshI == 0 ? myMesh : otherMeshes[meshI-1]);
+}
+
+
 void Foam::lduPrimitiveMesh::gather
 (
     const lduMesh& mesh,
     const labelList& procIDs,
-    PtrList<lduMesh>& procMeshes
+    PtrList<lduMesh>& otherMeshes
 )
 {
+    // Force calculation of schedule (since does parallel comms)
+    (void)mesh.lduAddr().patchSchedule();
+
+
     const label meshComm = mesh.comm();
 
-    lduInterfacePtrsList interfaces(mesh.interfaces());
-
-    const lduAddressing& addressing = mesh.lduAddr();
-
     if (Pstream::myProcNo(meshComm) == procIDs[0])
     {
-        procMeshes.setSize(procIDs.size());
+        otherMeshes.setSize(procIDs.size()-1);
-
-        // Master mesh (copy for now. TBD)
-        {
-            // Convert to processorGAMGInterfaces in new communicator
-            lduInterfacePtrsList newInterfaces(interfaces.size());
-            labelListList newPatchAddr(interfaces.size());
-
-            forAll(interfaces, intI)
-            {
-                if (interfaces.set(intI))
-                {
-                    const processorLduInterface& pldui =
-                        refCast<const processorLduInterface>(interfaces[intI]);
-
-                    newPatchAddr[intI] = interfaces[intI].faceCells();
-                    labelList faceRestrictAddressing
-                    (
-                        identity(interfaces[intI].faceCells().size())
-                    );
-
-                    newInterfaces.set
-                    (
-                        intI,
-                        new processorGAMGInterface
-                        (
-                            intI,
-                            newInterfaces,
-                            newPatchAddr[intI],
-                            faceRestrictAddressing,
-                            pldui.comm(),
-                            pldui.myProcNo(),
-                            pldui.neighbProcNo(),
-                            pldui.forwardT(),
-                            pldui.tag()
-                        )
-                    );
-                }
-            }
-
-
-            procMeshes.set
-            (
-                0,
-                new lduPrimitiveMesh
-                (
-                    addressing.size(),
-                    addressing.lowerAddr(),
-                    addressing.upperAddr(),
-                    newPatchAddr,
-                    newInterfaces,
-                    nonBlockingSchedule<processorGAMGInterface>
-                    (
-                        newInterfaces
-                    ),
-                    meshComm
-                )
-            );
-        }
 
         // Slave meshes
         for (label i = 1; i < procIDs.size(); i++)
         {
             //Pout<< "on master :"
-            //    << " receiving from slave " << procIDs[i]
+            //    << " receiving from slave " << procIDs[i] << endl;
-            //    << endl;
 
             IPstream fromSlave
             (
@@ -849,46 +819,35 @@ void Foam::lduPrimitiveMesh::gather
             label nCells = readLabel(fromSlave);
             labelList lowerAddr(fromSlave);
             labelList upperAddr(fromSlave);
-            labelListList patchAddr(fromSlave);
+            boolList validInterface(fromSlave);
-            labelList comm(fromSlave);
-            labelList myProcNo(fromSlave);
-            labelList neighbProcNo(fromSlave);
-            List<tensorField> forwardT(fromSlave);
-            labelList tag(fromSlave);
 
-            // Convert to processorGAMGInterfaces
+            // Construct GAMGInterfaces
-            lduInterfacePtrsList newInterfaces(patchAddr.size());
+            lduInterfacePtrsList newInterfaces(validInterface.size());
-            forAll(patchAddr, intI)
+            labelListList patchAddr(validInterface.size());
+            forAll(validInterface, intI)
             {
-                if (tag[intI] != -1)
+                if (validInterface[intI])
                 {
-                    labelList faceRestrictAddressing
+                    word coupleType(fromSlave);
-                    (
-                        identity(patchAddr[intI].size())
-                    );
 
                     newInterfaces.set
                     (
                         intI,
-                        new processorGAMGInterface
+                        GAMGInterface::New
                         (
+                            coupleType,
                             intI,
                             newInterfaces,
-                            patchAddr[intI],
+                            fromSlave
-                            faceRestrictAddressing,
+                        ).ptr()
-                            comm[intI],
-                            myProcNo[intI],
-                            neighbProcNo[intI],
-                            forwardT[intI],
-                            tag[intI]
-                        )
                     );
+                    patchAddr[intI] = newInterfaces[intI].faceCells();
                 }
             }
 
-            procMeshes.set
+            otherMeshes.set
             (
-                i,
+                i-1,
                 new lduPrimitiveMesh
                 (
                     nCells,
@@ -900,7 +859,8 @@ void Foam::lduPrimitiveMesh::gather
                     (
                         newInterfaces
                     ),
-                    meshComm
+                    meshComm,
+                    true
                 )
             );
         }
@@ -909,28 +869,12 @@ void Foam::lduPrimitiveMesh::gather
     {
         // Send to master
 
-        //- Extract info from processorGAMGInterface
+        const lduAddressing& addressing = mesh.lduAddr();
-        labelListList patchAddr(interfaces.size());
+        lduInterfacePtrsList interfaces(mesh.interfaces());
-        labelList comm(interfaces.size(), -1);
+        boolList validInterface(interfaces.size());
-        labelList myProcNo(interfaces.size(), -1);
-        labelList neighbProcNo(interfaces.size(), -1);
-        List<tensorField> forwardT(interfaces.size());
-        labelList tag(interfaces.size(), -1);
-
         forAll(interfaces, intI)
         {
-            if (interfaces.set(intI))
+            validInterface[intI] = interfaces.set(intI);
-            {
-                const processorLduInterface& pldui =
-                    refCast<const processorLduInterface>(interfaces[intI]);
-
-                patchAddr[intI] = interfaces[intI].faceCells();
-                comm[intI] = pldui.comm();
-                myProcNo[intI] = pldui.myProcNo();
-                neighbProcNo[intI] =  pldui.neighbProcNo();
-                forwardT[intI] =  pldui.forwardT();
-                tag[intI] =  pldui.tag();
-            }
         }
 
         OPstream toMaster
@@ -941,38 +885,34 @@ void Foam::lduPrimitiveMesh::gather
             Pstream::msgType(),
             meshComm
         );
+
+        //Pout<< "sent nCells:" << addressing.size()
+        //    << "sent lowerAddr:" << addressing.lowerAddr().size()
+        //    << "sent upperAddr:" << addressing.upperAddr()
+        //    << "sent validInterface:" << validInterface
+        //    << endl;
+
         toMaster
             << addressing.size()
             << addressing.lowerAddr()
             << addressing.upperAddr()
-            << patchAddr
+            << validInterface;
-            << comm
+
-            << myProcNo
+        forAll(interfaces, intI)
-            << neighbProcNo
+        {
-            << forwardT
+            if (interfaces.set(intI))
-            << tag;
+            {
+                const GAMGInterface& interface = refCast<const GAMGInterface>
+                (
+                    interfaces[intI]
+                );
+
+                toMaster << interface.type();
+                interface.write(toMaster);
+            }
+        }
     }
 }
 
 
-// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
-
-//Foam::Ostream& Foam::operator<<
-//(
-//    Ostream& os,
-//    const Foam::lduPrimitiveMesh& mesh
-//)
-//{
-//    os  << mesh.size()
-//        << mesh.lowerAddr_
-//        << mesh.upperAddr_
-//        << mesh.patchAddr_
-//        << mesh.interfaces_
-//        << mesh.patchSchedule_
-//        << mesh.comm_
-//
-//    return os;
-//}
-
-
 // ************************************************************************* //

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H

@@ -43,12 +43,6 @@ SourceFiles
 namespace Foam
 {
 
-//// Forward declaration of friend functions and operators
-//
-//class lduPrimitiveMesh;
-//
-//Ostream& operator<<(Ostream&, const lduPrimitiveMesh&);
-
 /*---------------------------------------------------------------------------*\
                        Class lduPrimitiveMesh Declaration
 \*---------------------------------------------------------------------------*/
@@ -132,9 +126,6 @@ public:
             bool reUse
         );
 
-        ////- Construct from Istream
-        //lduPrimitiveMesh(Istream& is);
-
         //- Construct by combining multiple meshes. The meshes come from
         //  processors procIDs but these are only used to detect the interfaces
         //  that are to be merged.
@@ -149,13 +140,14 @@ public:
         //      merged boundaries: the new internal face (on the owner) or
         //                         -faceI-1 (on the neighbour)
         //
-        // TBD: hardcoded for processorGAMGInterface. Should use some 'clone' or
+        //// TBD: hardcoded for processorGAMGInterface. Should use some
-        // New functionality instead.
+        //// 'clone' or New functionality instead.
         lduPrimitiveMesh
         (
             const label comm,
             const labelList& procIDs,
-            const PtrList<lduMesh>& procMeshes,
+            const lduMesh& myMesh,
+            const PtrList<lduMesh>& otherMeshes,
 
             labelList& cellOffsets,
             labelListList& faceMap,
@@ -219,19 +211,22 @@ public:
 
         // Helper
 
-            //- Gather meshes from processors. Hardcoded for
+            //- Select either mesh0 (meshI is 0) or otherMeshes[meshI-1]
-            //  processorGAMGInterface
+            static const lduMesh& mesh
+            (
+                const lduMesh& mesh0,
+                const PtrList<lduMesh>& otherMeshes,
+                const label meshI
+            );
+
+            //- Gather meshes from other processors onto procIDs[0].
+            //  Hardcoded for GAMGInterface.
             static void gather
             (
                 const lduMesh& mesh,
                 const labelList& procIDs,
-                PtrList<lduMesh>& procMeshes
+                PtrList<lduMesh>& otherMeshes
             );
-
-//
-//    // Ostream operator
-//
-//        friend Ostream& operator<<(Ostream&, const procLduMatrix&);
 };