From 1b19aceb584597bda327453740c34f62368625a9 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 3 Jun 2010 17:44:19 +0100
Subject: [PATCH] BUG: searchableCylinder nearest not correct for endcaps

---
 .../searchableSurface/searchableCylinder.C         | 55 +++++++++++++++++++---
 1 file changed, 49 insertions(+), 6 deletions(-)

diff --git a/src/meshTools/searchableSurface/searchableCylinder.C b/src/meshTools/searchableSurface/searchableCylinder.C
index 0df222cc..70d81e90 100644
--- a/src/meshTools/searchableSurface/searchableCylinder.C
+++ b/src/meshTools/searchableSurface/searchableCylinder.C
@@ -65,6 +65,7 @@ Foam::pointIndexHit Foam::searchableCylinder::findNearest
     v -= parallel*unitDir_;
     scalar magV = mag(v);
 
+
     if (parallel <= 0)
     {
         // nearest is at point1 end cap. Clip to radius.
@@ -74,13 +75,12 @@ Foam::pointIndexHit Foam::searchableCylinder::findNearest
         }
         else
         {
-            //info.setPoint(point1_ + min(magV, radius_)*v/magV);
-            info.setPoint(point1_ + radius_*(v/magV));
+            info.setPoint(point1_ + min(magV, radius_)*v/magV);
         }
     }
     else if (parallel >= magDir_)
     {
-        // nearest is at point2
+        // nearest is at point2 end cap
         if (magV < ROOTVSMALL)
         {
             info.setPoint(point2_);
@@ -92,14 +92,57 @@ Foam::pointIndexHit Foam::searchableCylinder::findNearest
     }
     else
     {
+        // inbetween endcaps. Might either be nearer endcaps or cylinder wall
+
         if (magV < ROOTVSMALL)
         {
-            info.setPoint(point1_ + parallel*unitDir_);
+            if (parallel < 0.5*magDir_)
+            {
+                info.setPoint(point1_);
+            }
+            else
+            {
+                info.setPoint(point2_);
+            }
         }
         else
         {
-            // Project onto radius
-            info.setPoint(point1_ + parallel*unitDir_ + radius_*v/magV);
+            vector n = v/magV;
+
+            // distance to endpoint: parallel or parallel-magDir
+            // distance to cylinder wall: magV-radius_
+
+            // Nearest cylinder point
+            point cylPt = sample + (radius_-magV)*n;
+
+            if (parallel < 0.5*magDir_)
+            {
+                // Project onto p1 endcap
+                point end1Pt = point1_ + min(magV, radius_)*n;
+
+                if (magSqr(sample-cylPt) < magSqr(sample-end1Pt))
+                {
+                    info.setPoint(cylPt);
+                }
+                else
+                {
+                    info.setPoint(end1Pt);
+                }
+            }
+            else
+            {
+                // Project onto p2 endcap
+                point end2Pt = point2_ + min(magV, radius_)*n;
+
+                if (magSqr(sample-cylPt) < magSqr(sample-end2Pt))
+                {
+                    info.setPoint(cylPt);
+                }
+                else
+                {
+                    info.setPoint(end2Pt);
+                }
+            }
         }
     }
 
-- 
2.11.4.GIT