From a33f1787b7a960300563c2d10afa3ecbb16b8475 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 21 Mar 2019 17:55:05 +0000
Subject: [PATCH] ENH: globalIndex: helper function to get remote values

---
 .../test/globalIndex/Test-globalIndex.C       | 49 +++++++++++
 .../polyMesh/globalMeshData/globalIndex.C     | 60 +++++++++++++
 .../polyMesh/globalMeshData/globalIndex.H     | 27 ++++++
 .../globalMeshData/globalIndexTemplates.C     | 84 +++++++++++++++++++
 4 files changed, 220 insertions(+)

diff --git a/applications/test/globalIndex/Test-globalIndex.C b/applications/test/globalIndex/Test-globalIndex.C
index b33fb21306..1587a1352e 100644
--- a/applications/test/globalIndex/Test-globalIndex.C
+++ b/applications/test/globalIndex/Test-globalIndex.C
@@ -36,6 +36,7 @@ Description
 #include "Time.H"
 #include "polyMesh.H"
 #include "IOstreams.H"
+#include "Random.H"
 
 using namespace Foam;
 
@@ -163,6 +164,54 @@ int main(int argc, char *argv[])
         }
     }
 
+
+
+    // Get a few cell indices
+    const label nTotalCells = globalNumbering.size();
+
+    Random rndGen(Pstream::myProcNo());
+    DynamicList<label> globalIDs;
+    for (label i = 0; i < 100; i++)
+    {
+        globalIDs.append(rndGen.position(0, nTotalCells-1));
+    }
+
+    // Get the cell centres at those cell indices
+    List<point> ccs;
+    globalNumbering.get
+    (
+        ccs,
+        globalIDs,
+        [&mesh](List<point>& ccs, const labelUList& localIds)
+        {
+            ccs = UIndirectList<point>(mesh.cellCentres(), localIds);
+        }
+    );
+
+    // Do the brute-force method as well : collect all cell centres on all
+    // processors
+    pointField allCcs(globalNumbering.size());
+    globalNumbering.gather
+    (
+        Pstream::worldComm,
+        Pstream::procID(Pstream::worldComm),
+        mesh.cellCentres(),
+        allCcs
+    );
+    Pstream::scatter(allCcs);
+
+    // Compare
+    forAll(ccs, i)
+    {
+        const point& cc = ccs[i];
+        const point& allCc = allCcs[globalIDs[i]];
+        if (cc != allCc)
+        {
+            FatalErrorInFunction << "Problem cc:" << cc
+                << " allCc:" << allCc << exit(FatalError);
+        }
+    }
+
     return 0;
 }
 
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C
index 6bfe208782..914853d02f 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C
@@ -37,6 +37,66 @@ Foam::globalIndex::globalIndex(Istream& is)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::globalIndex::bin
+(
+    const labelUList& offsets,
+    const labelUList& globalIds,
+    labelList& order,
+    CompactListList<label>& bins,
+    DynamicList<label>& validBins
+)
+{
+    sortedOrder(globalIds, order);
+
+    bins.m() = UIndirectList<label>(globalIds, order);
+
+    labelList& binOffsets = bins.offsets();
+    binOffsets.setSize(offsets.size());
+    binOffsets = 0;
+
+    validBins.clear();
+
+    if (globalIds.size())
+    {
+        const label id = bins.m()[0];
+        label proci = findLower(offsets, id+1);
+
+        validBins.append(proci);
+        label binSize = 1;
+
+        for (label i = 1; i < order.size(); i++)
+        {
+            const label id = bins.m()[i];
+
+            if (id < offsets[proci+1])
+            {
+                binSize++;
+            }
+            else
+            {
+                // Not local. Reset proci
+                label oldProci = proci;
+                proci = findLower(offsets, id+1);
+
+                // Set offsets
+                for (label j = oldProci+1; j < proci; ++j)
+                {
+                    binOffsets[j] = binOffsets[oldProci]+binSize;
+                }
+                binOffsets[proci] = i;
+                validBins.append(proci);
+                binSize = 1;
+            }
+        }
+
+        for (label j = proci+1; j < binOffsets.size(); ++j)
+        {
+            binOffsets[j] = binOffsets[proci]+binSize;
+        }
+    }
+}
+
+
 void Foam::globalIndex::reset
 (
     const label localSize,
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H
index bc2f791fc3..96863c042b 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H
@@ -44,6 +44,8 @@ SourceFiles
 #define globalIndex_H
 
 #include "Pstream.H"
+#include "CompactListList.H"
+#include "DynamicList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,6 +67,19 @@ Ostream& operator<<(Ostream& os, const globalIndex& gi);
 
 class globalIndex
 {
+    // Private Member Functions
+
+        //- Sort and bin. validBins contains bins with non-zero size.
+        static void bin
+        (
+            const labelUList& offsets,
+            const labelUList& globalIds,
+            labelList& order,
+            CompactListList<label>& sortedElems,
+            DynamicList<label>& validBins
+        );
+
+
     // Private data
 
         //- Start of proci. Size is nProcs()+1. (so like CompactListList)
@@ -363,6 +378,18 @@ public:
                     Pstream::commsTypes::nonBlocking
             ) const;
 
+            //- Get (potentially remote) data. Elements required given as
+            //  global indices
+            template<class Type, class CombineOp>
+            void get
+            (
+                List<Type>& allFld,
+                const labelUList& globalIds,
+                const CombineOp& cop,
+                const label comm = Pstream::worldComm,
+                const int tag = UPstream::msgType()
+            ) const;
+
 
     // IOstream Operators
 
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexTemplates.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexTemplates.C
index 8f3e3a9763..67eb12d02d 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexTemplates.C
@@ -461,4 +461,88 @@ void Foam::globalIndex::scatter
 }
 
 
+template<class Type, class CombineOp>
+void Foam::globalIndex::get
+(
+    List<Type>& allFld,
+    const labelUList& globalIds,
+    const CombineOp& cop,
+    const label comm,
+    const int tag
+) const
+{
+    allFld.setSize(globalIds.size());
+    if (globalIds.size())
+    {
+        // Sort according to processor
+        labelList order;
+        CompactListList<label> bins;
+        DynamicList<label> validBins(Pstream::nProcs());
+        bin
+        (
+            offsets(),
+            globalIds,
+            order,
+            bins,
+            validBins
+        );
+
+        // Send local indices to individual processors as local index
+        PstreamBuffers sendBufs(Pstream::commsTypes::nonBlocking, tag, comm);
+
+        for (const auto proci : validBins)
+        {
+            const labelUList& es = bins[proci];
+
+            labelList localIDs(es.size());
+            forAll(es, i)
+            {
+                localIDs[i] = toLocal(proci, es[i]);
+            }
+
+            UOPstream os(proci, sendBufs);
+            os << localIDs;
+        }
+        labelList recvSizes;
+        sendBufs.finishedSends(recvSizes);
+
+
+        PstreamBuffers returnBufs(Pstream::commsTypes::nonBlocking, tag, comm);
+
+        forAll(recvSizes, proci)
+        {
+            if (recvSizes[proci])
+            {
+                UIPstream is(proci, sendBufs);
+                labelList localIDs(is);
+
+                // Collect entries
+                List<Type> fld(localIDs.size());
+                cop(fld, localIDs);
+
+                UOPstream os(proci, returnBufs);
+                os << fld;
+            }
+        }
+        returnBufs.finishedSends();
+
+        // Slot back
+        for (const auto proci : validBins)
+        {
+            label start = bins.offsets()[proci];
+            const SubList<label> es
+            (
+                order,
+                bins.offsets()[proci+1]-start,  // start
+                start
+            );
+            UIPstream is(proci, returnBufs);
+            List<Type> fld(is);
+
+            UIndirectList<Type>(allFld, es) = fld;
+        }
+    }
+}
+
+
 // ************************************************************************* //