ENH: moved pointInTet() to inside() method in tetrahedron

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -174,10 +174,10 @@ public:
             ) const;
 
             //- Return nearest point to p on tetrahedron
-            inline pointHit nearestPoint
-            (
-                const point& p
-            ) const;
+            inline pointHit nearestPoint(const point& p) const;
+
+            //- Return true if point is inside tetrahedron
+            inline bool inside(const point& pt) const;
 
             //- Return (min)containment sphere, i.e. the smallest sphere with
             //  all points inside. Returns pointHit with:

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -380,6 +380,38 @@ inline Foam::pointHit Foam::tetrahedron<Point, PointRef>::nearestPoint
 }
 
 
+template<class Point, class PointRef>
+bool Foam::tetrahedron<Point, PointRef>::inside(const point& pt) const
+{
+    // For robustness, assuming that the point is in the tet unless
+    // "definitively" shown otherwise by obtaining a positive dot
+    // product greater than a tolerance of SMALL.
+
+    // The tet is defined: tet(Cc, tetBasePt, pA, pB) where the normal
+    // vectors and base points for the half-space planes are:
+    // area[0] = Sa();
+    // area[1] = Sb();
+    // area[2] = Sc();
+    // area[3] = Sd();
+    // planeBase[0] = tetBasePt = tet.b()
+    // planeBase[1] = ptA       = tet.c()
+    // planeBase[2] = tetBasePt = tet.b()
+    // planeBase[3] = tetBasePt = tet.b()
+
+    vector n = Sa();
+    n /= (Foam::mag(n) + VSMALL);
+
+    if (((pt - b_) & n) > SMALL)
+    {
+        return false;
+    }
+    else
+    {
+        return true;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class Point, class PointRef>