From 123f1a4fda4965a00dba61fdfcdfb20844cd0af2 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 4 Jun 2010 10:52:48 +0100
Subject: [PATCH] BUG: searchableCylinder : points on axis not handled
 correctly

---
 .../searchableSurface/searchableCylinder.C         | 82 ++++++++--------------
 1 file changed, 30 insertions(+), 52 deletions(-)

diff --git a/src/meshTools/searchableSurface/searchableCylinder.C b/src/meshTools/searchableSurface/searchableCylinder.C
index 70d81e90..71043688 100644
--- a/src/meshTools/searchableSurface/searchableCylinder.C
+++ b/src/meshTools/searchableSurface/searchableCylinder.C
@@ -61,87 +61,65 @@ Foam::pointIndexHit Foam::searchableCylinder::findNearest
     // Decompose sample-point1 into normal and parallel component
     scalar parallel = (v & unitDir_);
 
-    // Remove the parallel component
+    // Remove the parallel component and normalise
     v -= parallel*unitDir_;
     scalar magV = mag(v);
 
+    if (magV < ROOTVSMALL)
+    {
+        v = vector::zero;
+    }
+    else
+    {
+        v /= magV;
+    }
 
     if (parallel <= 0)
     {
         // nearest is at point1 end cap. Clip to radius.
-        if (magV < ROOTVSMALL)
-        {
-            info.setPoint(point1_);
-        }
-        else
-        {
-            info.setPoint(point1_ + min(magV, radius_)*v/magV);
-        }
+        info.setPoint(point1_ + min(magV, radius_)*v);
     }
     else if (parallel >= magDir_)
     {
-        // nearest is at point2 end cap
-        if (magV < ROOTVSMALL)
-        {
-            info.setPoint(point2_);
-        }
-        else
-        {
-            info.setPoint(point2_ + min(magV, radius_)*v/magV);
-        }
+        // nearest is at point2 end cap. Clip to radius.
+        info.setPoint(point2_ + min(magV, radius_)*v);
     }
     else
     {
         // inbetween endcaps. Might either be nearer endcaps or cylinder wall
 
-        if (magV < ROOTVSMALL)
+        // distance to endpoint: parallel or parallel-magDir
+        // distance to cylinder wall: magV-radius_
+
+        // Nearest cylinder point
+        point cylPt = sample + (radius_-magV)*v;
+
+        if (parallel < 0.5*magDir_)
         {
-            if (parallel < 0.5*magDir_)
+            // Project onto p1 endcap
+            point end1Pt = point1_ + min(magV, radius_)*v;
+
+            if (magSqr(sample-cylPt) < magSqr(sample-end1Pt))
             {
-                info.setPoint(point1_);
+                info.setPoint(cylPt);
             }
             else
             {
-                info.setPoint(point2_);
+                info.setPoint(end1Pt);
             }
         }
         else
         {
-            vector n = v/magV;
-
-            // distance to endpoint: parallel or parallel-magDir
-            // distance to cylinder wall: magV-radius_
+            // Project onto p2 endcap
+            point end2Pt = point2_ + min(magV, radius_)*v;
 
-            // Nearest cylinder point
-            point cylPt = sample + (radius_-magV)*n;
-
-            if (parallel < 0.5*magDir_)
+            if (magSqr(sample-cylPt) < magSqr(sample-end2Pt))
             {
-                // 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);
-                }
+                info.setPoint(cylPt);
             }
             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);
-                }
+                info.setPoint(end2Pt);
             }
         }
     }
-- 
2.11.4.GIT