From 838639cd5dc138c3d98fe3db76fc2d3dc067c476 Mon Sep 17 00:00:00 2001
From: laurence <laurence>
Date: Wed, 12 Jun 2013 13:08:15 +0100
Subject: [PATCH] ENH: Snapping to surfaces if point is far outside

---
 .../conformalVoronoiMeshCalcDualMesh.C        | 87 +++++++++++++++++++
 1 file changed, 87 insertions(+)

diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index 9cfa9e316c..c771256f4b 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -1787,6 +1787,11 @@ void Foam::conformalVoronoiMesh::indexDualVertices
         }
     }
 
+    OBJstream snapping1("snapToSurface1.obj");
+    OBJstream snapping2("snapToSurface2.obj");
+    OFstream tetToSnapTo("tetsToSnapTo.obj");
+    label offset = 0;
+
     for
     (
         Delaunay::Finite_cells_iterator cit = finite_cells_begin();
@@ -1892,6 +1897,88 @@ void Foam::conformalVoronoiMesh::indexDualVertices
                 }
             }
 
+            {
+                // Snapping points far outside
+                if (cit->boundaryDualVertex())
+                {
+                    pointFromPoint dual = cit->dual();
+
+                    pointIndexHit hitInfo;
+                    label surfHit;
+
+                    // Find nearest surface point
+                    geometryToConformTo_.findSurfaceNearest
+                    (
+                        dual,
+                        sqr(targetCellSize(dual)),
+                        hitInfo,
+                        surfHit
+                    );
+
+                    if (!hitInfo.hit())
+                    {
+                        // Project dual to nearest point on tet
+
+                        tetPointRef tet
+                        (
+                            topoint(cit->vertex(0)->point()),
+                            topoint(cit->vertex(1)->point()),
+                            topoint(cit->vertex(2)->point()),
+                            topoint(cit->vertex(3)->point())
+                        );
+
+                        pointFromPoint nearestPointOnTet =
+                            tet.nearestPoint(dual).rawPoint();
+
+                        // Get nearest point on surface from tet.
+                        geometryToConformTo_.findSurfaceNearest
+                        (
+                            nearestPointOnTet,
+                            sqr(targetCellSize(nearestPointOnTet)),
+                            hitInfo,
+                            surfHit
+                        );
+
+                        vector snapDir = nearestPointOnTet - dual;
+                        snapDir /= mag(snapDir) + SMALL;
+
+                        drawDelaunayCell(tetToSnapTo, cit, offset);
+                        offset += 1;
+
+                        vectorField norm(1);
+                        allGeometry_[surfHit].getNormal
+                        (
+                            List<pointIndexHit>(1, hitInfo),
+                            norm
+                        );
+                        norm[0] /= mag(norm[0]) + SMALL;
+
+                        if
+                        (
+                            hitInfo.hit()
+                         && (mag(snapDir & norm[0]) > 0.5)
+                        )
+                        {
+                            snapping1.write
+                            (
+                                linePointRef(dual, nearestPointOnTet)
+                            );
+
+                            snapping2.write
+                            (
+                                linePointRef
+                                (
+                                    nearestPointOnTet,
+                                    hitInfo.hitPoint()
+                                )
+                            );
+
+                            pts[cit->cellIndex()] = hitInfo.hitPoint();
+                        }
+                    }
+                }
+            }
+
             if (cit->boundaryDualVertex())
             {
                 if (cit->featureEdgeDualVertex())