Processor and GGI Block selective AMG: reference point

src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C

@@ -1034,13 +1034,15 @@ Foam::BlockMatrixSelection<Type>::restrictMatrix() const
             const lduInterface& fineInterface =
                 interfaceFields[intI].coupledInterface();
 
+            // Use filtered prolongation on master and slave side
+            // HJ, 7/Jun/2017
             coarseInterfaces.set
             (
                 intI,
                 SAMGInterface::New
                 (
                     coarseAddrPtr(),
-                    P,
+                    ownInterfaceProlongation[intI],
                     coarseInterfaces,
                     fineInterface,
                     nbrInterfaceProlongation[intI]

src/foam/matrices/lduMatrix/solvers/AMG/interfaces/AMGInterfaces/ggiAMGInterface/ggiAMGInterface.C

@@ -271,7 +271,7 @@ Foam::ggiAMGInterface::ggiAMGInterface
         return;
     }
 
-    // Continuing only with interfaces within the GGI comm only.
+    // Continuing with interfaces within the GGI comm only.
     // Note: on interfaces without the GGI comm, zone size will be zero
     // HJ, 11/Oct/2016
 

src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/SAMGInterface/SAMGInterface.H

@@ -63,10 +63,10 @@ class SAMGInterface
         //- Reference to ldu addressing
         const lduPrimitiveMesh& lduMesh_;
 
-        //- Reference to prolongation matrix
+        //- Reference to filtered (interface only) prolongation matrix
-        const crMatrix& prolongation_;
+        const crMatrix& interfaceProlongation_;
 
-        //- Reference to prolongation matrix
+        //- Reference to filtered neighbour prolongation matrix
         const crMatrix nbrInterfaceProlongation_;
 
 
@@ -125,14 +125,14 @@ public:
             lduInterface,
             (
                 const lduPrimitiveMesh& lduMesh,
-                const crMatrix& prolongation,
+                const crMatrix& interfaceProlongation,
                 const lduInterfacePtrsList& coarseInterfaces,
                 const lduInterface& fineInterface,
                 const crMatrix& nbrInterfaceProlongation
             ),
             (
                 lduMesh,
-                prolongation,
+                interfaceProlongation,
                 coarseInterfaces,
                 fineInterface,
                 nbrInterfaceProlongation
@@ -147,7 +147,7 @@ public:
         static autoPtr<SAMGInterface> New
         (
             const lduPrimitiveMesh& lduMesh,
-            const crMatrix& prolongation,
+            const crMatrix& interfaceProlongation,
             const lduInterfacePtrsList& coarseInterfaces,
             const lduInterface& fineInterface,
             const crMatrix& nbrInterfaceProlongation
@@ -161,12 +161,12 @@ public:
         SAMGInterface
         (
             const lduPrimitiveMesh& lduMesh,
-            const crMatrix& prolongation,
+            const crMatrix& interfaceProlongation,
             const crMatrix& nbrInterfaceProlongation
         )
         :
             lduMesh_(lduMesh),
-            prolongation_(prolongation),
+            interfaceProlongation_(interfaceProlongation),
             nbrInterfaceProlongation_(nbrInterfaceProlongation)
         {}
 
@@ -187,9 +187,9 @@ public:
             }
 
             //- Return reference to filtered prolongation matrix
-            const crMatrix& prolongation() const
+            const crMatrix& interfaceProlongation() const
             {
-                return prolongation_;
+                return interfaceProlongation_;
             }
 
             //- Return reference to neighbour's filtered prolongation matrix

src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/SAMGInterface/newSAMGInterface.C

@@ -32,7 +32,7 @@ License
 Foam::autoPtr<Foam::SAMGInterface> Foam::SAMGInterface::New
 (
     const lduPrimitiveMesh& lduMesh,
-    const crMatrix& prolongation,
+    const crMatrix& interfaceProlongation,
     const lduInterfacePtrsList& coarseInterfaces,
     const lduInterface& fineInterface,
     const crMatrix& nbrInterfaceProlongation
@@ -50,7 +50,7 @@ Foam::autoPtr<Foam::SAMGInterface> Foam::SAMGInterface::New
             "SAMGInterface::New\n"
             "(\n"
             "    const lduPrimitiveMesh& lduMesh,\n"
-            "    const crMatrix& prolongation,\n"
+            "    const crMatrix& interfaceProlongation,\n"
             "    const lduInterfacePtrsList& coarseInterfaces,\n"
             "    const lduInterface& fineInterface,\n"
             "    const crMatrix& nbrInterfaceProlongation\n"
@@ -66,7 +66,7 @@ Foam::autoPtr<Foam::SAMGInterface> Foam::SAMGInterface::New
         cstrIter()
         (
             lduMesh,
-            prolongation,
+            interfaceProlongation,
             coarseInterfaces,
             fineInterface,
             nbrInterfaceProlongation

src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/ggiSAMGInterface/ggiSAMGInterface.C

@@ -227,13 +227,13 @@ void Foam::ggiSAMGInterface::initFastReduce() const
 Foam::ggiSAMGInterface::ggiSAMGInterface
 (
     const lduPrimitiveMesh& lduMesh,
-    const crMatrix& prolongation,
+    const crMatrix& interfaceProlongation,
     const lduInterfacePtrsList& coarseInterfaces,
     const lduInterface& fineInterface,
     const crMatrix& nbrInterfaceProlongation
 )
 :
-    SAMGInterface(lduMesh, prolongation, nbrInterfaceProlongation),
+    SAMGInterface(lduMesh, interfaceProlongation, nbrInterfaceProlongation),
     fineGgiInterface_(refCast<const ggiLduInterface>(fineInterface)),
     zoneSize_(0),
     zoneAddressing_(),
@@ -275,7 +275,7 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
         return;
     }
 
-    // Continuing only with interfaces within the GGI comm only.
+    // Continuing with interfaces within the GGI comm only.
     // Note: on interfaces without the GGI comm, zone size will be zero
     // HJ, 11/Oct/2016
 
@@ -298,6 +298,10 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
         );
     }
 
+    // 
+    const crAddressing& prolongationCr = interfaceProlongation.crAddr();
+    const crAddressing& nbrExpandCr = nbrExpandProlongation.crAddr();
+
     // Create addressing for neighbour processors.  Note: expandAddrToZone will
     // expand the addressing to zone size.  HJ, 13/Jun/2016
     labelField neighbourExpandProc
@@ -379,11 +383,6 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
         label curMasterProc, curSlaveProc, curSide, nbrSide;
         long curMaster, curSlave;
 
-        // Save weights and columns of restriction and prolongation for master
-        // and slave
-        scalarField restrictionWeights, prolongationWeights;
-        labelList restrictionCol, prolongationCol;
-
         // ZONE
         forAll (fineZa, fineZaI)
         {
@@ -439,31 +438,36 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
                 if (fineGgiInterface_.master())
                 {
                     curSide = fineZaI;
+                    curMasterProc = Pstream::myProcNo();
                     nbrSide = nnI;
+                    curSlaveProc = neighbourExpandProc[nnI];
                 }
                 // Slave
                 else
                 {
                     curSide = nnI;
+                    curMasterProc = neighbourExpandProc[nnI];
                     nbrSide = fineZaI;
+                    curSlaveProc = Pstream::myProcNo();
                 }
 
                 // Triple product for only one row of prolongation and
                 // restriction - with included weights from GGI
                 for
                 (
-                    label indexR = prolongation.crAddr().rowStart()[curSide];
+                    label indexR = prolongationCr.rowStart()[curSide];
-                    indexR < prolongation.crAddr().rowStart()[curSide + 1];
+                    indexR < prolongationCr.rowStart()[curSide + 1];
                     indexR++
                 )
                 {
                     // Grab weight from restriction
-                    scalar rWeight = prolongation.coeffs()[indexR];
+                    scalar rWeight = interfaceProlongation.coeffs()[indexR];
 
+                    // HJ, replace nbrInterfaceProlongation with nbrExpandProlongation
                     for
                     (
-                        label indexP = nbrInterfaceProlongation.crAddr().rowStart()[nbrSide];
+                        label indexP = nbrExpandCr.rowStart()[nbrSide];
-                        indexP < nbrInterfaceProlongation.crAddr().rowStart()[nbrSide + 1];
+                        indexP < nbrExpandCr.rowStart()[nbrSide + 1];
                         indexP++
                     )
                     {
@@ -473,10 +477,10 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
                         // Code in the current master and slave - used for
                         // identifying the face
                         curMaster =
-                            prolongation.crAddr().column()[indexP]
+                            prolongationCr.column()[indexR]
                           + procOffset*curMasterProc;
                         curSlave =
-                            nbrInterfaceProlongation.crAddr().column()[indexP]
+                            nbrExpandCr.column()[indexP]
                           + procOffset*curSlaveProc;
 
                         if (neighboursTable.found(curMaster))
@@ -498,7 +502,8 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
                             DynamicList<DynamicList<label, 4>, 4>& curFaceFaces =
                                 faceFaceTable.find(curMaster)();
 
-                            DynamicList<DynamicList<label, 4>, 4>& curFaceFaceNbrs =
+                            DynamicList<DynamicList<label, 4>, 4>&
+                                curFaceFaceNbrs =
                                 faceFaceNbrTable.find(curMaster)();
 
                             DynamicList<DynamicList<scalar, 4>, 4>& curFaceWeights =
@@ -604,14 +609,234 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
                 }
             } // end for all current neighbours
         } // end for all fine faces
+        if (fineGgiInterface_.master())
+        {
+            Info<< "MASTER: ";
+        }
+        else
+        {
+            Info<< "SLAVE: ";
+        }
+        Info<< "Done fine level, 1. Created " << nAgglomPairs << " pairs and "
+            << nCoarseFaces << " faces" << endl;
     }
 //------------------------------------------------------------------------------
-//                               RUBBISH FROM HERE
+//                          FINE LEVEL - no GGI weights!
 //------------------------------------------------------------------------------
     else
     {
-        Info << "Not yet implementedi Copy from above and remove weights from GGI" << endl;
+        // Perform analysis only for local faces
+        // HJ, 22/Jun/2016
 
+        label curMasterProc, curSlaveProc, curSide, nbrSide;
+        long curMaster, curSlave;
+
+        // ZONE
+        forAll (fineZa, fineZaI)
+        {
+            // Get the local face (from zone) to analyse
+            const label ffI = fineZa[fineZaI];
+
+            curMaster = -1;
+            curMasterProc = -1;
+            curSlave = -1;
+            curSlaveProc = -1;
+
+            // Note.  Signalling in global clustering requires
+            // me to recognise clustering from separate
+            // processors as separate.  In the first phase,
+            // this will be used to recognise cluster from
+            // each processor as separate and in the second
+            // phase it will be used to filter local processor
+            // faces from the global patch.  Currently, I am
+            // calculating unique cluster index as:
+            //
+            // id = cluster + procOffset*myProcID
+            //
+            // With procOffset = 1 million, this should be
+            // sufficient for 2000 CPUs with 2 million coarse
+            // cells each.  For larger numbers, I need a
+            // larger max int, which can be changed on request
+            // HJ, 1/Apr/2009
+
+            // For this face - take the column of restriction (row of
+            // prolongation) which corresponds to the face, every element of
+            // this column will multiply the boundary coefficient, then take
+            // neighbour's row of prolongation which coressponds to the
+            // face, multiply with each of the weights and save
+
+            if (fineGgiInterface_.master())
+            {
+                curSide = fineZaI;
+                curMasterProc = Pstream::myProcNo();
+                nbrSide = ffI;
+                curSlaveProc = neighbourExpandProc[ffI];
+            }
+            // Slave
+            else
+            {
+                curSide = ffI;
+                curMasterProc = neighbourExpandProc[ffI];
+                nbrSide = fineZaI;
+                curSlaveProc = Pstream::myProcNo();
+            }
+
+            // Triple product for only one row of prolongation and
+            // restriction - with included weights from GGI
+            for
+            (
+                label indexR = prolongationCr.rowStart()[curSide];
+                indexR < prolongationCr.rowStart()[curSide + 1];
+                indexR++
+            )
+            {
+                // Grab weight from restriction
+                scalar rWeight = interfaceProlongation.coeffs()[indexR];
+
+                // HJ, replace nbrInterfaceProlongation with nbrExpandProlongation
+                for
+                (
+                    label indexP = nbrExpandCr.rowStart()[nbrSide];
+                    indexP < nbrExpandCr.rowStart()[nbrSide + 1];
+                    indexP++
+                )
+                {
+                    // Grab weight from prolongation
+                    scalar pWeight = nbrInterfaceProlongation.coeffs()[indexP];
+
+                    // Code in the current master and slave - used for
+                    // identifying the face
+                    curMaster = prolongationCr.column()[indexR]
+                        + procOffset*curMasterProc;
+                    curSlave = nbrExpandCr.column()[indexP]
+                        + procOffset*curSlaveProc;
+
+                    if (neighboursTable.found(curMaster))
+                    {
+                        // This contribution already exists - add the result
+                        // to the existing contribution
+
+                        // This master side face already exists
+
+                        // Check all current neighbours to see if the current
+                        // slave already exists.  If so, add the coefficient.
+
+                        DynamicList<long, 4>& curNbrs =
+                            neighboursTable.find(curMaster)();
+
+                        DynamicList<label, 4>& curNbrsProc =
+                            nbrsProcTable.find(curMaster)();
+
+                        DynamicList<DynamicList<label, 4>, 4>& curFaceFaces =
+                            faceFaceTable.find(curMaster)();
+
+                        DynamicList<DynamicList<label, 4>, 4>& curFaceFaceNbrs =
+                            faceFaceNbrTable.find(curMaster)();
+
+                        DynamicList<DynamicList<scalar, 4>, 4>& curFaceWeights =
+                            faceFaceWeightsTable.find(curMaster)();
+
+                        // Search for coded neighbour
+                        bool nbrFound = false;
+
+                        forAll (curNbrs, curNbrI)
+                        {
+                            // Check neighbour slave
+                            if (curNbrs[curNbrI] == curSlave)
+                            {
+                                nbrFound = true;
+                                curFaceFaces[curNbrI].append(ffI);
+                                curFaceFaceNbrs[curNbrI].append(0);
+                                curFaceWeights[curNbrI].append(pWeight*rWeight);
+
+                                // New agglomeration pair found in already
+                                // existing pair
+                                nAgglomPairs++;
+
+                                break;
+                            }
+                        }
+
+                        if (!nbrFound)
+                        {
+                            curNbrs.append(curSlave);
+                            curNbrsProc.append(curSlaveProc);
+
+                            DynamicList<label, 4> newFF;
+                            newFF.append(ffI);
+                            curFaceFaces.append(newFF);
+
+                            DynamicList<label, 4> newFNbr;
+                            newFNbr.append(0);
+                            curFaceFaceNbrs.append(newFNbr);
+
+                            DynamicList<scalar, 4> newFW;
+                            newFW.append(pWeight*rWeight);
+                            curFaceWeights.append(newFW);
+
+                            // New coarse face created for an existing master
+                            nCoarseFaces++;
+                            nAgglomPairs++;
+                        }
+                    }
+                    else
+                    {
+                        // This master has got no neighbours yet.
+                        // Add a neighbour, proc and a coefficient as a
+                        // new list, thus creating a new face
+
+                        DynamicList<long, 4> newNbrs;
+                        newNbrs.append(curSlave);
+                        neighboursTable.insert
+                        (
+                            curMaster,
+                            newNbrs
+                        );
+
+                        DynamicList<label, 4> newNbrsProc;
+                        newNbrsProc.append(curSlaveProc);
+                        nbrsProcTable.insert
+                        (
+                            curMaster,
+                            newNbrsProc
+                        );
+
+                        DynamicList<DynamicList<label, 4>, 4> newFF;
+                        newFF.append(DynamicList<label, 4>());
+                        newFF[0].append(ffI);
+                        faceFaceTable.insert
+                        (
+                            curMaster,
+                            newFF
+                        );
+
+                        DynamicList<DynamicList<label, 4>, 4> newFNbr;
+                        newFNbr.append(DynamicList<label, 4>());
+                        newFNbr[0].append(0);
+                        faceFaceNbrTable.insert
+                        (
+                            curMaster,
+                            newFNbr
+                        );
+
+                        DynamicList<DynamicList<scalar, 4>, 4> newFFWeights;
+                        newFFWeights.append(DynamicList<scalar, 4>());
+                        newFFWeights[0].append(pWeight*rWeight);
+                        faceFaceWeightsTable.insert
+                        (
+                            curMaster,
+                            newFFWeights
+                        );
+
+                        // New coarse face created for a new master
+                        nCoarseFaces++;
+                        nAgglomPairs++;
+                    }
+                }
+            }
+        } // end for all fine faces
+        Info<< "Done coarse level, 1. Created " << nAgglomPairs << " pairs and "
+            << nCoarseFaces << " faces" << endl;
     } // end of else in fine level (coarse level)
 
     // Since only local faces are analysed, lists can now be resized
@@ -669,7 +894,7 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
     // Sort makes sure the order is identical on both sides.
     // HJ, 20/Feb/2009 and 6/Jun/2016
     sort(contents);
-
+    Info<< "START MATRIX ASSEMBLY" << endl;
     // Note: Restriction is done on master side only because this is where
     // the local zone is created.  HJ, 1/Aug/2016
     if (master())
@@ -758,6 +983,8 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
                 nProcFaces++;
             }
         }
+        Info<< "MASTER ASSEMBLY: Created " << nAgglomPairs << " pairs and "
+            << nProcFaces << " master faces" << endl;
 
         // No need to resize arrays only local faces are used
         // HJ, 1/Aug/2016
@@ -910,7 +1137,11 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
                 nProcFaces++;
             }
         }
+        Info<< "SLAVE ASSEMBLY: Created " << nAgglomPairs << " pairs and "
+            << nProcFaces << " slave faces" << endl;
+
     }
+    Info<< "END MATRIX ASSEMBLY" << endl;
 }
 
 
@@ -1106,6 +1337,8 @@ void Foam::ggiSAMGInterface::expandAddrToZone(labelField& lf) const
 
 void Foam::ggiSAMGInterface::expandCrMatrixToZone(crMatrix&) const
 {
+    notImplemented("expandCrMatrixToZone");
+    // Code missing: collapse crMatrices into a zone crMatrix
     if (!localParallel())
     {
         notImplemented("expandCrMatrixToZone");

src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/ggiSAMGInterface/ggiSAMGInterface.H

@@ -127,7 +127,7 @@ public:
         ggiSAMGInterface
         (
             const lduPrimitiveMesh& lduMesh,
-            const crMatrix& prolongation,
+            const crMatrix& interfaceProlongation,
             const lduInterfacePtrsList& coarseInterfaces,
             const lduInterface& fineInterface,
             const crMatrix& nbrInterfaceProlongation
@@ -230,7 +230,7 @@ public:
             virtual void initProlongationTransfer
             (
                 const Pstream::commsTypes commsType,
-                const crMatrix& P
+                const crMatrix& filteredP
             ) const;
 
             //- Transfer and return prolongation matrix adjacent to
@@ -238,7 +238,7 @@ public:
             virtual autoPtr<crMatrix> prolongationTransfer
             (
                 const Pstream::commsTypes commsType,
-                const crMatrix& P
+                const crMatrix& filteredP
             ) const;
 
 

src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/processorSAMGInterface/processorSAMGInterface.C

@@ -45,71 +45,83 @@ namespace Foam
 Foam::processorSAMGInterface::processorSAMGInterface
 (
     const lduPrimitiveMesh& lduMesh,
-    const crMatrix& prolongation,
+    const crMatrix& interfaceProlongation,
     const lduInterfacePtrsList& coarseInterfaces,
     const lduInterface& fineInterface,
     const crMatrix& nbrInterfaceProlongation
 )
 :
-    SAMGInterface(lduMesh, prolongation, nbrInterfaceProlongation),
+    SAMGInterface(lduMesh, interfaceProlongation, nbrInterfaceProlongation),
     fineProcInterface_(refCast<const processorLduInterface>(fineInterface)),
     comm_(fineProcInterface_.comm()),
     tag_(fineProcInterface_.tag())
 {
-    /*
+    // MASTER processor
-    // On master:
+    if (myProcNo() < neighbProcNo())
-    //  - receive slave's filtered prolongation
+    {
-    //  - create a transpose of your local prolongation
+        // READ FIRST
-    //  - do the triple product: store the coefficients and the addresses
+//------------------------------------------------------------------------------
-    //    row = row of restriction, column = column of prolongation
+        // This code is written for the FILTERED local matrix. Why?
-    //  - the coefficients should be stored row by row, but with randomized
+        // Because the addressing is very natural. If this is not ok, filter the
-    //    columns
+        // big matrix here using fineInterface.faceCells(), transpose, and use
-    //  - send the addressing to slave to sort its coefficients to match
+        // the same names for crMatrix arrays (then, there is no need to change
-    //    the ordering on master (is it possible to sort at the same time?)
+        // the triple product code).
+//------------------------------------------------------------------------------
 
-    // Needed:
+        // Algorithm details:
-    //  - labelList masterOwner,
+        // Go through the triple product similar to SAMG Policy - collect
-    //  - labelList masterNeighbour,
+        // coarse addressing on master - send it to slave to sort its addressing
-    //  - label nCoarseCoeffs
+        // accordingly - HashTable containing labelPair(owner, neighbour) as key
+        // and coarse face index as entry
 
+        // The resulting contributions from triple product will be saved using
+        // the following arrays:
+        // - faceCells_: for coarse matrix, it tells us the owner of the
+        //               boundary coefficient (on local side)
+        // - fineAddressing_: saves the index of the fine boundary coefficient
+        //                    in the natural order of appearance on master in
+        //                    triple product
+        // - restrictAddressing_: saves the coarse face for which the
+        //                        contribution (coefficient) needs to sum-up
+        // - restrictWeights_: weights for each fine boundary coefficient in
+        //                     natural order of appearance on master
 
     // Question: Should we put the masterOwner and masterNeighbour into a linked
     // list which can be searched quickly? In this case we can create the
     // addressing on master, send and at the same time, calculate the coarse
     // coefficients on both sides and store into corresponding order
 
-    // Question: How can I save the addressing of the boundary coeffs for
+        // *Note: (A*B)^T = B^T*A^T, which is exactly what I have. Use the fact
-    // the next level?
+        // on the slave processor!!!
 
-    // Question: With such boundary coefficients, my vector-matrix product is a
+        // Addressing of coarse faces - for sorting faceCells_ on slave
-    // matrix-matrix product - change multiplication rule in SAMGInterfaceFields
+        HashTable<label, labelPair, Hash<labelPair> > coarseLevelTable;
-    // or account for it with addressing?
-//------------------------------------------------------------------------------
 
-    // Count coarse coefficients and create addressing on master
+        // Lists for saving addressing and weights - give size, but check at the
-    if (myProcNo() < neighbProcNo())
+        // end!
-    {
+        faceCells_.setSize(5*fineProcInterface_.interfaceSize());
-//------------------------------------------------------------------------------
+        fineAddressing_.setSize(5*fineProcInterface_.interfaceSize());
-//                              COUNT COARSE COEFFS
+        restrictWeights_.setSize(5*fineProcInterface_.interfaceSize());
-//------------------------------------------------------------------------------
+        restrictAddressing_.setSize(5*fineProcInterface_.interfaceSize());
-        // Filtered prolongation matrix from my side
-        const labelList& localOwner = prolongation.crAddr().column();
 
         // Filtered prolongation matrix from my side
-        tmp<crMatrix> tProlongationT(prolongation.T());
+        crMatrix prolongationT = interfaceProlongation.T();
-        const labelList& rRowStart = tProlongationT->crAddr().rowStart();
+        const labelList& rRowStart = prolongationT.crAddr().rowStart();
-        const labelList& rColumn = tProlongationT->crAddr().column();
+        const labelList& rColumn = prolongationT.crAddr().column();
-        const label rNRows = tProlongationT->crAddr().nRows();
+        const scalarField& rCoeffs = prolongationT.coeffs();
+        const label rNRows = prolongationT.crAddr().nRows();
 
         // Filtered prolongation matrix from neighbour - to obtain restriction,
         // make a transpose
         const labelList& pRowStart =
             nbrInterfaceProlongation.crAddr().rowStart();
         const labelList& pColumn = nbrInterfaceProlongation.crAddr().column();
+        const scalarField& pCoeffs = nbrInterfaceProlongation.coeffs();
         const label pNCols = nbrInterfaceProlongation.crAddr().nCols();
 
         labelList coeffMark(pNCols, -1);
         label nCoarseCoeffs = 0;
+        label nCoarseContribs = 0;
 
         // Row of R
         for (label ir = 0; ir < rNRows; ir++)
@@ -126,6 +138,9 @@ Foam::processorSAMGInterface::processorSAMGInterface
             )
             {
                 // Column of coefficient in R
+                // This is an important information: it tells us which
+                // boundary coefficient this restriction coeff multiplies! It is
+                // the index of the boundary coeff in the boundary array.
                 const label jr = rColumn[indexR];
 
                 for
@@ -135,114 +150,101 @@ Foam::processorSAMGInterface::processorSAMGInterface
                     indexP++
                 )
                 {
+                    // Column of P, goes into coarse address
                     label jp = pColumn[indexP];
 
-                    if (coeffMark[jp] == -1)
+                    // To which coeff do I add myself to? Does the address
-                    {
+                    // in this row already exist?
-                        // Found a new coarse coeff!
-                        coeffMark[jp] = nCoarseCoeffs;
-                        nCoarseCoeffs++;
-                    }
-                }
-            }
-        }
-
-//------------------------------------------------------------------------------
-//                      CREATE HASH TABLE FOR ADDRESSING
-//------------------------------------------------------------------------------
-
-        HashTable<label, labelPair, Hash<labelPair> > coarseLevel(nCoarseCoeffs);
-        faceCells_.setSize(nCoarseCoeffs);
-
-        // Reset for creating addressing
-        nCoarseCoeffs = 0;
-        coeffMark = -1;
-
-        // Row of R
-        for (label ir = 0; ir < rNRows; ir++)
-        {
-            // Restart coeffMark for each restriction row
-            coeffMark = -1;
-
-            // Row start addressing R
-            for
-            (
-                label indexR = rRowStart[ir];
-                indexR < rRowStart[ir + 1];
-                indexR++
-            )
-            {
-                // Column of coefficient in R
-                const label jr = rColumn[indexR];
-
-                // Grab column of the original prolongation, this is the
-                // faceCell_!
-                const label ja = localOwner[ir];
-
-                for
-                (
-                    label indexP = pRowStart[jr];
-                    indexP < pRowStart[jr + 1];
-                    indexP++
-                )
-                {
-                    label jp = pColumn[indexP];
-
                     label identify = coeffMark[jp];
 
+                    // If the coeff at this address doesn't exist
                     if (identify == -1)
                     {
                         // Found a new coarse coeff!
                         identify = nCoarseCoeffs;
                         coeffMark[jp] = nCoarseCoeffs;
 
-                        faceCells_[identify] = ja;
+                        faceCells_[identify] = ir;
 
-                        // Save address and coeff into HashTable
+                        // Save address and coeff into HashTable for sorting
-                        coarseLevel.insert(labelPair(ir, jp), nCoarseCoeffs);
+                        // faceCells on slave
+                        coarseLevelTable.insert
+                        (
+                            labelPair(ir, jp),
+                            nCoarseCoeffs
+                        );
                         nCoarseCoeffs++;
                     }
+
+                    restrictWeights_[nCoarseContribs] =
+                        rCoeffs[indexR]*pCoeffs[indexP];
+                    fineAddressing_[nCoarseContribs] = jr;
+                    restrictAddressing_[nCoarseContribs] = identify;
+
+                    nCoarseContribs++;
                 }
             }
         }
 
+        // Check for fixed lists
+        if (nCoarseContribs > 5*fineProcInterface_.interfaceSize())
+        {
+            FatalErrorIn("processorSAMGInterface::processorSAMGInterface(...)")
+                << "Coarse SAMG processor siginificantly bigger than fine: "
+                << "nCoarseFaces = " << nCoarseContribs
+                << " nFineFaces = " << fineProcInterface_.interfaceSize()
+                << abort(FatalError);
+        }
+
+        // Resize arrays to final size
+        faceCells_.setSize(nCoarseCoeffs);
+        fineAddressing_.setSize(nCoarseContribs);
+        restrictWeights_.setSize(nCoarseContribs);
+        restrictAddressing_.setSize(nCoarseContribs);
+
         // Send to slave
         OPstream toNbr(Pstream::blocking, neighbProcNo());
 
-        toNbr<< coarseLevel;
+        toNbr<< coarseLevelTable << nCoarseContribs;
-
     }
-    // Slave
+    // Slave processor
     else
     {
         // Slave must receive the hash table sent by the master to know how to
-        // sort its boundary coefficients
+        // sort its faceCells_
-
         IPstream fromNbr(Pstream::blocking, neighbProcNo());
 
-        masterCoarseLevel_ =
+        HashTable<label, labelPair, Hash<labelPair> > masterCoarseLevel =
             HashTable<label, labelPair, Hash<labelPair> >(fromNbr);
 
-        // Filtered prolongation matrix from my side
+        const label nCoarseContribs = readLabel(fromNbr);
-        const labelList& localOwner = prolongation.crAddr().column();
 
         // Filtered prolongation matrix from my side
-        tmp<crMatrix> tProlongationT(prolongation.T());
+        crMatrix prolongationT = interfaceProlongation.T();
-        const labelList& rRowStart = tProlongationT->crAddr().rowStart();
+        const labelList& rRowStart = prolongationT.crAddr().rowStart();
-        const labelList& rColumn = tProlongationT->crAddr().column();
+        const labelList& rColumn = prolongationT.crAddr().column();
-        const label rNRows = tProlongationT->crAddr().nRows();
+        const scalarField& rCoeffs = prolongationT.coeffs();
+        const label rNRows = prolongationT.crAddr().nRows();
 
         // Filtered prolongation matrix from neighbour - to obtain restriction,
         // make
         // a transpose
-        const labelList& pRowStart = nbrInterfaceProlongation.crAddr().rowStart();
+        const labelList& pRowStart =
+            nbrInterfaceProlongation.crAddr().rowStart();
         const labelList& pColumn = nbrInterfaceProlongation.crAddr().column();
+        const scalarField& pCoeffs = nbrInterfaceProlongation.coeffs();
         const label pNCols = nbrInterfaceProlongation.crAddr().nCols();
 
-        faceCells_.setSize(masterCoarseLevel_.size());
+        faceCells_.setSize(masterCoarseLevel.size(), -1);
+        restrictWeights_.setSize(nCoarseContribs);
+        fineAddressing_.setSize(nCoarseContribs);
+        restrictAddressing_.setSize(nCoarseContribs);
 
         labelList coeffMark(pNCols, -1);
+        label nCoarseEntries = 0;
 
+        // This loop is only for sorting the faceCells_ on slave side to match
+        // the order on master
         // Row of R
         for (label ir = 0; ir < rNRows; ir++)
         {
@@ -259,11 +261,6 @@ Foam::processorSAMGInterface::processorSAMGInterface
             {
                 // Column of coefficient in R
                 const label jr = rColumn[indexR];
-
-                // Grab column of the original prolongation -  this is the
-                // faceCell_!!!
-                label ja = localOwner[ir];
-
                 for
                 (
                     label indexP = pRowStart[jr];
@@ -271,306 +268,34 @@ Foam::processorSAMGInterface::processorSAMGInterface
                     indexP++
                 )
                 {
+                    // Column of P, into coarse address
                     label jp = pColumn[indexP];
 
+                    // Array for marking new contributions in the row ir
                     label identify = coeffMark[jp];
 
                     if (identify == -1)
                     {
-                        label address = masterCoarseLevel_[labelPair(jp,ja)];
+                        label address = masterCoarseLevel[labelPair(jp, ir)];
 
                         // Found a new coarse coeff!
                         identify = address;
 
                         coeffMark[jp] = address;
-                        faceCells_[identify] = ja;
+                        faceCells_[identify] = ir;
-                    }
-                }
-            }
-        }
                     }
 
+                    restrictWeights_[nCoarseEntries] =
+                        rCoeffs[indexR]*pCoeffs[indexP];
+                    fineAddressing_[nCoarseEntries] = jr;
+                    restrictAddressing_[nCoarseEntries] = identify;
 
-    // Analyse the local and neighbour row label:
+                    nCoarseEntries++;
-    //  local coarse, remote coarse = regular coarse face
-    //  local coarse, remote fine = local expanded face: receive prolonged data
-    //  local fine, remote coarse = neighbour expanded face: send prolonged data
-    //  local fine, remote fine = local and neighbour expanded face: send and
-    //                            receive
-
-    // Algorithm
-    // Go through local row labels and examine cases
-    // 1) coarse local and coarse remote:
-    //    create coarse processor face and set faceCells.  Set weight to 1
-    // 2) coarse local and fine remote:
-    //    will receive coarse neighbours and weights from opposite side
-    // 3) fine local and coarse remote:
-    //    assemble local coarse neighbours and weights from the prolongation
-    // 4) fine local and fine remote - ignore
-    //
-    // On completion of selection:
-    //    send and receive local coarse neighbours for fine local/coarse remote
-    //    sort coarse equations by increasing coarse index from master side
-    //    create faceCells, fineAddressing and fineWeights
-
-    // First collect and communicate internal neighbours from the local fine
-    // side (ie neighbour processor cell is coarse)
-    HashTable<labelList, label, Hash<label> > neighboursFromLocalFine
-    (
-        Foam::max(128, fineProcInterface_.interfaceSize()/4)
-    );
-
-    HashTable<scalarField, label, Hash<label> > weightsFromLocalFine
-    (
-        Foam::max(128, fineProcInterface_.interfaceSize()/4)
-    );
-
-    // For FINE to FINE communication - local boundary coefficient multiplyed by
-    // local weight - send to the other side for triple product
-    HashTable<labelList, label, Hash<label> > neighboursFromLocalFineToFine
-    (
-        Foam::max(128, fineProcInterface_.interfaceSize()/4)
-    );
-
-    HashTable<scalarField, label, Hash<label> > weightsFromLocalFineToFine
-    (
-        Foam::max(128, fineProcInterface_.interfaceSize()/4)
-    );
-
-
-    // Get access to the prolongation addressing and coefficients
-    const labelList& pRowStart = prolongation.crAddr().rowStart();
-    const labelList& pColumn = prolongation.crAddr().column();
-    const scalarField& pCoeffs = prolongation.coeffs();
-
-    // Get fine faceCells
-    const labelList& fineFaceCells = fineInterface.faceCells();
-
-    // Collect local fine to neighbour coarse connections for communication
-    // TU, 19 May 2017: and fine to fine (communication through triple product)
-    forAll (localRowLabel, faceI)
-    {
-        if
-        (
-            localRowLabel[faceI] < 0
-         && neighbourRowLabel[faceI] >= 0
-        )
-        {
-            // Found local fine to neighbour coarse interface
-
-            // Collect local coarse neighbours and weights from the
-            // prolongation matrix to send to other processor
-            const label curStart = pRowStart[fineFaceCells[faceI]];
-            const label curEnd = pRowStart[fineFaceCells[faceI] + 1];
-            const label nCoarse = curEnd - curStart;
-
-            labelList nbrs(nCoarse);
-            scalarField weights(nCoarse);
-
-            forAll (nbrs, i)
-            {
-                nbrs[i] = pColumn[curStart + i];
-                weights[i] = pCoeffs[curStart + i];
-            }
-
-            // Insert neighbours under remote coarse index
-            neighboursFromLocalFine.insert(neighbourRowLabel[faceI], nbrs);
-            weightsFromLocalFine.insert(neighbourRowLabel[faceI], weights);
-        }
-        // FINE to FINE communication - TU, May 2017
-        if
-        (
-            localRowLabel[faceI] < 0
-         && neighbourRowLabel[faceI] < 0
-        )
-        {
-            // Found local FINE to neighbour FINE interface
-
-            // First, multiply all local interface boundary coefficients with
-            // local prolongation weights - you will send this to the other side
-            // in localBoundaryProlongation
-
-            // Collect local coarse neighbours and weights from the
-            // prolongation matrix to send to other processor
-            const label curStart = pRowStart[fineFaceCells[faceI]];
-            const label curEnd = pRowStart[fineFaceCells[faceI] + 1];
-            const label nCoeffs = curEnd - curStart;
-
-            labelList nbrs(nCoeffs);
-            scalarField weights(nCoeffs);
-
-            forAll (nbrs, i)
-            {
-                nbrs[i] = pColumn[curStart + i];
-                weights[i] = pCoeffs[curStart + i];
-            }
-
-            // Insert neighbours under remote coarse index
-            neighboursFromLocalFineToFine.insert(neighbourRowLabel[faceI], nbrs);
-            weightsFromLocalFineToFine.insert(neighbourRowLabel[faceI], weights);
-        }
-    }
-
-    // Receive remote prolongation data
-    HashTable<labelList, label, Hash<label> > neighboursFromRemoteFine;
-    HashTable<scalarField, label, Hash<label> > weightsFromRemoteFine;
-    HashTable<labelList, label, Hash<label> > neighboursFromRemoteFineToFine;
-    HashTable<scalarField, label, Hash<label> > weightsFromRemoteFineToFine;
-
-    // Send and receive the addressing from the other side
-    {
-        OPstream toNbr(Pstream::blocking, neighbProcNo());
-        toNbr<< neighboursFromLocalFine << weightsFromLocalFine;
-    }
-
-    {
-        IPstream fromNbr(Pstream::blocking, neighbProcNo());
-
-        neighboursFromRemoteFine =
-            HashTable<labelList, label, Hash<label> >(fromNbr);
-
-        weightsFromRemoteFine =
-            HashTable<scalarField, label, Hash<label> >(fromNbr);
-
-        neighboursFromRemoteFineToFine =
-            HashTable<labelList, label, Hash<label> >(fromNbr);
-
-        weightsFromRemoteFineToFine =
-            HashTable<scalarField, label, Hash<label> >(fromNbr);
-
-    }
-
-    // Assemble connectivity
-
-    // Resize arrays to size of fine interface
-    // Note: it is technically possible to have MORE coarse processor faces
-    // so a check will be made in the end
-    faceCells_.setSize(5*fineProcInterface_.interfaceSize());
-    fineAddressing_.setSize(5*fineProcInterface_.interfaceSize());
-    fineWeights_.setSize(5*fineProcInterface_.interfaceSize());
-
-    // Count coarse faces
-    label nCoarseFaces = 0;
-
-    // Collect coarse-to-fine connections
-    forAll (localRowLabel, faceI)
-    {
-        if
-        (
-            localRowLabel[faceI] >= 0
-         && neighbourRowLabel[faceI] >= 0
-        )
-        {
-            // Found local coarse to neighbour coarse face
-
-            // Create new coarse face
-            faceCells_[nCoarseFaces] = localRowLabel[faceI];
-            fineAddressing_[nCoarseFaces] = faceI;
-            fineWeights_[nCoarseFaces] = 1;
-            nCoarseFaces++;
-        }
-        else if
-        (
-            localRowLabel[faceI] < 0
-         && neighbourRowLabel[faceI] >= 0
-        )
-        {
-            // Found local fine to neighbour coarse face
-
-            // Pick up local prolongation coarse entries and for all
-            // add a new face with appropriate weight
-            // Note: faceCells changes due to (internal) local coarse cells
-            const labelList& curLocalCoarseNbrs =
-                neighboursFromLocalFine[neighbourRowLabel[faceI]];
-
-            const scalarField& curLocalCoarseWeights =
-                weightsFromLocalFine[neighbourRowLabel[faceI]];
-
-            forAll (curLocalCoarseNbrs, curNbrI)
-            {
-                // Create new coarse face
-                faceCells_[nCoarseFaces] = curLocalCoarseNbrs[curNbrI];
-                fineAddressing_[nCoarseFaces] = faceI;
-                fineWeights_[nCoarseFaces] = curLocalCoarseWeights[curNbrI];
-                nCoarseFaces++;
-            }
-        }
-        else if
-        (
-            localRowLabel[faceI] >= 0
-         && neighbourRowLabel[faceI] < 0
-        )
-        {
-            // Found local coarse to neighbour fine face
-
-            // Pick up neighbour prolongation coarse entries and for all
-            // add a new face with appropriate weight
-            // Note: faceCells remains the same: single local coarse cell
-            const labelList& curNbrCoarseNbrs =
-                neighboursFromRemoteFine[localRowLabel[faceI]];
-
-            const scalarField& curNbrCoarseWeights =
-                weightsFromRemoteFine[localRowLabel[faceI]];
-
-            forAll (curNbrCoarseNbrs, curNbrI)
-            {
-                // Create new coarse face
-                faceCells_[nCoarseFaces] = localRowLabel[faceI];
-                fineAddressing_[nCoarseFaces] = faceI;
-                fineWeights_[nCoarseFaces] = curNbrCoarseWeights[curNbrI];
-                nCoarseFaces++;
-            }
-        }
-        else
-        {
-            // Fine to fine.
-            // Establish connection as in triple product: multiply restriction
-            // matrix on my side, send to the other side and multiply with
-            // prolongation there
-
-            // Get fine faceCells
-            const scalarField& fineLocalCoeffs =
-            // Pick up local prolongation coarse entries and for all
-            // add a new face with appropriate weight
-            // Note: faceCells changes, but how???
-
-            const labelList& curLocalCoarseNbrs =
-                neighboursFromLocalFineToFine[neighbourRowLabel[faceI]];
-
-            const scalarField& curLocalCoarseWeights =
-                weightsFromLocalFineToFine[neighbourRowLabel[faceI]];
-
-            forAll (curLocalCoarseNbrs, curNbrI)
-            {
-                // Create new coarse face
-                faceCells_[nCoarseFaces] = curLocalCoarseNbrs[curNbrI];
-                fineAddressing_[nCoarseFaces] = faceI;
-                fineWeights_[nCoarseFaces] =
-                    curLocalCoarseWeights[curNbrI]*fineLocalCoeffs[faceI];
-                nCoarseFaces++;
-
-                // Store fineFaceCells, recording the fine index of prolonged
-                // cell locally and all coarse neighbours on the other side
-                // HJ, HERE!!!
                 }
             }
         }
 
-    // Check for fixed lists
-    if (nCoarseFaces > 5*fineProcInterface_.interfaceSize())
-    {
-        FatalErrorIn("processorSAMGInterface::processorSAMGInterface(...)")
-            << "Coarse SAMG processor siginificantly bigger than fine: "
-            << "nCoarseFaces = " << nCoarseFaces
-            << " nFineFaces = " << fineProcInterface_.interfaceSize()
-            << abort(FatalError);
     }
-
-    // Resize arrays to final size
-    faceCells_.setSize(nCoarseFaces);
-    fineAddressing_.setSize(nCoarseFaces);
-    fineWeights_.setSize(nCoarseFaces);
-*/
 }
 
 

src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/processorSAMGInterface/processorSAMGInterface.H

@@ -87,7 +87,7 @@ public:
         processorSAMGInterface
         (
             const lduPrimitiveMesh& lduMesh,
-            const crMatrix& prolongation,
+            const crMatrix& interfaceProlongation,
             const lduInterfacePtrsList& coarseInterfaces,
             const lduInterface& fineInterface,
             const crMatrix& nbrInterfaceProlongation
@@ -165,7 +165,7 @@ public:
             virtual void initProlongationTransfer
             (
                 const Pstream::commsTypes commsType,
-                const crMatrix& P
+                const crMatrix& filteredP
             ) const;
 
             //- Transfer and return prolongation matrix adjacent to
@@ -173,7 +173,7 @@ public:
             virtual autoPtr<crMatrix> prolongationTransfer
             (
                 const Pstream::commsTypes commsType,
-                const crMatrix& P
+                const crMatrix& filteredP
             ) const;