ENH: triangle: added nearest to line

src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,6 +42,7 @@ SourceFiles
 #include "Random.H"
 #include "FixedList.H"
 #include "UList.H"
+#include "linePointRef.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -222,6 +223,15 @@ public:
                 label& nearLabel
             ) const;
 
+            //- Return nearest point to line on triangle. Returns hit if
+            //  point is inside triangle. Sets edgePoint to point on edge
+            //  (hit if nearest is inside line)
+            inline pointHit nearestPoint
+            (
+                const linePointRef& edge,
+                pointHit& edgePoint
+            ) const;
+
 
     // IOstream operators
 

src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -659,6 +659,132 @@ inline bool Foam::triangle<Point, PointRef>::classify
 }
 
 
+template<class Point, class PointRef>
+inline Foam::pointHit Foam::triangle<Point, PointRef>::nearestPoint
+(
+    const linePointRef& ln,
+    pointHit& lnInfo
+) const
+{
+    vector q = ln.vec();
+    pointHit triInfo
+    (
+        triangle<Point, PointRef>::intersection
+        (
+            ln.start(),
+            q,
+            intersection::FULL_RAY
+        )
+    );
+
+    if (triInfo.hit())
+    {
+        // Line hits triangle. Find point on line.
+        if (triInfo.distance() > 1)
+        {
+            // Hit beyond endpoint
+            lnInfo.setMiss(true);
+            lnInfo.setPoint(ln.end());
+            scalar dist = Foam::mag(triInfo.hitPoint()-lnInfo.missPoint());
+            lnInfo.setDistance(dist);
+            triInfo.setMiss(true);
+            triInfo.setDistance(dist);
+        }
+        else if (triInfo.distance() < 0)
+        {
+            // Hit beyond startpoint
+            lnInfo.setMiss(true);
+            lnInfo.setPoint(ln.start());
+            scalar dist = Foam::mag(triInfo.hitPoint()-lnInfo.missPoint());
+            lnInfo.setDistance(dist);
+            triInfo.setMiss(true);
+            triInfo.setDistance(dist);
+        }
+        else
+        {
+            // Hit on line
+            lnInfo.setHit();
+            lnInfo.setPoint(triInfo.hitPoint());
+            lnInfo.setDistance(0.0);
+            triInfo.setDistance(0.0);
+        }
+    }
+    else
+    {
+        // Line skips triangle. See which triangle edge it gets closest to
+
+        point nearestEdgePoint;
+        point nearestLinePoint;
+        label minEdgeIndex = 0;
+        scalar minDist = ln.nearestDist
+        (
+            linePointRef(a_, b_),
+            nearestLinePoint,
+            nearestEdgePoint
+        );
+
+        {
+            point linePoint;
+            point triEdgePoint;
+            scalar dist = ln.nearestDist
+            (
+                linePointRef(b_, c_),
+                linePoint,
+                triEdgePoint
+            );
+            if (dist < minDist)
+            {
+                minDist = dist;
+                nearestEdgePoint = triEdgePoint;
+                nearestLinePoint = linePoint;
+                minEdgeIndex = 1;
+            }
+        }
+
+        {
+            point linePoint;
+            point triEdgePoint;
+            scalar dist = ln.nearestDist
+            (
+                linePointRef(c_, a_),
+                linePoint,
+                triEdgePoint
+            );
+            if (dist < minDist)
+            {
+                minDist = dist;
+                nearestEdgePoint = triEdgePoint;
+                nearestLinePoint = linePoint;
+                minEdgeIndex = 2;
+            }
+        }
+
+        lnInfo.setDistance(minDist);
+        triInfo.setDistance(minDist);
+        triInfo.setMiss(false);
+        triInfo.setPoint(nearestEdgePoint);
+
+        // Convert point on line to pointHit
+        if (Foam::mag(nearestLinePoint-ln.start()) < SMALL)
+        {
+            lnInfo.setMiss(true);
+            lnInfo.setPoint(ln.start());
+        }
+        else if (Foam::mag(nearestLinePoint-ln.end()) < SMALL)
+        {
+            lnInfo.setMiss(true);
+            lnInfo.setPoint(ln.end());
+        }
+        else
+        {
+            lnInfo.setHit();
+            lnInfo.setPoint(nearestLinePoint);
+        }
+    }
+    return triInfo;
+}
+
+
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class Point, class PointRef>