Merge remote-tracking branch 'origin/feature/procAgglom'

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

@@ -126,8 +126,7 @@ void Foam::GAMGAgglomeration::restrictField
         const List<int>& procIDs = agglomProcIDs(coarseLevelIndex);
         const labelList& offsets = cellOffsets(coarseLevelIndex);
 
-        const globalIndex cellOffsetter(offsets);
+        globalIndex::gather(offsets, fineComm, procIDs, cf);
-        cellOffsetter.gather(fineComm, procIDs, cf);
     }
 }
 
@@ -192,19 +191,10 @@ void Foam::GAMGAgglomeration::prolongField
         const List<int>& procIDs = agglomProcIDs(coarseLevelIndex);
         const labelList& offsets = cellOffsets(coarseLevelIndex);
 
-        const globalIndex cellOffsetter(offsets);
-
         label localSize = nCells_[levelIndex];
 
-
         Field<Type> allCf(localSize);
-        cellOffsetter.scatter
+        globalIndex::scatter(offsets, coarseComm, procIDs, cf, allCf);
-        (
-            coarseComm,
-            procIDs,
-            cf,
-            allCf
-        );
 
         forAll(fineToCoarse, i)
         {

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/GAMGProcAgglomeration/GAMGProcAgglomeration.C

@@ -25,6 +25,7 @@ License
 
 #include "GAMGProcAgglomeration.H"
 #include "GAMGAgglomeration.H"
+#include "lduMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -114,6 +115,152 @@ void Foam::GAMGProcAgglomeration::printStats
 }
 
 
+Foam::labelListList Foam::GAMGProcAgglomeration::globalCellCells
+(
+    const lduMesh& mesh
+)
+{
+    const lduAddressing& addr = mesh.lduAddr();
+    lduInterfacePtrsList interfaces = mesh.interfaces();
+
+    const label myProcID = Pstream::myProcNo(mesh.comm());
+
+    globalIndex globalNumbering
+    (
+        addr.size(),
+        Pstream::msgType(),
+        mesh.comm(),
+        Pstream::parRun()
+    );
+
+    labelList globalIndices(addr.size());
+    forAll(globalIndices, cellI)
+    {
+        globalIndices[cellI] = globalNumbering.toGlobal(myProcID, cellI);
+    }
+
+
+    // Get the interface cells
+    PtrList<labelList> nbrGlobalCells(interfaces.size());
+    {
+        // Initialise transfer of restrict addressing on the interface
+        forAll(interfaces, inti)
+        {
+            if (interfaces.set(inti))
+            {
+                interfaces[inti].initInternalFieldTransfer
+                (
+                    Pstream::nonBlocking,
+                    globalIndices
+                );
+            }
+        }
+
+        if (Pstream::parRun())
+        {
+            Pstream::waitRequests();
+        }
+
+        forAll(interfaces, inti)
+        {
+            if (interfaces.set(inti))
+            {
+                nbrGlobalCells.set
+                (
+                    inti,
+                    new labelList
+                    (
+                        interfaces[inti].internalFieldTransfer
+                        (
+                            Pstream::nonBlocking,
+                            globalIndices
+                        )
+                    )
+                );
+            }
+        }
+    }
+
+
+    // Scan the neighbour list to find out how many times the cell
+    // appears as a neighbour of the face. Done this way to avoid guessing
+    // and resizing list
+    labelList nNbrs(addr.size(), 1);
+
+    const labelUList& nbr = addr.upperAddr();
+    const labelUList& own = addr.lowerAddr();
+
+    {
+        forAll(nbr, faceI)
+        {
+            nNbrs[nbr[faceI]]++;
+            nNbrs[own[faceI]]++;
+        }
+
+        forAll(interfaces, inti)
+        {
+            if (interfaces.set(inti))
+            {
+                const labelUList& faceCells = interfaces[inti].faceCells();
+
+                forAll(faceCells, i)
+                {
+                    nNbrs[faceCells[i]]++;
+                }
+            }
+        }
+    }
+
+
+    // Create cell-cells addressing
+    labelListList cellCells(addr.size());
+
+    forAll(cellCells, cellI)
+    {
+        cellCells[cellI].setSize(nNbrs[cellI], -1);
+    }
+
+    // Reset the list of number of neighbours to zero
+    nNbrs = 0;
+
+    // Scatter the neighbour faces
+    forAll(nbr, faceI)
+    {
+        label c0 = own[faceI];
+        label c1 = nbr[faceI];
+
+        cellCells[c0][nNbrs[c0]++] = globalIndices[c1];
+        cellCells[c1][nNbrs[c1]++] = globalIndices[c0];
+    }
+    forAll(interfaces, inti)
+    {
+        if (interfaces.set(inti))
+        {
+            const labelUList& faceCells = interfaces[inti].faceCells();
+
+            forAll(faceCells, i)
+            {
+                label c0 = faceCells[i];
+                cellCells[c0][nNbrs[c0]++] = nbrGlobalCells[inti][i];
+            }
+        }
+    }
+
+    forAll(cellCells, cellI)
+    {
+        Foam::stableSort(cellCells[cellI]);
+    }
+
+    // Replace the initial element (always -1) with the local cell
+    forAll(cellCells, cellI)
+    {
+        cellCells[cellI][0] = globalIndices[cellI];
+    }
+
+    return cellCells;
+}
+
+
 bool Foam::GAMGProcAgglomeration::agglomerate
 (
     const label fineLevelIndex,
@@ -243,7 +390,4 @@ Foam::GAMGProcAgglomeration::~GAMGProcAgglomeration()
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/GAMGProcAgglomeration/GAMGProcAgglomeration.H

@@ -44,6 +44,7 @@ namespace Foam
 {
 
 class GAMGAgglomeration;
+class lduMesh;
 
 /*---------------------------------------------------------------------------*\
                     Class GAMGProcAgglomeration Declaration
@@ -74,6 +75,9 @@ protected:
             const label procAgglomComm
         );
 
+        //- Debug: calculate global cell-cells
+        static labelListList globalCellCells(const lduMesh&);
+
 private:
 
     // Private data

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/processorGAMGInterface/processorGAMGInterface.H

@@ -29,7 +29,6 @@ Description
 
 SourceFiles
     processorGAMGInterface.C
-    processorGAMGInterfaceTemplates.C
 
 \*---------------------------------------------------------------------------*/
 

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C

@@ -360,6 +360,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
           + procMesh.lduAddr().size();
     }
 
+
     // Faces initially get added in order, sorted later
     labelList internalFaceOffsets(nMeshes+1);
     internalFaceOffsets[0] = 0;
@@ -818,9 +819,37 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     {
         const labelPairList& elems = iter();
 
-        // Sort processors in increasing order
+        // Sort processors in increasing order so both sides walk through in
-        labelList order(identity(elems.size()));
+        // same order.
-        Foam::sort(order, procLess(elems));
+        labelPairList procPairs(elems.size());
+        forAll(elems, i)
+        {
+            const labelPair& elem = elems[i];
+            label procMeshI = elem[0];
+            label interfaceI = elem[1];
+            const lduInterfacePtrsList interfaces = mesh
+            (
+                myMesh,
+                otherMeshes,
+                procMeshI
+            ).interfaces();
+
+            const processorLduInterface& pldui =
+                refCast<const processorLduInterface>
+                (
+                    interfaces[interfaceI]
+                );
+            label myProcNo = pldui.myProcNo();
+            label nbrProcNo = pldui.neighbProcNo();
+            procPairs[i] = labelPair
+            (
+                min(myProcNo, nbrProcNo),
+                max(myProcNo, nbrProcNo)
+            );
+        }
+
+        labelList order;
+        sortedOrder(procPairs, order);
 
         // Count
         label n = 0;
@@ -837,11 +866,6 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                 procMeshI
             ).interfaces();
 
-            //Pout<< "    adding interface:" << " procMeshI:" << procMeshI
-            //    << " interface:" << interfaceI
-            //    << " size:" << interfaces[interfaceI].faceCells().size()
-            //    << endl;
-
             n += interfaces[interfaceI].faceCells().size();
         }
 
@@ -868,6 +892,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
 
             const labelUList& l = interfaces[interfaceI].faceCells();
             bfMap.setSize(l.size());
+
             forAll(l, faceI)
             {
                 allFaceCells[n] = cellOffsets[procMeshI]+l[faceI];
@@ -947,7 +972,10 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
 
     patchSchedule_ = nonBlockingSchedule<processorGAMGInterface>(interfaces_);
 
+    if (debug)
+    {
         checkUpperTriangular(cellOffsets.last(), lowerAddr_, upperAddr_);
+    }
 }
 
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H

@@ -150,6 +150,19 @@ public:
 
         // Other
 
+            //- Collect data in processor order on master (== procIDs[0]).
+            //  Needs offsets only on master.
+            template<class Type>
+            static void gather
+            (
+                const labelUList& offsets,
+                const label comm,
+                const labelList& procIDs,
+                const UList<Type>& fld,
+                List<Type>& allFld,
+                const int tag = UPstream::msgType()
+            );
+
             //- Collect data in processor order on master (== procIDs[0]).
             //  Needs offsets only on master.
             template<class Type>
@@ -160,7 +173,22 @@ public:
                 const UList<Type>& fld,
                 List<Type>& allFld,
                 const int tag = UPstream::msgType()
-            ) const;
+            ) const
+            {
+                gather(offsets_, comm, procIDs, fld, allFld, tag);
+            }
+
+            //- Inplace collect data in processor order on master
+            //  (== procIDs[0]). Needs offsets only on master.
+            template<class Type>
+            static void gather
+            (
+                const labelUList& offsets,
+                const label comm,
+                const labelList& procIDs,
+                List<Type>& fld,
+                const int tag = UPstream::msgType()
+            );
 
             //- Inplace collect data in processor order on master
             //  (== procIDs[0]). Needs offsets only on master.
@@ -171,7 +199,22 @@ public:
                 const labelList& procIDs,
                 List<Type>& fld,
                 const int tag = UPstream::msgType()
-            ) const;
+            ) const
+            {
+                gather(offsets_, comm, procIDs, fld, tag);
+            }
+
+            //- Distribute data in processor order. Requires fld to be sized!
+            template<class Type>
+            static void scatter
+            (
+                const labelUList& offsets,
+                const label comm,
+                const labelList& procIDs,
+                const UList<Type>& allFld,
+                UList<Type>& fld,
+                const int tag = UPstream::msgType()
+            );
 
             //- Distribute data in processor order. Requires fld to be sized!
             template<class Type>
@@ -182,7 +225,10 @@ public:
                 const UList<Type>& allFld,
                 UList<Type>& fld,
                 const int tag = UPstream::msgType()
-            ) const;
+            ) const
+            {
+                scatter(offsets_, comm, procIDs, allFld, fld, tag);
+            }
 
 
     // IOstream Operators

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexTemplates.C

@@ -30,28 +30,24 @@ License
 template<class Type>
 void Foam::globalIndex::gather
 (
+    const labelUList& offsets,
     const label comm,
     const labelList& procIDs,
     const UList<Type>& fld,
     List<Type>& allFld,
     const int tag
-) const
+)
 {
     if (Pstream::myProcNo(comm) == procIDs[0])
     {
-        allFld.setSize(size());
+        allFld.setSize(offsets.last());
 
         // Assign my local data
         SubList<Type>(allFld, fld.size(), 0).assign(fld);
 
         for (label i = 1; i < procIDs.size(); i++)
         {
-            SubList<Type> procSlot
+            SubList<Type> procSlot(allFld, offsets[i+1]-offsets[i], offsets[i]);
-            (
-                allFld,
-                offsets_[i+1]-offsets_[i],
-                offsets_[i]
-            );
 
             if (contiguous<Type>())
             {
@@ -112,27 +108,23 @@ void Foam::globalIndex::gather
 template<class Type>
 void Foam::globalIndex::gather
 (
+    const labelUList& offsets,
     const label comm,
     const labelList& procIDs,
     List<Type>& fld,
     const int tag
-) const
+)
 {
     if (Pstream::myProcNo(comm) == procIDs[0])
     {
-        List<Type> allFld(size());
+        List<Type> allFld(offsets.last());
 
         // Assign my local data
         SubList<Type>(allFld, fld.size(), 0).assign(fld);
 
         for (label i = 1; i < procIDs.size(); i++)
         {
-            SubList<Type> procSlot
+            SubList<Type> procSlot(allFld, offsets[i+1]-offsets[i], offsets[i]);
-            (
-                allFld,
-                offsets_[i+1]-offsets_[i],
-                offsets_[i]
-            );
 
             if (contiguous<Type>())
             {
@@ -195,31 +187,25 @@ void Foam::globalIndex::gather
 template<class Type>
 void Foam::globalIndex::scatter
 (
+    const labelUList& offsets,
     const label comm,
     const labelList& procIDs,
     const UList<Type>& allFld,
     UList<Type>& fld,
     const int tag
-) const
+)
 {
     if (Pstream::myProcNo(comm) == procIDs[0])
     {
-        fld.assign
+        fld.assign(SubList<Type>(allFld, offsets[1]-offsets[0]));
-        (
-            SubList<Type>
-            (
-                allFld,
-                offsets_[1]-offsets_[0]
-            )
-        );
 
         for (label i = 1; i < procIDs.size(); i++)
         {
             const SubList<Type> procSlot
             (
                 allFld,
-                offsets_[i+1]-offsets_[i],
+                offsets[i+1]-offsets[i],
-                offsets_[i]
+                offsets[i]
             );
 
             if (contiguous<Type>())