Remove the coincident vertices guards from the Tetrahedron volume function.

[spline3.git] / src / Tetrahedron.hs

diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs
index a6d69a2c68af2c8051d3fab76b9bcfd4e3201bfd..2a4227a7f1bbee2d10a99baa66a07904e1e16c15 100644 (file)
--- a/src/Tetrahedron.hs
+++ b/src/Tetrahedron.hs
@@ -5,7 +5,9 @@ module Tetrahedron (
   b1, -- Cube test
   b2, -- Cube test
   b3, -- Cube test
+  barycenter,
   c,
+  contains_point,
   polynomial,
   tetrahedron_properties,
   tetrahedron_tests,
@@ -30,7 +32,6 @@ import FunctionValues (FunctionValues(..), empty_values)
 import Misc (factorial)
 import Point (Point(..), scale)
 import RealFunction (RealFunction, cmult, fexp)
-import ThreeDimensional (ThreeDimensional(..))
 
 data Tetrahedron =
   Tetrahedron { function_values :: FunctionValues,
@@ -67,34 +68,40 @@ instance Show Tetrahedron where
              "  v3: " ++ (show (v3 t)) ++ "\n"
 
 
-instance ThreeDimensional Tetrahedron where
-    center (Tetrahedron _ v0' v1' v2' v3' _) =
-        (v0' + v1' + v2' + v3') `scale` (1/4)
-
-    -- contains_point is only used in tests.
-    contains_point t p0 =
-      b0_unscaled `nearly_ge` 0 &&
-      b1_unscaled `nearly_ge` 0 &&
-      b2_unscaled `nearly_ge` 0 &&
-      b3_unscaled `nearly_ge` 0
+-- | Find the barycenter of the given tetrahedron.
+--   We just average the four vertices.
+barycenter :: Tetrahedron -> Point
+barycenter (Tetrahedron _ v0' v1' v2' v3' _) =
+  (v0' + v1' + v2' + v3') `scale` (1/4)
+
+-- | A point is internal to a tetrahedron if all of its barycentric
+--   coordinates with respect to that tetrahedron are non-negative.
+contains_point :: Tetrahedron -> Point -> Bool
+contains_point t p0 =
+  b0_unscaled `nearly_ge` 0 &&
+  b1_unscaled `nearly_ge` 0 &&
+  b2_unscaled `nearly_ge` 0 &&
+  b3_unscaled `nearly_ge` 0
+  where
+    -- Drop the useless division and volume calculation that we
+    -- would do if we used the regular b0,..b3 functions.
+    b0_unscaled :: Double
+    b0_unscaled = volume inner_tetrahedron
       where
-        -- Drop the useless division and volume calculation that we
-        -- would do if we used the regular b0,..b3 functions.
-        b0_unscaled :: Double
-        b0_unscaled = volume inner_tetrahedron
-          where inner_tetrahedron = t { v0 = p0 }
+        inner_tetrahedron = t { v0 = p0 }
+
+    b1_unscaled :: Double
+    b1_unscaled = volume inner_tetrahedron
+      where inner_tetrahedron = t { v1 = p0 }
 
-        b1_unscaled :: Double
-        b1_unscaled = volume inner_tetrahedron
-          where inner_tetrahedron = t { v1 = p0 }
+    b2_unscaled :: Double
+    b2_unscaled = volume inner_tetrahedron
+      where inner_tetrahedron = t { v2 = p0 }
 
-        b2_unscaled :: Double
-        b2_unscaled = volume inner_tetrahedron
-          where inner_tetrahedron = t { v2 = p0 }
+    b3_unscaled :: Double
+    b3_unscaled = volume inner_tetrahedron
+      where inner_tetrahedron = t { v3 = p0 }
 
-        b3_unscaled :: Double
-        b3_unscaled = volume inner_tetrahedron
-          where inner_tetrahedron = t { v3 = p0 }
 
 {-# INLINE polynomial #-}
 polynomial :: Tetrahedron -> (RealFunction Point)
@@ -306,20 +313,8 @@ det p0 p1 p2 p3 =
 --   page 436.
 {-# INLINE volume #-}
 volume :: Tetrahedron -> Double
-volume t
-       | v0' == v1' = 0
-       | v0' == v2' = 0
-       | v0' == v3' = 0
-       | v1' == v2' = 0
-       | v1' == v3' = 0
-       | v2' == v3' = 0
-       | otherwise = (1/6)*(det v0' v1' v2' v3')
-  where
-    v0' = v0 t
-    v1' = v1 t
-    v2' = v2 t
-    v3' = v3 t
-
+volume (Tetrahedron _ v0' v1' v2' v3' _) =
+  (1/6)*(det v0' v1' v2' v3')
 
 -- | The barycentric coordinates of a point with respect to v0.
 {-# INLINE b0 #-}