From 4528eb8a69582bf9007db7296d1a337ffa39aa73 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Wed, 8 Jun 2011 17:26:32 +0100
Subject: [PATCH] ENH: backgroundMeshDecompostion distributed queries.

---
 .../backgroundMeshDecomposition.C             | 418 +++++++++++++++++-
 .../backgroundMeshDecomposition.H             |  38 +-
 .../conformalVoronoiMesh.C                    |  50 ++-
 .../conformalVoronoiMesh.H                    |  16 +-
 .../conformalVoronoiMeshConformToSurface.C    |   5 +-
 5 files changed, 512 insertions(+), 15 deletions(-)

diff --git a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index 41a37cbb17..2e512fbe8c 100644
--- a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -353,6 +353,8 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
         cellWeights.write();
         mesh_.write();
     }
+
+    buildPatchAndTree();
 }
 
 
@@ -616,6 +618,12 @@ void Foam::backgroundMeshDecomposition::buildPatchAndTree()
         )
     );
 
+    // Give the bounds of every processor to every other processor
+    allBackgroundMeshBounds_[Pstream::myProcNo()] = overallBb;
+
+    Pstream::gatherList(allBackgroundMeshBounds_);
+    Pstream::scatterList(allBackgroundMeshBounds_);
+
     if (debug)
     {
         OFstream fStr
@@ -659,6 +667,92 @@ void Foam::backgroundMeshDecomposition::buildPatchAndTree()
 }
 
 
+Foam::autoPtr<Foam::mapDistribute> Foam::backgroundMeshDecomposition::buildMap
+(
+    const List<label>& toProc
+) const
+{
+    // Determine send map
+    // ~~~~~~~~~~~~~~~~~~
+
+    // 1. Count
+    labelList nSend(Pstream::nProcs(), 0);
+
+    forAll(toProc, i)
+    {
+        label procI = toProc[i];
+
+        nSend[procI]++;
+    }
+
+    // Send over how many I need to receive
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelListList sendSizes(Pstream::nProcs());
+
+    sendSizes[Pstream::myProcNo()] = nSend;
+
+    combineReduce(sendSizes, UPstream::listEq());
+
+    // 2. Size sendMap
+    labelListList sendMap(Pstream::nProcs());
+
+    forAll(nSend, procI)
+    {
+        sendMap[procI].setSize(nSend[procI]);
+
+        nSend[procI] = 0;
+    }
+
+    // 3. Fill sendMap
+    forAll(toProc, i)
+    {
+        label procI = toProc[i];
+
+        sendMap[procI][nSend[procI]++] = i;
+    }
+
+    // Determine receive map
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    labelListList constructMap(Pstream::nProcs());
+
+    // Local transfers first
+    constructMap[Pstream::myProcNo()] = identity
+    (
+        sendMap[Pstream::myProcNo()].size()
+    );
+
+    label constructSize = constructMap[Pstream::myProcNo()].size();
+
+    forAll(constructMap, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            label nRecv = sendSizes[procI][Pstream::myProcNo()];
+
+            constructMap[procI].setSize(nRecv);
+
+            for (label i = 0; i < nRecv; i++)
+            {
+                constructMap[procI][i] = constructSize++;
+            }
+        }
+    }
+
+    return autoPtr<mapDistribute>
+    (
+        new mapDistribute
+        (
+            constructSize,
+            sendMap.xfer(),
+            constructMap.xfer()
+        )
+    );
+}
+
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
@@ -687,6 +781,7 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
     ),
     boundaryFacesPtr_(),
     bFTreePtr_(),
+    allBackgroundMeshBounds_(Pstream::nProcs()),
     decomposeDict_
     (
         IOobject
@@ -738,8 +833,6 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
     Info<< nl << "Building initial background mesh decomposition" << endl;
 
     initialRefinement();
-
-    buildPatchAndTree();
 }
 
 
@@ -963,7 +1056,7 @@ Foam::backgroundMeshDecomposition::distribute
 }
 
 
-bool Foam::backgroundMeshDecomposition::positionInThisDomain
+bool Foam::backgroundMeshDecomposition::positionOnThisProcessor
 (
     const point& pt
 ) const
@@ -974,4 +1067,323 @@ bool Foam::backgroundMeshDecomposition::positionInThisDomain
 }
 
 
+Foam::labelList Foam::backgroundMeshDecomposition::processorPosition
+(
+    const pointField& pts
+) const
+{
+    DynamicList<label> toCandidateProc;
+    DynamicList<point> testPoints;
+    labelList ptBlockStart(pts.size(), -1);
+    labelList ptBlockSize(pts.size(), -1);
+
+    label nTotalCandidates = 0;
+
+    forAll(pts, pI)
+    {
+        const point& pt = pts[pI];
+
+        label nCandidates = 0;
+
+        forAll(allBackgroundMeshBounds_, procI)
+        {
+            if (allBackgroundMeshBounds_[procI].contains(pt))
+            {
+                toCandidateProc.append(procI);
+                testPoints.append(pt);
+
+                nCandidates++;
+            }
+        }
+
+        ptBlockStart[pI] = nTotalCandidates;
+        ptBlockSize[pI] = nCandidates;
+
+        nTotalCandidates += nCandidates;
+    }
+
+    // Needed for reverseDistribute
+    label preDistributionToCandidateProcSize = toCandidateProc.size();
+
+    autoPtr<mapDistribute> map(buildMap(toCandidateProc));
+
+    map().distribute(testPoints);
+
+    List<bool> pointOnCandidate(testPoints.size(), false);
+
+    // Test candidate points on candidate processors
+    forAll(testPoints, tPI)
+    {
+        pointOnCandidate[tPI] = positionOnThisProcessor(testPoints[tPI]);
+    }
+
+    map().reverseDistribute
+    (
+        preDistributionToCandidateProcSize,
+        pointOnCandidate
+    );
+
+    labelList ptProc(pts.size(), -1);
+
+    forAll(pts, pI)
+    {
+        // Extract the sub list of results for this point
+
+        SubList<bool> ptProcResults
+        (
+            pointOnCandidate,
+            ptBlockSize[pI],
+            ptBlockStart[pI]
+        );
+
+        forAll(ptProcResults, pPRI)
+        {
+            if (ptProcResults[pPRI])
+            {
+                ptProc[pI] = toCandidateProc[ptBlockStart[pI] + pPRI];
+
+                break;
+            }
+        }
+
+        if (ptProc[pI] < 0)
+        {
+            // No processor was found
+            FatalErrorIn
+            (
+                "Foam::labelList"
+                "Foam::backgroundMeshDecomposition::processorPosition"
+                "("
+                    "const pointField& pts"
+                ") const"
+            )
+                << "The position " << pts[pI]
+                << " was not found in any part of the background mesh."
+                << exit(FatalError);
+        }
+    }
+
+    return ptProc;
+}
+
+
+Foam::List<Foam::List<Foam::pointIndexHit> >
+Foam::backgroundMeshDecomposition::intersectsProcessors
+(
+    const pointField& starts,
+    const pointField& ends,
+    bool includeOwnProcessor
+) const
+{
+    DynamicList<label> toCandidateProc;
+    DynamicList<point> testStarts;
+    DynamicList<point> testEnds;
+    labelList segmentBlockStart(starts.size(), -1);
+    labelList segmentBlockSize(starts.size(), -1);
+
+    label nTotalCandidates = 0;
+
+    forAll(starts, sI)
+    {
+        const point& s = starts[sI];
+        const point& e = ends[sI];
+
+        // Dummy point for treeBoundBox::intersects
+        point p(vector::zero);
+
+        label nCandidates = 0;
+
+        forAll(allBackgroundMeshBounds_, procI)
+        {
+            // It is assumed that the sphere in question overlaps the source
+            // processor, so don't test it, unless includeOwnProcessor is true
+            if
+            (
+                (includeOwnProcessor || procI != Pstream::myProcNo())
+              && allBackgroundMeshBounds_[procI].intersects(s, e, p)
+            )
+            {
+                toCandidateProc.append(procI);
+                testStarts.append(s);
+                testEnds.append(e);
+
+                nCandidates++;
+            }
+        }
+
+        segmentBlockStart[sI] = nTotalCandidates;
+        segmentBlockSize[sI] = nCandidates;
+
+        nTotalCandidates += nCandidates;
+    }
+
+    // Needed for reverseDistribute
+    label preDistributionToCandidateProcSize = toCandidateProc.size();
+
+    autoPtr<mapDistribute> map(buildMap(toCandidateProc));
+
+    map().distribute(testStarts);
+    map().distribute(testEnds);
+
+    List<pointIndexHit> segmentIntersectsCandidate(testStarts.size());
+
+    // Test candidate segments on candidate processors
+    forAll(testStarts, sI)
+    {
+        const point& s = testStarts[sI];
+        const point& e = testEnds[sI];
+
+        // If the sphere finds some elements of the patch, then it overlaps
+        segmentIntersectsCandidate[sI] = bFTreePtr_().findLine(s, e);
+    }
+
+    map().reverseDistribute
+    (
+        preDistributionToCandidateProcSize,
+        segmentIntersectsCandidate
+    );
+
+    List<List<pointIndexHit> > segmentHitProcs(starts.size());
+
+    // Working storage for assessing processors
+    DynamicList<pointIndexHit> tmpProcHits;
+
+    forAll(starts, sI)
+    {
+        tmpProcHits.clear();
+
+        // Extract the sub list of results for this point
+
+        SubList<pointIndexHit> segmentProcResults
+        (
+            segmentIntersectsCandidate,
+            segmentBlockSize[sI],
+            segmentBlockStart[sI]
+        );
+
+        forAll(segmentProcResults, sPRI)
+        {
+            if (segmentProcResults[sPRI].hit())
+            {
+                tmpProcHits.append(segmentProcResults[sPRI]);
+
+                tmpProcHits.last().setIndex
+                (
+                    toCandidateProc[segmentBlockStart[sI] + sPRI]
+                );
+            }
+        }
+
+        segmentHitProcs[sI] = tmpProcHits;
+    }
+
+    return segmentHitProcs;
+}
+
+
+Foam::labelListList Foam::backgroundMeshDecomposition::overlapsProcessors
+(
+    const pointField& centres,
+    const scalarField& radiusSqrs,
+    bool includeOwnProcessor
+) const
+{
+    DynamicList<label> toCandidateProc;
+    DynamicList<point> testCentres;
+    DynamicList<scalar> testRadiusSqrs;
+    labelList sphereBlockStart(centres.size(), -1);
+    labelList sphereBlockSize(centres.size(), -1);
+
+    label nTotalCandidates = 0;
+
+    forAll(centres, sI)
+    {
+        const point& c = centres[sI];
+        scalar rSqr = radiusSqrs[sI];
+
+        label nCandidates = 0;
+
+        forAll(allBackgroundMeshBounds_, procI)
+        {
+            // It is assumed that the sphere in question overlaps the source
+            // processor, so don't test it, unless includeOwnProcessor is true
+            if
+            (
+                (includeOwnProcessor || procI != Pstream::myProcNo())
+             && allBackgroundMeshBounds_[procI].overlaps(c, rSqr)
+            )
+            {
+                toCandidateProc.append(procI);
+                testCentres.append(c);
+                testRadiusSqrs.append(rSqr);
+
+                nCandidates++;
+            }
+        }
+
+        sphereBlockStart[sI] = nTotalCandidates;
+        sphereBlockSize[sI] = nCandidates;
+
+        nTotalCandidates += nCandidates;
+    }
+
+    // Needed for reverseDistribute
+    label preDistributionToCandidateProcSize = toCandidateProc.size();
+
+    autoPtr<mapDistribute> map(buildMap(toCandidateProc));
+
+    map().distribute(testCentres);
+    map().distribute(testRadiusSqrs);
+
+    List<bool> sphereOverlapsCandidate(testCentres.size(), false);
+
+    // Test candidate spheres on candidate processors
+    forAll(testCentres, sI)
+    {
+        const point& c = testCentres[sI];
+        scalar rSqr = testRadiusSqrs[sI];
+
+        // If the sphere finds some elements of the patch, then it overlaps
+        sphereOverlapsCandidate[sI] = !bFTreePtr_().findSphere(c, rSqr).empty();
+    }
+
+    map().reverseDistribute
+    (
+        preDistributionToCandidateProcSize,
+        sphereOverlapsCandidate
+    );
+
+    labelListList sphereProcs(centres.size());
+
+    // Working storage for assessing processors
+    DynamicList<label> tmpProcs;
+
+    forAll(centres, sI)
+    {
+        tmpProcs.clear();
+
+        // Extract the sub list of results for this point
+
+        SubList<bool> sphereProcResults
+        (
+            sphereOverlapsCandidate,
+            sphereBlockSize[sI],
+            sphereBlockStart[sI]
+        );
+
+        forAll(sphereProcResults, sPRI)
+        {
+            if (sphereProcResults[sPRI])
+            {
+                tmpProcs.append(toCandidateProc[sphereBlockStart[sI] + sPRI]);
+            }
+        }
+
+        sphereProcs[sI] = tmpProcs;
+    }
+
+    return sphereProcs;
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
index f6a0b448a6..8a970b02f4 100644
--- a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
+++ b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
@@ -68,7 +68,7 @@ SourceFiles
 #include "treeBoundBox.H"
 #include "primitivePatch.H"
 #include "face.H"
-#include "List.H"
+#include "labelList.H"
 #include "pointField.H"
 #include "indexedOctree.H"
 #include "treeDataPrimitivePatch.H"
@@ -111,6 +111,9 @@ class backgroundMeshDecomposition
         //- Search tree for the boundaryFaces_ patch
         autoPtr<indexedOctree<treeDataBPatch> > bFTreePtr_;
 
+        //- The bounds of all background meshes on all processors
+        treeBoundBoxList allBackgroundMeshBounds_;
+
         //- Decomposition dictionary
         IOdictionary decomposeDict_;
 
@@ -167,6 +170,9 @@ class backgroundMeshDecomposition
         //- Build the surface patch and search tree
         void buildPatchAndTree();
 
+        //- Build a mapDistribute for the supplied destination processor data
+        autoPtr<mapDistribute> buildMap(const List<label>& toProc) const;
+
         //- Disallow default bitwise copy construct
         backgroundMeshDecomposition(const backgroundMeshDecomposition&);
 
@@ -205,7 +211,35 @@ public:
         );
 
         //- Is the given position inside the domain of this decomposition
-        bool positionInThisDomain(const point& pt) const;
+        bool positionOnThisProcessor(const point& pt) const;
+
+        //- What processor is the given position on?
+        labelList processorPosition(const pointField& pts) const;
+
+        //- Which processors are intersected by the line segment, returns all
+        //  processors whose boundary patch is intersected by the sphere.  By
+        //  default this does not return the processor that the query is
+        //  launched from, it is assumed that the point is on that processor.
+        //  The index data member of the pointIndexHit is replaced with the
+        //  processor index.
+        List<List<pointIndexHit> > intersectsProcessors
+        (
+            const pointField& starts,
+            const pointField& ends,
+            bool includeOwnProcessor = false
+        ) const;
+
+        //- Which processors overlap the given sphere, returns all processors
+        //  whose boundary patch is touched by the sphere or whom the sphere is
+        //  inside.  By default this does not return the processor that the
+        //  query is launched from, it is assumed that the point is on that
+        //  processor.
+        labelListList overlapsProcessors
+        (
+            const pointField& centres,
+            const scalarField& radiusSqrs,
+            bool includeOwnProcessor = false
+        ) const;
 
         // Access
 
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index a1dcbbec6d..679cf1f188 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -1837,28 +1837,62 @@ bool Foam::conformalVoronoiMesh::positionOnThisProc
 {
     if (Pstream::parRun())
     {
-        return decomposition_().positionInThisDomain(pt);
+        return decomposition_().positionOnThisProcessor(pt);
     }
 
     return true;
 }
 
 
-Foam::label Foam::conformalVoronoiMesh::positionProc
+Foam::labelList Foam::conformalVoronoiMesh::positionProc
 (
-    const Foam::point& pt
+    const Foam::pointField& pts
 ) const
 {
-    Pout<< "Dummy positionProc" << endl;
-
     if (!Pstream::parRun())
     {
-        return -1;
+        return labelList(pts.size(), -1);
     }
-    else
+
+    List<List<pointIndexHit> > inter = decomposition_().intersectsProcessors
+    (
+        pts,
+        pts + vector::one,
+        false
+    );
+
+    return decomposition_().processorPosition(pts);
+}
+
+
+Foam::List<Foam::List<Foam::pointIndexHit> >
+Foam::conformalVoronoiMesh::intersectsProc
+(
+    const pointField& starts,
+    const pointField& ends
+) const
+{
+    if (!Pstream::parRun())
     {
-        return Pstream::myProcNo();
+        return List<List<pointIndexHit> >(starts.size());
     }
+
+    return decomposition_().intersectsProcessors(starts, ends, false);
+}
+
+
+Foam::labelListList Foam::conformalVoronoiMesh::overlapsProc
+(
+    const pointField& centres,
+    const scalarField& radiusSqrs
+) const
+{
+    if (!Pstream::parRun())
+    {
+        return labelListList(centres.size(), labelList(0));
+    }
+
+    return decomposition_().overlapsProcessors(centres, radiusSqrs, false);
 }
 
 
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 4a18ddfc23..bda1b28917 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -794,7 +794,21 @@ public:
         bool positionOnThisProc(const Foam::point& pt) const;
 
         //- Which processor's domain handles this point
-        label positionProc(const Foam::point& pt) const;
+        labelList positionProc(const Foam::pointField& pts) const;
+
+        //- Which other processors does each line segment intersect
+        List<List<pointIndexHit> > intersectsProc
+        (
+            const pointField& starts,
+            const pointField& ends
+        ) const;
+
+        //- Which other processors does each sphere overlap
+        labelListList overlapsProc
+        (
+            const pointField& centres,
+            const scalarField& radiusSqrs
+        ) const;
 
         //- Conversion functions between point (OpenFOAM) and Point (CGAL)
 
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index d0772251d7..8efc2485d7 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -1214,7 +1214,10 @@ void Foam::conformalVoronoiMesh::parallelInterfaceInfluence
         //     const treeBoundBox& procBb =
         //         geometryToConformTo_.processorMeshBounds()[procI];
 
-        //     if (procBb.overlaps(circumcentre, circumradiusSqr))
+        //     // Extend the circumsphere slightly to ensure that any floating
+        //     // point error in its calculation or the definition of the bounds
+        //     // of each processor are masked.
+        //     if (procBb.overlaps(circumcentre, sqr(1.01)*circumradiusSqr))
         //     {
         //         toProc[procI] = true;
         //     }