From 1eb8b8ef124f06519f4370e49a93c57325771a14 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Mon, 13 Jun 2011 15:44:29 +0100
Subject: [PATCH] ENH: ...nearest... query, remove duplicate map build.

---
 .../backgroundMeshDecomposition.C             | 176 +++++++++++++++++-
 .../backgroundMeshDecomposition.H             |   7 +
 .../conformalVoronoiMesh.C                    |   7 -
 .../conformalVoronoiMesh.H                    |   8 +-
 .../conformalVoronoiMeshConformToSurface.C    |  87 +--------
 5 files changed, 184 insertions(+), 101 deletions(-)

diff --git a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index a5667885ad..e89096b70e 100644
--- a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -711,6 +711,17 @@ void Foam::backgroundMeshDecomposition::buildPatchAndTree()
     Pstream::gatherList(allBackgroundMeshBounds_);
     Pstream::scatterList(allBackgroundMeshBounds_);
 
+    point bbMin(GREAT, GREAT, GREAT);
+    point bbMax(-GREAT, -GREAT, -GREAT);
+
+    forAll(allBackgroundMeshBounds_, procI)
+    {
+        bbMin = min(bbMin, allBackgroundMeshBounds_[procI].min());
+        bbMax = max(bbMax, allBackgroundMeshBounds_[procI].max());
+    }
+
+    globalBackgroundBounds_ = treeBoundBox(bbMin, bbMax);
+
     if (debug)
     {
         OFstream fStr
@@ -783,6 +794,7 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
     boundaryFacesPtr_(),
     bFTreePtr_(),
     allBackgroundMeshBounds_(Pstream::nProcs()),
+    globalBackgroundBounds_(),
     decomposeDict_
     (
         IOobject
@@ -1167,6 +1179,9 @@ Foam::labelList Foam::backgroundMeshDecomposition::processorPosition
 
     labelList ptProc(pts.size(), -1);
 
+    DynamicList<label> failedPointIndices;
+    DynamicList<point> failedPoints;
+
     forAll(pts, pI)
     {
         // Extract the sub list of results for this point
@@ -1190,25 +1205,178 @@ Foam::labelList Foam::backgroundMeshDecomposition::processorPosition
 
         if (ptProc[pI] < 0)
         {
-            // No processor was found
+            if (!globalBackgroundBounds_.contains(pts[pI]))
+            {
+                FatalErrorIn
+                (
+                    "Foam::labelList"
+                    "Foam::backgroundMeshDecomposition::processorPosition"
+                    "("
+                        "const List<point>& pts"
+                    ") const"
+                )
+                    << "The position " << pts[pI]
+                    << " is not in any part of the background mesh.  "
+                    << "A background mesh with a wider margin around "
+                    << "the geometry may help."
+                    << exit(FatalError);
+            }
+
+            if (debug)
+            {
+                WarningIn
+                (
+                    "Foam::labelList"
+                    "Foam::backgroundMeshDecomposition::processorPosition"
+                    "("
+                        "const List<point>& pts"
+                    ") const"
+                )
+                    << "The position " << pts[pI]
+                    << " was not found in the background mesh, finding nearest."
+                    << endl;
+            }
+
+            failedPointIndices.append(pI);
+            failedPoints.append(pts[pI]);
+        }
+    }
+
+    labelList ptNearestProc(processorNearestPosition(failedPoints));
+
+    forAll(failedPoints, fPI)
+    {
+        label pI = failedPointIndices[fPI];
+
+        ptProc[pI] = ptNearestProc[fPI];
+    }
+
+    return ptProc;
+}
+
+
+Foam::labelList Foam::backgroundMeshDecomposition::processorNearestPosition
+(
+    const List<point>& 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)
+        {
+            // Candidate points may lie just outside a processor box, increase
+            // test range by using overlaps rather than contains
+            if (allBackgroundMeshBounds_[procI].overlaps(pt, sqr(SMALL*100)))
+            {
+                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<scalar> distanceSqrToCandidate(testPoints.size(), sqr(GREAT));
+
+    // Test candidate points on candidate processors
+    forAll(testPoints, tPI)
+    {
+        pointIndexHit info = bFTreePtr_().findNearest
+        (
+            testPoints[tPI],
+            sqr(GREAT)
+        );
+
+        if (info.hit())
+        {
+            distanceSqrToCandidate[tPI] = magSqr
+            (
+                testPoints[tPI] - info.hitPoint()
+            );
+        }
+    }
+
+    map().reverseDistribute
+    (
+        preDistributionToCandidateProcSize,
+        distanceSqrToCandidate
+    );
+
+    labelList ptNearestProc(pts.size(), -1);
+
+    forAll(pts, pI)
+    {
+        // Extract the sub list of results for this point
+
+        SubList<scalar> ptNearestProcResults
+        (
+            distanceSqrToCandidate,
+            ptBlockSize[pI],
+            ptBlockStart[pI]
+        );
+
+        scalar nearestProcDistSqr = GREAT;
+
+        forAll(ptNearestProcResults, pPRI)
+        {
+            if (ptNearestProcResults[pPRI] < nearestProcDistSqr)
+            {
+                nearestProcDistSqr = ptNearestProcResults[pPRI];
+
+                ptNearestProc[pI] = toCandidateProc[ptBlockStart[pI] + pPRI];
+            }
+        }
+
+        if (debug)
+        {
+            Pout<< pts[pI] << " nearestProcDistSqr " << nearestProcDistSqr
+                << " ptNearestProc[pI] " << ptNearestProc[pI] << endl;
+        }
+
+        if (ptNearestProc[pI] < 0)
+        {
             FatalErrorIn
             (
                 "Foam::labelList"
-                "Foam::backgroundMeshDecomposition::processorPosition"
+                "Foam::backgroundMeshDecomposition::processorNearestPosition"
                 "("
                     "const List<point>& pts"
                 ") const"
             )
                 << "The position " << pts[pI]
-                << " was not found in any part of the background mesh."
+                << " did not find a nearest point on the background mesh."
                 << exit(FatalError);
         }
     }
 
-    return ptProc;
+    return ptNearestProc;
 }
 
 
+
 Foam::List<Foam::List<Foam::pointIndexHit> >
 Foam::backgroundMeshDecomposition::intersectsProcessors
 (
diff --git a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
index 757beb2afe..a9462bf515 100644
--- a/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
+++ b/src/mesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
@@ -114,6 +114,10 @@ class backgroundMeshDecomposition
         //- The bounds of all background meshes on all processors
         treeBoundBoxList allBackgroundMeshBounds_;
 
+        //- The overall bounds of all of the background meshes, used to test if
+        //  a point that is not found on any processor is in the domain at all
+        treeBoundBox globalBackgroundBounds_;
+
         //- Decomposition dictionary
         IOdictionary decomposeDict_;
 
@@ -230,6 +234,9 @@ public:
         //- What processor is the given position on?
         labelList processorPosition(const List<point>& pts) const;
 
+        //- What is the nearest processor to the given position?
+        labelList processorNearestPosition(const List<point>& 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
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 092ec2dc0f..73433b7ecb 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -1992,13 +1992,6 @@ Foam::labelList Foam::conformalVoronoiMesh::positionProc
         return labelList(pts.size(), -1);
     }
 
-    List<List<pointIndexHit> > inter = decomposition_().intersectsProcessors
-    (
-        pts,
-        pts + vector::one,
-        false
-    );
-
     return decomposition_().processorPosition(pts);
 }
 
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index f280516d52..90467dd0b8 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -486,12 +486,6 @@ private:
             const word& outputName
         );
 
-        //- Build map to refer vertices to target processors
-        mapDistribute buildReferringMap
-        (
-            const DynamicList<label>& targetProcessor
-        ) const;
-
         //- Refer vertices to their required processors
         void referVertices
         (
@@ -673,7 +667,7 @@ private:
         ) const;
 
         //- Collapse a face to an edge, updating the point and point
-        //- map.  Returns the collapse mode that was applied.
+        // map.  Returns the collapse mode that was applied.
         faceCollapseMode collapseFace
         (
             const face& f,
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index 38ecdf6e0b..283ba501c9 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -1065,88 +1065,6 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
 }
 
 
-Foam::mapDistribute Foam::conformalVoronoiMesh::buildReferringMap
-(
-    const DynamicList<label>& targetProcessor
-) const
-{
-    // Determine send map
-    // ~~~~~~~~~~~~~~~~~~
-
-    // 1. Count
-    labelList nSend(Pstream::nProcs(), 0);
-
-    forAll(targetProcessor, i)
-    {
-        label procI = targetProcessor[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(targetProcessor, i)
-    {
-        label procI = targetProcessor[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 mapDistribute
-    (
-        constructSize,
-        sendMap.xfer(),
-        constructMap.xfer()
-    );
-}
-
-
 void Foam::conformalVoronoiMesh::referVertices
 (
     const DynamicList<label>& targetProcessor,
@@ -1166,7 +1084,10 @@ void Foam::conformalVoronoiMesh::referVertices
         );
     }
 
-    mapDistribute pointMap = buildReferringMap(targetProcessor);
+    mapDistribute pointMap = backgroundMeshDecomposition::buildMap
+    (
+        targetProcessor
+    );
 
     label totalVertices = parallelPoints.size();