From 09a44ad6064fc910132d1485426a20a50cc0371c Mon Sep 17 00:00:00 2001
From: laurence <laurence>
Date: Tue, 30 Jul 2013 12:12:18 +0100
Subject: [PATCH] ENH: foamyHexMesh: Check moved internal points are inside the
 geometry before inserting them

---
 .../conformalVoronoiMesh.C                    | 114 ++++++------------
 .../conformalVoronoiMesh.H                    |   3 +
 .../conformalVoronoiMeshI.H                   |  18 +++
 .../conformalVoronoiMeshIO.C                  |   4 +
 4 files changed, 60 insertions(+), 79 deletions(-)

diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 14d1042778..d87d1ed97b 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -562,7 +562,7 @@ void Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()
     {
         if (vit->internalOrBoundaryPoint())
         {
-            pointFromPoint pt = topoint(vit->point());
+            const pointFromPoint pt = topoint(vit->point());
 
             cellShapeControls().cellSizeAndAlignment
             (
@@ -1135,10 +1135,12 @@ void Foam::conformalVoronoiMesh::move()
                      && pointToBeRetained[vB->index()] == true
                     )
                     {
-                        pointsToInsert.append
-                        (
-                            toPoint(0.5*(dVA + dVB))
-                        );
+                        Foam::point pt(0.5*(dVA + dVB));
+
+                        if (internalPointIsInside(pt))
+                        {
+                            pointsToInsert.append(toPoint(pt));
+                        }
                     }
                 }
 
@@ -1274,7 +1276,10 @@ void Foam::conformalVoronoiMesh::move()
                             if (positionOnThisProc(newPt))
                             {
                                 // Prevent insertions spanning surfaces
-                                pointsToInsert.append(toPoint(newPt));
+                                if (internalPointIsInside(newPt))
+                                {
+                                    pointsToInsert.append(toPoint(newPt));
+                                }
                             }
                         }
                     }
@@ -1310,10 +1315,12 @@ void Foam::conformalVoronoiMesh::move()
                              && pointToBeRetained[vB->index()] == true
                             )
                             {
-                                pointsToInsert.append
-                                (
-                                    toPoint(0.5*(dVA + dVB))
-                                );
+                                Foam::point pt(0.5*(dVA + dVB));
+
+                                if (internalPointIsInside(pt))
+                                {
+                                    pointsToInsert.append(toPoint(pt));
+                                }
                             }
                         }
 
@@ -1428,21 +1435,31 @@ void Foam::conformalVoronoiMesh::move()
                 // 14/1/2011.
                 // Only necessary if using an exact constructions kernel
                 // (extended precision)
-
-                pointsToInsert.append
+                Foam::point pt
                 (
-                    toPoint
-                    (
-                        topoint(vit->point())
-                      + displacementAccumulator[vit->index()]
-                    )
+                    topoint(vit->point())
+                  + displacementAccumulator[vit->index()]
                 );
+
+                if (internalPointIsInside(pt))
+                {
+                    pointsToInsert.append(toPoint(pt));
+                }
             }
         }
     }
 
     pointsToInsert.shrink();
 
+    Info<< indent
+        << returnReduce
+           (
+               pointToBeRetained.count() - pointsToInsert.size(),
+               sumOp<label>()
+           )
+        << " internal points were inserted outside the domain. "
+        << "They have been removed." << endl;
+
     // Save displacements to file.
     if (foamyHexMeshControls().objOutput() && time().outputTime())
     {
@@ -1477,73 +1494,12 @@ void Foam::conformalVoronoiMesh::move()
 
     timeCheck("Displacement calculated");
 
-    Info<< nl << "Inserting displaced tessellation" << endl;
+    Info<< nl<< "Inserting displaced tessellation" << endl;
 
     insertFeaturePoints(true);
 
     insertInternalPoints(pointsToInsert, true);
 
-    {
-        // Remove internal points that have been inserted outside the surface.
-        label internalPtIsOutside = 0;
-
-        autoPtr<OBJstream> str;
-
-        if (foamyHexMeshControls().objOutput() && time().outputTime())
-        {
-            str.set
-            (
-                new OBJstream(time().path()/"internalPointsOutsideDomain.obj")
-            );
-        }
-
-        DynamicList<Vertex_handle> pointsToRemove;
-
-        for
-        (
-            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-            vit != finite_vertices_end();
-            ++vit
-        )
-        {
-            if
-            (
-                (vit->internalPoint() || vit->internalBoundaryPoint())
-             //&& !vit->referred()
-            )
-            {
-                const pointFromPoint pt = topoint(vit->point());
-
-                bool inside = geometryToConformTo_.inside(pt);
-
-                if
-                (
-                    !inside
-                 || !geometryToConformTo_.globalBounds().contains(pt)
-                )
-                {
-                    if
-                    (
-                        foamyHexMeshControls().objOutput()
-                     && time().outputTime()
-                    )
-                    {
-                        str().write(topoint(vit->point()));
-                    }
-
-                    pointsToRemove.append(vit);
-                    internalPtIsOutside++;
-                }
-            }
-        }
-
-        remove(pointsToRemove.begin(), pointsToRemove.end());
-
-        Info<< "    " << returnReduce(internalPtIsOutside, sumOp<label>())
-            << " internal points were inserted outside the domain. "
-            << "They have been removed." << endl;
-    }
-
     // Fix points that have not been significantly displaced
 //    for
 //    (
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 6190b0418e..7e78d35557 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -244,6 +244,9 @@ private:
             DynamicList<Vb>& pts
         ) const;
 
+        //- Check internal point is completely inside the meshable region
+        inline bool internalPointIsInside(const Foam::point& pt) const;
+
         inline bool isPointPair
         (
             const Vertex_handle& vA,
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
index def7120d55..64534cb839 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
@@ -325,6 +325,24 @@ inline void Foam::conformalVoronoiMesh::createBafflePointPair
 }
 
 
+inline bool Foam::conformalVoronoiMesh::internalPointIsInside
+(
+    const Foam::point& pt
+) const
+{
+    if
+    (
+        !geometryToConformTo_.inside(pt)
+     || !geometryToConformTo_.globalBounds().contains(pt)
+    )
+    {
+        return false;
+    }
+
+    return true;
+}
+
+
 inline bool Foam::conformalVoronoiMesh::isPointPair
 (
     const Vertex_handle& vA,
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 9c33e9f5ab..bb34cbd4a6 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -674,6 +674,8 @@ void Foam::conformalVoronoiMesh::reorderProcessorPatches
 
     pBufs.finishedSends();
 
+    Info<< incrIndent << indent << "    Face ordering initialised..." << endl;
+
     // Receive and calculate ordering
     bool anyChanged = false;
 
@@ -731,6 +733,8 @@ void Foam::conformalVoronoiMesh::reorderProcessorPatches
         }
     }
 
+    Info<< incrIndent << indent << "    Faces matched." << endl;
+
     reduce(anyChanged, orOp<bool>());
 
     if (anyChanged)