BUG: searchableCylinder nearest not correct for endcaps
authormattijs <mattijs>
Thu, 3 Jun 2010 16:44:19 +0000 (3 17:44 +0100)
committermattijs <mattijs>
Thu, 3 Jun 2010 16:44:19 +0000 (3 17:44 +0100)
src/meshTools/searchableSurface/searchableCylinder.C

index 0df222c..70d81e9 100644 (file)
@@ -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);
+                }
+            }
         }
     }