From: Michael Orlitzky <michael@orlitzky.com>
Date: Sun, 6 Nov 2011 23:16:54 +0000 (-0500)
Subject: Drop the no-longer-useful ThreeDimensional class.
X-Git-Tag: 0.0.1~71
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=5f01596d42cca3ec2b8236d697adb468cfcdb055;p=spline3.git

Drop the no-longer-useful ThreeDimensional class.
---

diff --git a/src/Cube.hs b/src/Cube.hs
index 9f93143..d863c29 100644
--- a/src/Cube.hs
+++ b/src/Cube.hs
@@ -24,12 +24,11 @@ import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
 
 import Cardinal
 import Comparisons ((~=), (~~=))
-import qualified Face (Face(Face, v0, v1, v2, v3))
+import qualified Face (Face(..), center)
 import FunctionValues (FunctionValues, eval, rotate)
 import Misc (all_equal, disjoint)
 import Point (Point(..), dot)
-import Tetrahedron (Tetrahedron(..), c, volume)
-import ThreeDimensional
+import Tetrahedron (Tetrahedron(..), barycenter, c, volume)
 
 data Cube = Cube { h  :: !Double,
                    i  :: !Int,
@@ -122,30 +121,20 @@ zmax cube = (k' + 1/2)*delta
       k' = fromIntegral (k cube) :: Double
       delta = h cube
 
-instance ThreeDimensional Cube where
-    -- | The center of Cube_ijk coincides with v_ijk at
-    --   (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
-    center cube = Point x y z
-           where
-             delta = h cube
-             i' = fromIntegral (i cube) :: Double
-             j' = fromIntegral (j cube) :: Double
-             k' = fromIntegral (k cube) :: Double
-             x = delta * i'
-             y = delta * j'
-             z = delta * k'
-
-    -- | It's easy to tell if a point is within a cube; just make sure
-    --   that it falls on the proper side of each of the cube's faces.
-    contains_point cube (Point x y z)
-        | x < (xmin cube) = False
-        | x > (xmax cube) = False
-        | y < (ymin cube) = False
-        | y > (ymax cube) = False
-        | z < (zmin cube) = False
-        | z > (zmax cube) = False
-        | otherwise = True
 
+-- | The center of Cube_ijk coincides with v_ijk at
+--   (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
+center :: Cube -> Point
+center cube =
+  Point x y z
+  where
+    delta = h cube
+    i' = fromIntegral (i cube) :: Double
+    j' = fromIntegral (j cube) :: Double
+    k' = fromIntegral (k cube) :: Double
+    x = delta * i'
+    y = delta * j'
+    z = delta * k'
 
 
 -- Face stuff.
@@ -230,7 +219,7 @@ tetrahedron cube 0 =
     where
       v0' = center cube
       ff  = front_face cube
-      v1' = center ff
+      v1' = Face.center ff
       v2' = Face.v0 ff
       v3' = Face.v1 ff
       vol = tetrahedra_volume cube
@@ -240,7 +229,7 @@ tetrahedron cube 1 =
     where
       v0' = center cube
       ff  = front_face cube
-      v1' = center ff
+      v1' = Face.center ff
       v2' = Face.v1 ff
       v3' = Face.v2 ff
       fv' = rotate ccwx (fv cube)
@@ -251,7 +240,7 @@ tetrahedron cube 2 =
     where
       v0' = center cube
       ff  = front_face cube
-      v1' = center ff
+      v1' = Face.center ff
       v2' = Face.v2 ff
       v3' = Face.v3 ff
       fv' = rotate ccwx $ rotate ccwx $ fv cube
@@ -262,7 +251,7 @@ tetrahedron cube 3 =
     where
       v0' = center cube
       ff  = front_face cube
-      v1' = center ff
+      v1' = Face.center ff
       v2' = Face.v3 ff
       v3' = Face.v0 ff
       fv' = rotate cwx (fv cube)
@@ -273,7 +262,7 @@ tetrahedron cube 4 =
     where
       v0' = center cube
       tf  = top_face cube
-      v1' = center tf
+      v1' = Face.center tf
       v2' = Face.v0 tf
       v3' = Face.v1 tf
       fv' = rotate cwy (fv cube)
@@ -284,7 +273,7 @@ tetrahedron cube 5 =
     where
       v0' = center cube
       tf  = top_face cube
-      v1' = center tf
+      v1' = Face.center tf
       v2' = Face.v1 tf
       v3' = Face.v2 tf
       fv' = rotate cwy $ rotate cwz $ fv cube
@@ -295,7 +284,7 @@ tetrahedron cube 6 =
     where
       v0' = center cube
       tf  = top_face cube
-      v1' = center tf
+      v1' = Face.center tf
       v2' = Face.v2 tf
       v3' = Face.v3 tf
       fv' = rotate cwy $ rotate cwz
@@ -308,7 +297,7 @@ tetrahedron cube 7 =
     where
       v0' = center cube
       tf  = top_face cube
-      v1' = center tf
+      v1' = Face.center tf
       v2' = Face.v3 tf
       v3' = Face.v0 tf
       fv' = rotate cwy $ rotate ccwz $ fv cube
@@ -319,7 +308,7 @@ tetrahedron cube 8 =
     where
       v0' = center cube
       bf  = back_face cube
-      v1' = center bf
+      v1' = Face.center bf
       v2' = Face.v0 bf
       v3' = Face.v1 bf
       fv' = rotate cwy $ rotate cwy $ fv cube
@@ -330,7 +319,7 @@ tetrahedron cube 9 =
     where
       v0' = center cube
       bf  = back_face cube
-      v1' = center bf
+      v1' = Face.center bf
       v2' = Face.v1 bf
       v3' = Face.v2 bf
       fv' = rotate cwy $ rotate cwy
@@ -343,7 +332,7 @@ tetrahedron cube 10 =
     where
       v0' = center cube
       bf  = back_face cube
-      v1' = center bf
+      v1' = Face.center bf
       v2' = Face.v2 bf
       v3' = Face.v3 bf
       fv' = rotate cwy $ rotate cwy
@@ -358,7 +347,7 @@ tetrahedron cube 11 =
     where
       v0' = center cube
       bf  = back_face cube
-      v1' = center bf
+      v1' = Face.center bf
       v2' = Face.v3 bf
       v3' = Face.v0 bf
       fv' = rotate cwy $ rotate cwy
@@ -371,7 +360,7 @@ tetrahedron cube 12 =
     where
       v0' = center cube
       df  = down_face cube
-      v1' = center df
+      v1' = Face.center df
       v2' = Face.v0 df
       v3' = Face.v1 df
       fv' = rotate ccwy $ fv cube
@@ -382,7 +371,7 @@ tetrahedron cube 13 =
     where
       v0' = center cube
       df  = down_face cube
-      v1' = center df
+      v1' = Face.center df
       v2' = Face.v1 df
       v3' = Face.v2 df
       fv' = rotate ccwy $ rotate ccwz $ fv cube
@@ -393,7 +382,7 @@ tetrahedron cube 14 =
     where
       v0' = center cube
       df  = down_face cube
-      v1' = center df
+      v1' = Face.center df
       v2' = Face.v2 df
       v3' = Face.v3 df
       fv' = rotate ccwy $ rotate ccwz
@@ -406,7 +395,7 @@ tetrahedron cube 15 =
     where
       v0' = center cube
       df  = down_face cube
-      v1' = center df
+      v1' = Face.center df
       v2' = Face.v3 df
       v3' = Face.v0 df
       fv' = rotate ccwy $ rotate cwz $ fv cube
@@ -417,7 +406,7 @@ tetrahedron cube 16 =
     where
       v0' = center cube
       rf  = right_face cube
-      v1' = center rf
+      v1' = Face.center rf
       v2' = Face.v0 rf
       v3' = Face.v1 rf
       fv' = rotate ccwz $ fv cube
@@ -428,7 +417,7 @@ tetrahedron cube 17 =
     where
       v0' = center cube
       rf  = right_face cube
-      v1' = center rf
+      v1' = Face.center rf
       v2' = Face.v1 rf
       v3' = Face.v2 rf
       fv' = rotate ccwz $ rotate cwy $ fv cube
@@ -439,7 +428,7 @@ tetrahedron cube 18 =
     where
       v0' = center cube
       rf  = right_face cube
-      v1' = center rf
+      v1' = Face.center rf
       v2' = Face.v2 rf
       v3' = Face.v3 rf
       fv' = rotate ccwz $ rotate cwy
@@ -452,7 +441,7 @@ tetrahedron cube 19 =
     where
       v0' = center cube
       rf  = right_face cube
-      v1' = center rf
+      v1' = Face.center rf
       v2' = Face.v3 rf
       v3' = Face.v0 rf
       fv' = rotate ccwz $ rotate ccwy
@@ -464,7 +453,7 @@ tetrahedron cube 20 =
     where
       v0' = center cube
       lf  = left_face cube
-      v1' = center lf
+      v1' = Face.center lf
       v2' = Face.v0 lf
       v3' = Face.v1 lf
       fv' = rotate cwz $ fv cube
@@ -475,7 +464,7 @@ tetrahedron cube 21 =
     where
       v0' = center cube
       lf  = left_face cube
-      v1' = center lf
+      v1' = Face.center lf
       v2' = Face.v1 lf
       v3' = Face.v2 lf
       fv' = rotate cwz $ rotate ccwy $ fv cube
@@ -486,7 +475,7 @@ tetrahedron cube 22 =
     where
       v0' = center cube
       lf  = left_face cube
-      v1' = center lf
+      v1' = Face.center lf
       v2' = Face.v2 lf
       v3' = Face.v3 lf
       fv' = rotate cwz $ rotate ccwy
@@ -499,7 +488,7 @@ tetrahedron cube 23 =
     where
       v0' = center cube
       lf  = left_face cube
-      v1' = center lf
+      v1' = Face.center lf
       v2' = Face.v3 lf
       v3' = Face.v0 lf
       fv' = rotate cwz $ rotate cwy
@@ -626,6 +615,7 @@ find_containing_tetrahedron cube p =
     top_half = in_top_half cube p
     left_half = in_left_half cube p
 
+    candidates :: V.Vector Tetrahedron
     candidates =
       if front_half then
 
@@ -656,10 +646,20 @@ find_containing_tetrahedron cube p =
     -- Use the dot product instead of Euclidean distance here to save
     -- a sqrt(). So, "distances" below really means "distances
     -- squared."
-    distances = V.map ((dot p) . center) candidates
+    distances :: V.Vector Double
+    distances = V.map ((dot p) . barycenter) candidates
+
+    shortest_distance :: Double
     shortest_distance = V.minimum distances
+
+    -- Compute the index of the tetrahedron with the center closest to
+    -- p. This is a bad algorithm, but don't change it! If you make it
+    -- smarter by finding the index of shortest_distance in distances
+    -- (this should give the same answer and avoids recomputing the
+    -- dot product), the program gets slower. Seriously!
+    lucky_idx :: Maybe Int
     lucky_idx = V.findIndex
-                  (\t -> (center t) `dot` p == shortest_distance)
+                  (\t -> (barycenter t) `dot` p == shortest_distance)
                   candidates
 
 
@@ -716,7 +716,8 @@ prop_all_volumes_exact cube =
     where
       delta = h cube
 
--- | All tetrahedron should have their v0 located at the center of the cube.
+-- | All tetrahedron should have their v0 located at the center of the
+--   cube.
 prop_v0_all_equal :: Cube -> Bool
 prop_v0_all_equal cube = (v0 t0) == (v0 t1)
     where
diff --git a/src/Face.hs b/src/Face.hs
index bc316a7..785acaf 100644
--- a/src/Face.hs
+++ b/src/Face.hs
@@ -1,11 +1,12 @@
 -- | The Face module just contains the definition of the 'Face' data
 --   type and its two typeclass instances.
-module Face
-  ( Face(..) )
+module Face (
+  Face(..),
+  center
+  )
 where
 
 import Point
-import ThreeDimensional
 
 data Face = Face { v0 :: !Point,
                    v1 :: !Point,
@@ -14,22 +15,17 @@ data Face = Face { v0 :: !Point,
           deriving (Eq)
 
 instance Show Face where
-    show f = "Face:\n" ++
-             "  v0: " ++ (show (v0 f)) ++ "\n" ++
-             "  v1: " ++ (show (v1 f)) ++ "\n" ++
-             "  v2: " ++ (show (v2 f)) ++ "\n" ++
-             "  v3: " ++ (show (v3 f)) ++ "\n"
+    show (Face v0' v1' v2' v3') =
+      "Face:\n" ++
+        "  v0: " ++ (show v0') ++ "\n" ++
+        "  v1: " ++ (show v1') ++ "\n" ++
+        "  v2: " ++ (show v2') ++ "\n" ++
+        "  v3: " ++ (show v3') ++ "\n"
 
--- | The 'Face' type is an instance of 'ThreeDimensional' so that we
---   can call the 'center' function on it. This is useful because the
---   center of a face is always a vertex of a tetrahedron.
-instance ThreeDimensional Face where
-    -- | Since a face is square, we can just average the four vertices
-    --   to find the center.
-    center f = ((v0 f) + (v1 f) + (v2 f) + (v3 f)) `scale` (1/4)
-
-    -- | It's possible to implement this, but it hasn't been done
-    --   yet. A face will contain a point if the point lies in the same
-    --   plane as the vertices of the face, and if it falls on the
-    --   correct side of the four sides of the face.
-    contains_point _ _ = False
+-- | Returns the center of the given face. Since a face is just
+--   square, we can average the four vertices to find its center. This
+--   is useful because the center of a face is always a vertex of a
+--   tetrahedron.
+center :: Face -> Point
+center (Face v0' v1' v2' v3') =
+  (v0' + v1' + v2' + v3') `scale` (1/4)
diff --git a/src/Grid.hs b/src/Grid.hs
index a1058d0..94b9e1d 100644
--- a/src/Grid.hs
+++ b/src/Grid.hs
@@ -32,8 +32,12 @@ import Examples (trilinear, trilinear9x9x9, zeros, naturals_1d)
 import FunctionValues (make_values, value_at)
 import Point (Point(..))
 import ScaleFactor (ScaleFactor)
-import Tetrahedron (Tetrahedron, c, polynomial, v0, v1, v2, v3)
-import ThreeDimensional (ThreeDimensional(..))
+import Tetrahedron (
+  Tetrahedron(v0,v1,v2,v3),
+  c,
+  contains_point,
+  polynomial,
+  )
 import Values (Values3D, dims, empty3d, zoom_shape)
 
 
diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs
index a6d69a2..7595729 100644
--- 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 }
-
-        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 }
-
-        b3_unscaled :: Double
-        b3_unscaled = volume inner_tetrahedron
-          where inner_tetrahedron = t { v3 = p0 }
+        inner_tetrahedron = t { v0 = 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 }
+
+    b3_unscaled :: Double
+    b3_unscaled = volume inner_tetrahedron
+      where inner_tetrahedron = t { v3 = p0 }
+
 
 {-# INLINE polynomial #-}
 polynomial :: Tetrahedron -> (RealFunction Point)
diff --git a/src/ThreeDimensional.hs b/src/ThreeDimensional.hs
deleted file mode 100644
index 4b74f18..0000000
--- a/src/ThreeDimensional.hs
+++ /dev/null
@@ -1,9 +0,0 @@
-module ThreeDimensional
-  ( ThreeDimensional(..) )
-where
-
-import Point(Point)
-
-class ThreeDimensional a where
-    center :: a -> Point
-    contains_point :: a -> Point -> Bool