From d90fa36e96991894d5dfedf5b0a11b47dd45008f Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Fri, 18 Mar 2011 15:08:15 +0000
Subject: [PATCH] ENH: Fixed parallel interface identification, restrict
 surface intersection test to box.

---
 .../conformalVoronoiMesh.C                    |   4 +-
 .../conformalVoronoiMeshConformToSurface.C    | 142 ++++++++++++++----
 .../uniformGrid/uniformGrid.C                 |   6 +-
 3 files changed, 114 insertions(+), 38 deletions(-)

diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 72886df34b..40137ddd84 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -493,13 +493,13 @@ void Foam::conformalVoronoiMesh::createFeaturePoints()
     pointField farPts = geometryToConformTo_.globalBounds().points();
 
     // Shift corners of bounds relative to origin
-    farPts -= geometryToConformTo_.bounds().midpoint();
+    farPts -= geometryToConformTo_.globalBounds().midpoint();
 
     // Scale the box up
     farPts *= 2.0;
 
     // Shift corners of bounds back to be relative to midpoint
-    farPts += geometryToConformTo_.bounds().midpoint();
+    farPts += geometryToConformTo_.globalBounds().midpoint();
 
     forAll(farPts, fPI)
     {
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index ad7347edd2..048b10b865 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -385,17 +385,17 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
                     //                 DynamicList<label>& iVertexToProc =
                     //                     verticesToProc[ivIndex];
 
-                                    // if
-                                    // (
-                                    //     findIndex
-                                    //     (
-                                    //         iVertexToProc,
-                                    //         toProc
-                                    //     ) == -1
-                                    // )
-                                    // {
-                                    //     iVertexToProc.append(toProc);
-                                    // }
+                    //                 if
+                    //                 (
+                    //                     findIndex
+                    //                     (
+                    //                         iVertexToProc,
+                    //                         toProc
+                    //                     ) == -1
+                    //                 )
+                    //                 {
+                    //                     iVertexToProc.append(toProc);
+                    //                 }
                     //             }
                     //         }
                     //     }
@@ -451,11 +451,14 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
                 }
             }
 
-            writePoints
-            (
-                "parallelInterfacePointsToSend_" + name(nIter) + ".obj",
-                parallelInterfacePoints
-            );
+            if (cvMeshControls().objOutput())
+            {
+                writePoints
+                (
+                    "parallelInterfacePointsToSend_" + name(nIter) + ".obj",
+                    parallelInterfacePoints
+                );
+            }
 
             // Determine send map
             // ~~~~~~~~~~~~~~~~~~
@@ -540,11 +543,14 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
 
             pointMap.distribute(parallelInterfaceIndices);
 
-            writePoints
-            (
-                "parallelInterfacePointsReceived_"+ name(nIter) + ".obj",
-                parallelInterfacePoints
-            );
+            if (cvMeshControls().objOutput())
+            {
+                writePoints
+                (
+                    "parallelInterfacePointsReceived_"+ name(nIter) + ".obj",
+                    parallelInterfacePoints
+                );
+            }
 
             for (label domain = 0; domain < Pstream::nProcs(); domain++)
             {
@@ -903,6 +909,12 @@ Foam::conformalVoronoiMesh::parallelInterfaceIntersection
     std::list<Facet> facets;
     incident_facets(vit, std::back_inserter(facets));
 
+    // Pout<< "# " << label(vit->index()) << " facets size "
+    //     << label(facets.size())
+    //     << endl;
+
+    // bool anyCuts = false;
+
     DynamicList<label> procs;
 
     for
@@ -927,20 +939,72 @@ Foam::conformalVoronoiMesh::parallelInterfaceIntersection
         Foam::point boxPt(vector::one*GREAT);
         direction ptOnFace = -1;
 
-        geometryToConformTo_.bounds().intersects
+        // bool cuts =
+        //     (
+        //         geometryToConformTo_.bounds().contains(dE0)
+        //      && !geometryToConformTo_.bounds().contains(dE1)
+        //     )
+        //  || (
+        //         !geometryToConformTo_.bounds().contains(dE0)
+        //      && geometryToConformTo_.bounds().contains(dE1)
+        //     );
+
+        // anyCuts =
+        //      !geometryToConformTo_.bounds().contains(dE1)
+        //   ||
+        //      !geometryToConformTo_.bounds().contains(dE0);
+
+        Foam::point start = dE0;
+        Foam::point end   = dE1;
+
+        bool intersects = geometryToConformTo_.bounds().intersects
         (
-            dE0,
-            dE1 - dE0,
-            dE0,
-            dE1,
+            start,
+            end - start,
+            start,
+            end,
             boxPt,
             ptOnFace
         );
 
+        if (intersects && ptOnFace == 0)
+        {
+            // If the box is intersected, but doesn't return the appropriate
+            // bits, then this means that the start point was inside the box,
+            // so reverse the direction of the query.
+
+            start = dE1;
+            end   = dE0;
+
+            geometryToConformTo_.bounds().intersects
+            (
+                start,
+                end - start,
+                start,
+                end,
+                boxPt,
+                ptOnFace
+            );
+        }
+
+        // Pout<< "# " << label(vit->index())
+        //     << " Dual edge cuts " << cuts
+        //     << " intersects " << intersects
+        //     << " ptOnFace " << ptOnFace
+        //     << endl;
+
+        // meshTools::writeOBJ(Pout, dE0);
+        // meshTools::writeOBJ(Pout, dE1);
+
+        // Pout << "l dE0 dE1" << endl;
+
+        // If an intersection with the box has been found, then it is the "end"
+        // point that is inside the box, so test the line from end to
+        // boxPt to see if the surface is in-between.
         if
         (
             ptOnFace > 0
-         && !geometryToConformTo_.findSurfaceAnyIntersection(dE0, dE1)
+         && !geometryToConformTo_.findSurfaceAnyIntersection(end, boxPt)
         )
         {
             // A surface penetration is not found, and a bounds penetration is
@@ -960,12 +1024,19 @@ Foam::conformalVoronoiMesh::parallelInterfaceIntersection
             // Pout<< "Parallel box penetration, face " << ptOnFace
             //     << " proc " << target << endl;
             // meshTools::writeOBJ(Pout, dE0);
-            // meshTools::writeOBJ(Pout, dE1);
             // meshTools::writeOBJ(Pout, boxPt);
-            // Pout << "l dE0 dE1" << endl;
+            // meshTools::writeOBJ(Pout, dE1);
+            // Pout << "l dE0 boxPt dE1" << endl;
         }
     }
 
+    // if (anyCuts && procs.empty())
+    // {
+    //     Pout<< "# " << label(vit->index())
+    //         << " anyCuts && procs.empty()"
+    //         << endl;
+    // }
+
     return procs;
 }
 
@@ -977,6 +1048,11 @@ Foam::label Foam::conformalVoronoiMesh::targetProc
 {
     // Hard coded to 8 proc, bound-box octant split:
 
+    if (ptOnFace == 0)
+    {
+        return -1;
+    }
+
     label faceIndex = 0;
 
     // From: http://graphics.stanford.edu/~seander/bithacks.html
@@ -987,10 +1063,10 @@ Foam::label Foam::conformalVoronoiMesh::targetProc
 
     label procArray[8][6] =
     {
-        {-1,  1, -1,  2, -1, 4},
-        { 0, -1, -1,  3, -1, 5},
-        {-1,  3,  0, -1, -1, 6},
-        { 2, -1,  1, -1, -1, 7},
+        {-1,  1, -1,  2, -1,  4},
+        { 0, -1, -1,  3, -1,  5},
+        {-1,  3,  0, -1, -1,  6},
+        { 2, -1,  1, -1, -1,  7},
         {-1,  5, -1,  6,  0, -1},
         { 4, -1, -1,  7,  1, -1},
         {-1,  7,  4, -1,  2, -1},
diff --git a/src/mesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C b/src/mesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C
index 552e67c2ca..d5895dc955 100644
--- a/src/mesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C
+++ b/src/mesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C
@@ -100,9 +100,9 @@ std::vector<Vb::Point> uniformGrid::initialPoints() const
             {
                 point p
                 (
-                    x0 + i*delta.x(),
-                    y0 + j*delta.y(),
-                    z0 + k*delta.z()
+                    x0 + (i + 0.5)*delta.x(),
+                    y0 + (j + 0.5)*delta.y(),
+                    z0 + (k + 0.5)*delta.z()
                 );
 
                 if (randomiseInitialGrid_)