From f37e7ff7d95ba8570f3055f5966eea34d88c2c05 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Mon, 4 Apr 2011 17:19:45 +0100
Subject: [PATCH] ENH: 2 or 8 procs.  pointFile only inserts on one proc.

---
 .../conformalVoronoiMeshConformToSurface.C    | 287 ++++++++++++------
 .../conformalVoronoiMeshIO.C                  |  12 +-
 .../conformationSurfaces.C                    |  70 ++++-
 .../conformationSurfaces.H                    |   3 +
 .../conformationSurfacesI.H                   |   7 +
 .../initialPointsMethod/pointFile/pointFile.C |  58 ++++
 6 files changed, 324 insertions(+), 113 deletions(-)

diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index 0f8af79be3..0f3e2d3846 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -563,123 +563,206 @@ void Foam::conformalVoronoiMesh::buildParallelInterface
         return;
     }
 
+    {
+        // Update the bounds of the domain
+
+        Info<< "USING SINGLE PROCESSOR PER DOMAIN" << endl;
+
+        DynamicList<Foam::point> parallelAllPoints;
+        DynamicList<label> targetProcessor;
+        DynamicList<label> parallelAllIndices;
+
+        Foam::point minPt(VGREAT, VGREAT, VGREAT);
+        Foam::point maxPt(-VGREAT, -VGREAT, -VGREAT);
+
+        for
+        (
+            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+            vit != finite_vertices_end();
+            vit++
+        )
+        {
+            if (vit->real())
+            {
+                Foam::point v = topoint(vit->point());
+
+                minPt = Foam::min(minPt, v);
+                maxPt = Foam::max(maxPt, v);
+            }
+        }
+
+        Pout<< "Before processorDomains update"
+            << geometryToConformTo_.processorDomains()[Pstream::myProcNo()]
+            << endl;
+
+        geometryToConformTo_.processorDomains()[Pstream::myProcNo()][0] =
+            treeBoundBox(minPt, maxPt);
+
+        Pout<< "After processorDomains update"
+            << geometryToConformTo_.processorDomains()[Pstream::myProcNo()]
+            << endl;
+
+        {
+            OFstream str
+            (
+                runTime_.path()
+                /"procDomainUpdated_"
+              + name(Pstream::myProcNo())
+              + "_bounds.obj"
+            );
+
+            Pout<< "Writing " << str.name() << endl;
+
+            pointField bbPoints
+            (
+                geometryToConformTo_.processorDomains()
+                [
+                    Pstream::myProcNo()
+                ][0].points()
+            );
+
+            forAll(bbPoints, i)
+            {
+                meshTools::writeOBJ(str, bbPoints[i]);
+            }
+
+            forAll(treeBoundBox::faces, i)
+            {
+                const face& f = treeBoundBox::faces[i];
+
+                str << "f"
+                    << ' ' << f[0] + 1
+                    << ' ' << f[1] + 1
+                    << ' ' << f[2] + 1
+                    << ' ' << f[3] + 1
+                    << nl;
+            }
+        }
+
+        Pstream::gatherList(geometryToConformTo_.processorDomains());
+
+        Pstream::scatterList(geometryToConformTo_.processorDomains());
+    }
+
     boolList sendToProc(Pstream::nProcs(), false);
 
-    // {
-    //     // Refer all points to all processors
+    bool allPointReferral = false;
 
-    //     DynamicList<Foam::point> parallelAllPoints;
-    //     DynamicList<label> targetProcessor;
-    //     DynamicList<label> parallelAllIndices;
+    if (allPointReferral)
+    {
+        // Refer all points to all processors
 
-    //     for
-    //     (
-    //         Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-    //         vit != finite_vertices_end();
-    //         vit++
-    //     )
-    //     {
-    //         if (!vit->farPoint())
-    //         {
-    //             sendToProc = true;
+        DynamicList<Foam::point> parallelAllPoints;
+        DynamicList<label> targetProcessor;
+        DynamicList<label> parallelAllIndices;
 
-    //             sendToProc[Pstream::myProcNo()] = false;
+        for
+        (
+            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+            vit != finite_vertices_end();
+            vit++
+        )
+        {
+            if (!vit->farPoint())
+            {
+                sendToProc = true;
 
-    //             forAll(sendToProc, procI)
-    //             {
-    //                 if (sendToProc[procI])
-    //                 {
-    //                     label vIndex = vit->index();
+                sendToProc[Pstream::myProcNo()] = false;
 
-    //                     // Using the hashSet to ensure that each vertex is
-    //                     // only referred once to each processor
-    //                     if (!referralVertices[procI].found(vIndex))
-    //                     {
-    //                         referralVertices[procI].insert(vIndex);
+                forAll(sendToProc, procI)
+                {
+                    if (sendToProc[procI])
+                    {
+                        label vIndex = vit->index();
 
-    //                         parallelAllPoints.append
-    //                         (
-    //                             topoint(vit->point())
-    //                         );
+                        // Using the hashSet to ensure that each vertex is
+                        // only referred once to each processor
+                        if (!referralVertices[procI].found(vIndex))
+                        {
+                            referralVertices[procI].insert(vIndex);
 
-    //                         targetProcessor.append(procI);
+                            parallelAllPoints.append
+                            (
+                                topoint(vit->point())
+                            );
 
-    //                         if (vit->internalOrBoundaryPoint())
-    //                         {
-    //                             parallelAllIndices.append(vIndex);
-    //                         }
-    //                         else
-    //                         {
-    //                             parallelAllIndices.append(-vIndex);
-    //                         }
-    //                     }
-    //                 }
-    //             }
-    //         }
-    //     }
+                            targetProcessor.append(procI);
 
-    //     if (cvMeshControls().objOutput())
-    //     {
-    //         writePoints
-    //         (
-    //             "parallelAllPointsToSend.obj",
-    //             parallelAllPoints
-    //         );
-    //     }
+                            if (vit->internalOrBoundaryPoint())
+                            {
+                                parallelAllIndices.append(vIndex);
+                            }
+                            else
+                            {
+                                parallelAllIndices.append(-vIndex);
+                            }
+                        }
+                    }
+                }
+            }
+        }
 
-    //     mapDistribute pointMap = buildReferringMap(targetProcessor);
+        if (cvMeshControls().objOutput())
+        {
+            writePoints
+            (
+                "parallelAllPointsToSend.obj",
+                parallelAllPoints
+            );
+        }
 
-    //     label totalAllVertices = parallelAllPoints.size();
+        mapDistribute pointMap = buildReferringMap(targetProcessor);
 
-    //     reduce(totalAllVertices, sumOp<label>());
+        label totalAllVertices = parallelAllPoints.size();
 
-    //     pointMap.distribute(parallelAllPoints);
+        reduce(totalAllVertices, sumOp<label>());
 
-    //     pointMap.distribute(parallelAllIndices);
+        pointMap.distribute(parallelAllPoints);
 
-    //     if (cvMeshControls().objOutput())
-    //     {
-    //         writePoints
-    //         (
-    //             "parallelAllPointsReceived_" + outputName + ".obj",
-    //             parallelAllPoints
-    //         );
-    //     }
+        pointMap.distribute(parallelAllIndices);
 
-    //     for (label procI = 0; procI < Pstream::nProcs(); procI++)
-    //     {
-    //         const labelList& constructMap = pointMap.constructMap()[procI];
+        if (cvMeshControls().objOutput())
+        {
+            writePoints
+            (
+                "parallelAllPointsReceived_" + outputName + ".obj",
+                parallelAllPoints
+            );
+        }
 
-    //         if (constructMap.size())
-    //         {
-    //             forAll(constructMap, i)
-    //             {
-    //                 label origIndex =
-    //                     parallelAllIndices[constructMap[i]];
+        for (label procI = 0; procI < Pstream::nProcs(); procI++)
+        {
+            const labelList& constructMap = pointMap.constructMap()[procI];
 
-    //                 if (!receivedVertices[procI].found(origIndex))
-    //                 {
-    //                     // For the initial referred vertices, the original
-    //                     // processor is the one that is sending it.
-    //                     label encodedProcI = -(procI + 1);
+            if (constructMap.size())
+            {
+                forAll(constructMap, i)
+                {
+                    label origIndex =
+                        parallelAllIndices[constructMap[i]];
 
-    //                     insertPoint
-    //                     (
-    //                         parallelAllPoints[constructMap[i]],
-    //                         origIndex,
-    //                         encodedProcI
-    //                     );
+                    if (!receivedVertices[procI].found(origIndex))
+                    {
+                        // For the initial referred vertices, the original
+                        // processor is the one that is sending it.
+                        label encodedProcI = -(procI + 1);
 
-    //                     receivedVertices[procI].insert(origIndex);
-    //                 }
-    //             }
-    //         }
-    //     }
+                        insertPoint
+                        (
+                            parallelAllPoints[constructMap[i]],
+                            origIndex,
+                            encodedProcI
+                        );
 
-    //     Info<< "totalAllVertices "
-    //         << totalAllVertices << endl;
-    // }
+                        receivedVertices[procI].insert(origIndex);
+                    }
+                }
+            }
+        }
 
+        Info<< "totalAllVertices "
+            << totalAllVertices << endl;
+    }
 
     if (initialEdgeReferrral)
     {
@@ -779,6 +862,12 @@ void Foam::conformalVoronoiMesh::buildParallelInterface
                         // processor is the one that is sending it.
                         label encodedProcI = -(procI + 1);
 
+                        // Pout<< "Insert "
+                        //     << parallelIntersectionPoints[constructMap[i]]
+                        //     << " " << origIndex
+                        //     << " " << procI
+                        //     << endl;
+
                         insertPoint
                         (
                             parallelIntersectionPoints[constructMap[i]],
@@ -893,6 +982,12 @@ void Foam::conformalVoronoiMesh::buildParallelInterface
                         // processor is the one that is sending it.
                         label encodedProcI = -(procI + 1);
 
+                        // Pout<< "Insert "
+                        //     << parallelInfluencePoints[constructMap[i]]
+                        //     << " " << origIndex
+                        //     << " " << procI
+                        //     << endl;
+
                         insertPoint
                         (
                             parallelInfluencePoints[constructMap[i]],
@@ -1112,14 +1207,8 @@ void Foam::conformalVoronoiMesh::parallelInterfaceInfluence
         );
 
         // Pout<< nl << "# circumradius " << sqrt(circumradiusSqr) << endl;
-
         // drawDelaunayCell(Pout, cit);
 
-        // meshTools::writeOBJ(Pout, topoint(cit->vertex(0)->point()));
-        // meshTools::writeOBJ(Pout, circumcentre);
-
-        // Pout<< "l cr0 cr1" << endl;
-
         forAll(geometryToConformTo_.processorDomains(), procI)
         {
             if (procI == Pstream::myProcNo())
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 10ed5b05a4..bc2d0f0043 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -85,10 +85,14 @@ void Foam::conformalVoronoiMesh::drawDelaunayCell
     const Delaunay::Finite_cells_iterator& cit
 ) const
 {
-    meshTools::writeOBJ(os, topoint(cit->vertex(0)->point()));
-    meshTools::writeOBJ(os, topoint(cit->vertex(1)->point()));
-    meshTools::writeOBJ(os, topoint(cit->vertex(2)->point()));
-    meshTools::writeOBJ(os, topoint(cit->vertex(3)->point()));
+    for (int i = 0; i < 4; i++)
+    {
+        os  << "# index type: "
+            << label(cit->vertex(i)->index()) << " "
+            << label(cit->vertex(i)->type()) << endl;
+
+        meshTools::writeOBJ(os, topoint(cit->vertex(i)->point()));
+    }
 
     os  << "f 1 3 2" << nl
         << "f 2 3 4" << nl
diff --git a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index 149e7c5189..4f0b6c24bc 100644
--- a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -246,9 +246,45 @@ Foam::conformationSurfaces::conformationSurfaces
 
     if (Pstream::parRun())
     {
-        label procLimit = 8;
+        processorDomains_.setSize(Pstream::nProcs());
 
-        if (Pstream::nProcs() != procLimit)
+        if (Pstream::nProcs() == 2)
+        {
+            processorDomains_[Pstream::myProcNo()] = treeBoundBoxList(4);
+
+            forAll(processorDomains_[Pstream::myProcNo()], pDI)
+            {
+                processorDomains_[Pstream::myProcNo()][pDI] =
+                    globalBounds_.subBbox
+                    (
+                        direction(Pstream::myProcNo()*4 + pDI)
+                    );
+            }
+
+            DynamicList<Foam::point> allBbPoints;
+
+            forAll(processorDomains_[Pstream::myProcNo()], pDI)
+            {
+                allBbPoints.append
+                (
+                    processorDomains_[Pstream::myProcNo()][pDI].points()
+                );
+            }
+
+            bounds_ = treeBoundBox(allBbPoints);
+
+            // Replace 4 bound boxes with one
+            processorDomains_[Pstream::myProcNo()] =
+                treeBoundBoxList(1, bounds_);
+        }
+        else if (Pstream::nProcs() == 8)
+        {
+            bounds_ = globalBounds_.subBbox(direction(Pstream::myProcNo()));
+
+            processorDomains_[Pstream::myProcNo()] =
+                treeBoundBoxList(1, bounds_);
+        }
+        else
         {
             FatalErrorIn
             (
@@ -259,16 +295,10 @@ Foam::conformationSurfaces::conformationSurfaces
                     "const dictionary& surfaceConformationDict"
                 ")"
             )
-                << "Hard coded to " << procLimit << " << procs"
+                << "Hard coded to " << "2 or 8" << " procs"
                 << exit(FatalError);
         }
 
-        bounds_ = globalBounds_.subBbox(direction(Pstream::myProcNo()));
-
-        processorDomains_.setSize(Pstream::nProcs());
-
-        processorDomains_[Pstream::myProcNo()] = treeBoundBoxList(1, bounds_);
-
         Pstream::gatherList(processorDomains_);
 
         Pstream::scatterList(processorDomains_);
@@ -379,7 +409,27 @@ bool Foam::conformationSurfaces::positionOnThisProc(const point& pt) const
     // This is likely to give problems when a point is on the boundary between
     // two processors.
 
-    return bounds_.contains(pt);
+    if (Pstream::parRun())
+    {
+        const treeBoundBoxList& procBbs =
+            processorDomains_[Pstream::myProcNo()];
+
+        forAll(procBbs, pBI)
+        {
+            const treeBoundBox& procBb = procBbs[pBI];
+
+            if (procBb.contains(pt))
+            {
+                return true;
+            }
+        }
+    }
+    else
+    {
+        return bounds_.contains(pt);
+    }
+
+    return false;
 }
 
 
diff --git a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
index 4c553407ce..3664186533 100644
--- a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
+++ b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
@@ -166,6 +166,9 @@ public:
             //- Return spanSqr
             inline scalar spanMagSqr() const;
 
+            //- Return non-const access to the processorDomains
+            inline List<treeBoundBoxList>& processorDomains();
+
             //- Return the processorDomains
             inline const List<treeBoundBoxList>& processorDomains() const;
 
diff --git a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfacesI.H b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfacesI.H
index b964e8b85c..b35e45c5d1 100644
--- a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfacesI.H
+++ b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfacesI.H
@@ -74,6 +74,13 @@ inline Foam::scalar Foam::conformationSurfaces::spanMagSqr() const
 }
 
 
+inline Foam::List<Foam::treeBoundBoxList>&
+Foam::conformationSurfaces::processorDomains()
+{
+    return processorDomains_;
+}
+
+
 inline const Foam::List<Foam::treeBoundBoxList>&
 Foam::conformationSurfaces::processorDomains() const
 {
diff --git a/src/mesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C b/src/mesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
index 8755348645..75fb9e6b1a 100644
--- a/src/mesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
+++ b/src/mesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
@@ -85,6 +85,64 @@ std::vector<Vb::Point> pointFile::initialPoints() const
         );
     }
 
+    if (Pstream::parRun())
+    {
+        List<boolList> allProcPt(Pstream::nProcs());
+
+        allProcPt[Pstream::myProcNo()] = procPt;
+
+        Pstream::gatherList(allProcPt);
+
+        Pstream::scatterList(allProcPt);
+
+        // forAll(procPt, ptI)
+        // {
+        //     label nProcsTrue = 0;
+
+        //     forAll(allProcPt, procI)
+        //     {
+        //         if (allProcPt[procI][ptI])
+        //         {
+        //             nProcsTrue++;
+        //         }
+        //     }
+
+        //     if (nProcsTrue != 1)
+        //     {
+        //         Info<< points[ptI] << " " << nProcsTrue << " true" << endl;
+
+        //         if (nProcsTrue > 1)
+        //         {
+        //             forAll(allProcPt, procI)
+        //             {
+        //                 allProcPt[procI][ptI] = false;
+        //             }
+        //         }
+        //     }
+        // }
+
+        forAll(procPt, ptI)
+        {
+            bool foundAlready = false;
+
+            forAll(allProcPt, procI)
+            {
+                // If a processor with a lower index has found this point to
+                // insert already, defer to it and don't insert.
+                if (foundAlready)
+                {
+                    allProcPt[procI][ptI] = false;
+                }
+                else if (allProcPt[procI][ptI])
+                {
+                    foundAlready = true;
+                }
+            }
+        }
+
+        procPt = allProcPt[Pstream::myProcNo()];
+    }
+
     inplaceSubset(procPt, points);
 
     std::vector<Vb::Point> initialPoints;