From: Michael Orlitzky <michael@orlitzky.com>
Date: Wed, 4 May 2011 00:57:37 +0000 (-0400)
Subject: Begin overhauling the program to handle other tetrahedra. Main is
X-Git-Tag: 0.0.1~343
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=6fb9ab6b6068870323e996da931fc04c7710e3e4;p=spline3.git

Begin overhauling the program to handle other tetrahedra. Main is
pretty much disabled at this point, waiting for me to implement the
tetrahedron[0-23] functions.
---

diff --git a/src/Cardinal.hs b/src/Cardinal.hs
index e4e9464..f62896f 100644
--- a/src/Cardinal.hs
+++ b/src/Cardinal.hs
@@ -5,19 +5,46 @@ data Cardinal = F
               | B
               | L
               | R
-              | T
               | D
+              | T
+              | FL
+              | FR
+              | FD
+              | FT
+              | BL
+              | BR
+              | BD
+              | BT
+              | LD
+              | LT
+              | RD
+              | RT
+              | FLD
+              | FLT
+              | FRD
+              | FRT
+              | BLD
+              | BLT
+              | BRD
+              | BRT
+              | I
+              | Scalar Double
               | Sum Cardinal Cardinal
               | Difference Cardinal Cardinal
               | Product Cardinal Cardinal
-              | ScalarProduct Double Cardinal
+              | Quotient Cardinal Cardinal
                 deriving (Show, Eq)
 
 instance Num Cardinal where
     x + y = Sum x y
     x - y = Difference x y
-    x * y = Sum x y
-    negate x = ScalarProduct (-1) x
+    x * y = Product x y
+    negate x = Product (Scalar (-1)) x
     abs x = x
     signum x = x
-    fromInteger _ = F -- Whatever.
+    fromInteger x = Scalar (fromIntegral x)
+
+instance Fractional Cardinal where
+    x / y   = Quotient x y
+    recip x = Quotient (Scalar 1) x
+    fromRational x = Scalar (fromRational x)
diff --git a/src/Cube.hs b/src/Cube.hs
index 3011281..23cbe4d 100644
--- a/src/Cube.hs
+++ b/src/Cube.hs
@@ -1,59 +1,42 @@
 module Cube
 where
 
-import Grid
+import Face
+import FunctionValues
+--import Grid
 import Point
 import ThreeDimensional
 
-class Gridded a where
-    back :: a -> Cube
-    down :: a -> Cube
-    front :: a -> Cube
-    left :: a -> Cube
-    right :: a -> Cube
-    top :: a -> Cube
-
-
-data Cube = Cube { grid :: Grid,
+data Cube = Cube { h :: Double,
                    i :: Int,
                    j :: Int,
                    k :: Int,
-                   datum :: Double }
+                   fv :: FunctionValues }
             deriving (Eq)
 
 
 instance Show Cube where
     show c =
         "Cube_" ++ (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k c)) ++
-        " (Grid: " ++ (show (grid c)) ++ ")" ++
         " (Center: " ++ (show (center c)) ++ ")" ++
         " (xmin: " ++ (show (xmin c)) ++ ")" ++
         " (xmax: " ++ (show (xmax c)) ++ ")" ++
         " (ymin: " ++ (show (ymin c)) ++ ")" ++
         " (ymax: " ++ (show (ymax c)) ++ ")" ++
         " (zmin: " ++ (show (zmin c)) ++ ")" ++
-        " (zmax: " ++ (show (zmax c)) ++ ")" ++
-        " (datum: " ++ (show (datum c)) ++ ")\n\n"
+        " (zmax: " ++ (show (zmax c)) ++ ")"
 
 empty_cube :: Cube
-empty_cube = Cube empty_grid 0 0 0 0
+empty_cube = Cube 0 0 0 0 empty_values
 
 
-instance Gridded Cube where
-    back c = cube_at (grid c) ((i c) + 1) (j c) (k c)
-    down c = cube_at (grid c) (i c) (j c) ((k c) - 1)
-    front c = cube_at (grid c) ((i c) - 1) (j c) (k c)
-    left c = cube_at (grid c) (i c) ((j c) - 1) (k c)
-    right c = cube_at (grid c) (i c) ((j c) + 1) (k c)
-    top c = cube_at (grid c) (i c) (j c) ((k c) + 1)
-
 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
 --   p. 76.
 xmin :: Cube -> Double
 xmin c = (2*i' - 1)*delta / 2
     where
       i' = fromIntegral (i c) :: Double
-      delta = h (grid c)
+      delta = h c
 
 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
 --   p. 76.
@@ -61,7 +44,7 @@ xmax :: Cube -> Double
 xmax c = (2*i' + 1)*delta / 2
     where
       i' = fromIntegral (i c) :: Double
-      delta = h (grid c)
+      delta = h c
 
 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
 --   p. 76.
@@ -69,7 +52,7 @@ ymin :: Cube -> Double
 ymin c = (2*j' - 1)*delta / 2
     where
       j' = fromIntegral (j c) :: Double
-      delta = h (grid c)
+      delta = h c
 
 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
 --   p. 76.
@@ -77,7 +60,7 @@ ymax :: Cube -> Double
 ymax c = (2*j' + 1)*delta / 2
     where
       j' = fromIntegral (j c) :: Double
-      delta = h (grid c)
+      delta = h c
 
 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
 --   p. 76.
@@ -85,7 +68,7 @@ zmin :: Cube -> Double
 zmin c = (2*k' - 1)*delta / 2
     where
       k' = fromIntegral (k c) :: Double
-      delta = h (grid c)
+      delta = h c
 
 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
 --   p. 76.
@@ -93,14 +76,14 @@ zmax :: Cube -> Double
 zmax c = (2*k' + 1)*delta / 2
     where
       k' = fromIntegral (k c) :: Double
-      delta = h (grid c)
+      delta = h c
 
 instance ThreeDimensional Cube where
     -- | The center of Cube_ijk coincides with v_ijk at
     --   (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
     center c = (x, y, z)
            where
-             delta = h (grid c)
+             delta = h c
              i' = fromIntegral (i c) :: Double
              j' = fromIntegral (j c) :: Double
              k' = fromIntegral (k c) :: Double
@@ -118,73 +101,120 @@ instance ThreeDimensional Cube where
         | otherwise = True
 
 
-instance Num Cube where
-    (Cube g1 i1 j1 k1 d1) + (Cube _ i2 j2 k2 d2) =
-        Cube g1 (i1 + i2) (j1 + j2) (k1 + k2) (d1 + d2)
+-- instance Num Cube where
+--     (Cube g1 i1 j1 k1 d1) + (Cube _ i2 j2 k2 d2) =
+--         Cube g1 (i1 + i2) (j1 + j2) (k1 + k2) (d1 + d2)
 
-    (Cube g1 i1 j1 k1 d1) - (Cube _ i2 j2 k2 d2) =
-        Cube g1 (i1 - i2) (j1 - j2) (k1 - k2) (d1 - d2)
+--     (Cube g1 i1 j1 k1 d1) - (Cube _ i2 j2 k2 d2) =
+--         Cube g1 (i1 - i2) (j1 - j2) (k1 - k2) (d1 - d2)
 
-    (Cube g1 i1 j1 k1 d1) * (Cube _ i2 j2 k2 d2) =
-        Cube g1 (i1 * i2) (j1 * j2) (k1 * k2) (d1 * d2)
+--     (Cube g1 i1 j1 k1 d1) * (Cube _ i2 j2 k2 d2) =
+--         Cube g1 (i1 * i2) (j1 * j2) (k1 * k2) (d1 * d2)
 
-    abs (Cube g1 i1 j1 k1 d1) =
-        Cube g1 (abs i1) (abs j1) (abs k1) (abs d1)
+--     abs (Cube g1 i1 j1 k1 d1) =
+--         Cube g1 (abs i1) (abs j1) (abs k1) (abs d1)
 
-    signum (Cube g1 i1 j1 k1 d1) =
-        Cube g1 (signum i1) (signum j1) (signum k1) (signum d1)
+--     signum (Cube g1 i1 j1 k1 d1) =
+--         Cube g1 (signum i1) (signum j1) (signum k1) (signum d1)
 
-    fromInteger x = empty_cube { datum = (fromIntegral x) }
+--     fromInteger x = empty_cube { datum = (fromIntegral x) }
 
-instance Fractional Cube where
-    (Cube g1 i1 j1 k1 d1) / (Cube _ _ _ _ d2) =
-        Cube g1 i1 j1 k1 (d1 / d2)
+-- instance Fractional Cube where
+--     (Cube g1 i1 j1 k1 d1) / (Cube _ _ _ _ d2) =
+--         Cube g1 i1 j1 k1 (d1 / d2)
 
-    recip (Cube g1 i1 j1 k1 d1) =
-        Cube g1 i1 j1 k1 (recip d1)
+--     recip (Cube g1 i1 j1 k1 d1) =
+--         Cube g1 i1 j1 k1 (recip d1)
 
-    fromRational q = empty_cube { datum = fromRational q }
+--     fromRational q = empty_cube { datum = fromRational q }
 
--- | Constructs a cube, switching the x and z axes.
-reverse_cube :: Grid -> Int -> Int -> Int -> Double -> Cube
-reverse_cube g k' j' i' = Cube g i' j' k'
 
 
 -- | Return the cube corresponding to the grid point i,j,k. The list
 --   of cubes is stored as [z][y][x] but we'll be requesting it by
 --   [x][y][z] so we flip the indices in the last line.
-cube_at :: Grid -> Int -> Int -> Int -> Cube
-cube_at g i' j' k'
-        | i' >= length (function_values g) = Cube g i' j' k' 0
-        | i' < 0                          = Cube g i' j' k' 0
-        | j' >= length ((function_values g) !! i') = Cube g i' j' k' 0
-        | j' < 0                                  = Cube g i' j' k' 0
-        | k' >= length (((function_values g) !! i') !! j') = Cube g i' j' k' 0
-        | k' < 0                                          = Cube g i' j' k' 0
-        | otherwise =
-            (((cubes g) !! k') !! j') !! i'
-
-
--- These next three functions basically form a 'for' loop, looping
--- through the xs, ys, and zs in that order.
-
--- | The cubes_from_values function will return a function that takes
---   a list of values and returns a list of cubes. It could just as
---   well be written to take the values as a parameter; the omission
---   of the last parameter is known as an eta reduce.
-cubes_from_values :: Grid -> Int -> Int -> ([Double] -> [Cube])
-cubes_from_values g i' j' =
-    zipWith (reverse_cube g i' j') [0..]
-
--- | The cubes_from_planes function will return a function that takes
---   a list of planes and returns a list of cubes. It could just as
---   well be written to take the planes as a parameter; the omission
---   of the last parameter is known as an eta reduce.
-cubes_from_planes :: Grid -> Int -> ([[Double]] -> [[Cube]])
-cubes_from_planes g i' =
-    zipWith (cubes_from_values g i') [0..]
-
--- | Takes a grid as an argument, and returns a three-dimensional list
---   of cubes centered on its grid points.
-cubes :: Grid -> [[[Cube]]]
-cubes g = zipWith (cubes_from_planes g) [0..] (function_values g)
+-- cube_at :: Grid -> Int -> Int -> Int -> Cube
+-- cube_at g i' j' k'
+--         | i' >= length (function_values g) = Cube g i' j' k' 0
+--         | i' < 0                          = Cube g i' j' k' 0
+--         | j' >= length ((function_values g) !! i') = Cube g i' j' k' 0
+--         | j' < 0                                  = Cube g i' j' k' 0
+--         | k' >= length (((function_values g) !! i') !! j') = Cube g i' j' k' 0
+--         | k' < 0                                          = Cube g i' j' k' 0
+--         | otherwise =
+--             (((cubes g) !! k') !! j') !! i'
+
+
+
+
+
+
+-- Face stuff.
+
+-- | The top (in the direction of z) face of the cube.
+top_face :: Cube -> Face
+top_face c = Face v0' v1' v2' v3'
+    where
+      delta = (1/2)*(h c)
+      v0' = (center c) + (-delta, delta, delta)
+      v1' = (center c) + (delta, delta, delta)
+      v2' = (center c) + (delta, -delta, delta)
+      v3' = (center c) + (-delta, -delta, delta)
+
+
+
+-- | The back (in the direction of x) face of the cube.
+back_face :: Cube -> Face
+back_face c = Face v0' v1' v2' v3'
+    where
+      delta = (1/2)*(h c)
+      v0' = (center c) + (delta, delta, delta)
+      v1' = (center c) + (delta, delta, -delta)
+      v2' = (center c) + (delta, -delta, -delta)
+      v3' = (center c) + (delta, -delta, delta)
+
+
+-- The bottom face (in the direction of -z) of the cube.
+down_face :: Cube -> Face
+down_face c = Face v0' v1' v2' v3'
+    where
+      delta = (1/2)*(h c)
+      v0' = (center c) + (delta, delta, -delta)
+      v1' = (center c) + (-delta, delta, -delta)
+      v2' = (center c) + (-delta, -delta, -delta)
+      v3' = (center c) + (delta, -delta, -delta)
+
+
+
+-- | The front (in the direction of -x) face of the cube.
+front_face :: Cube -> Face
+front_face c = Face v0' v1' v2' v3'
+    where
+      delta = (1/2)*(h c)
+      v0' = (center c) + (-delta, delta, -delta)
+      v1' = (center c) + (-delta, delta, delta)
+      v2' = (center c) + (-delta, -delta, delta)
+      v3' = (center c) + (-delta, -delta, -delta)
+
+
+-- | The left (in the direction of -y) face of the cube.
+left_face :: Cube -> Face
+left_face c = Face v0' v1' v2' v3'
+    where
+      delta = (1/2)*(h c)
+      v0' = (center c) + (-delta, -delta, delta)
+      v1' = (center c) + (delta, -delta, delta)
+      v2' = (center c) + (delta, -delta, -delta)
+      v3' = (center c) + (-delta, -delta, -delta)
+
+
+-- | The right (in the direction of y) face of the cube.
+right_face :: Cube -> Face
+right_face c = Face v0' v1' v2' v3'
+    where
+      delta = (1/2)*(h c)
+      v0' = (center c) + (-delta, delta, -delta)
+      v1' = (center c) + (delta, delta, -delta)
+      v2' = (center c) + (delta, delta, delta)
+      v3' = (center c) + (-delta, delta, delta)
+
diff --git a/src/Face.hs b/src/Face.hs
index 6f7e99b..bc1dfee 100644
--- a/src/Face.hs
+++ b/src/Face.hs
@@ -1,160 +1,83 @@
 module Face
 where
 
-import Cube
-import Grid
 import Point
-import Tetrahedron hiding (c, cube, v0, v1, v2, v3)
 import ThreeDimensional
 
-data Face = Face { cube :: Cube,
-                   v0 :: Point,
+data Face = Face { v0 :: Point,
                    v1 :: Point,
                    v2 :: Point,
                    v3 :: Point }
           deriving (Eq)
 
 instance Show Face where
-    show f = "Face (Cube_" ++ (show i') ++ "," ++ (show j') ++ "," ++
-             (show k') ++ ") " ++ "(v0: " ++ (show (v0 f)) ++ ") (v1: " ++
-             (show (v1 f)) ++ ") (v2: " ++ (show (v2 f)) ++ ") (v3: " ++
-             (show (v3 f)) ++ ")\n\n"
-             where
-               i' = i (cube f)
-               j' = j (cube f)
-               k' = k (cube f)
+    show f = "Face:\n" ++
+             "  v0: " ++ (show (v0 f)) ++ "\n" ++
+             "  v1: " ++ (show (v1 f)) ++ "\n" ++
+             "  v2: " ++ (show (v2 f)) ++ "\n" ++
+             "  v3: " ++ (show (v3 f)) ++ "\n"
 
 instance ThreeDimensional Face where
     center f = ((v0 f) + (v1 f) + (v2 f) + (v3 f)) `scale` (1/4)
     -- Too lazy to implement this right now.
     contains_point _ _ = False
 
--- | The top (in the direction of z) face of the cube.
-face0 :: Cube -> Face
-face0 c = Face c v0' v1' v2' v3'
-    where
-      g = grid c
-      delta = (1/2)*(h g)
-      v0' = (center c) + (-delta, delta, delta)
-      v1' = (center c) + (delta, delta, delta)
-      v2' = (center c) + (delta, -delta, delta)
-      v3' = (center c) + (-delta, -delta, delta)
 
--- | The back (in the direction of x) face of the cube.
-face1 :: Cube -> Face
-face1 c = Face c v0' v1' v2' v3'
-    where
-      g = grid c
-      delta = (1/2)*(h g)
-      v0' = (center c) + (delta, delta, delta)
-      v1' = (center c) + (delta, delta, -delta)
-      v2' = (center c) + (delta, -delta, -delta)
-      v3' = (center c) + (delta, -delta, delta)
-
-
--- The bottom face (in the direction of -z) of the cube.
-face2 :: Cube -> Face
-face2 c = Face c v0' v1' v2' v3'
-    where
-      g = grid c
-      delta = (1/2)*(h g)
-      v0' = (center c) + (delta, delta, -delta)
-      v1' = (center c) + (-delta, delta, -delta)
-      v2' = (center c) + (-delta, -delta, -delta)
-      v3' = (center c) + (delta, -delta, -delta)
-
-
--- | The front (in the direction of -x) face of the cube.
-face3 :: Cube -> Face
-face3 c = Face c v0' v1' v2' v3'
-    where
-      g = grid c
-      delta = (1/2)*(h g)
-      v0' = (center c) + (-delta, delta, -delta)
-      v1' = (center c) + (-delta, delta, delta)
-      v2' = (center c) + (-delta, -delta, delta)
-      v3' = (center c) + (-delta, -delta, -delta)
-
-
--- | The left (in the direction of -y) face of the cube.
-face4 :: Cube -> Face
-face4 c = Face c v0' v1' v2' v3'
-    where
-      g = grid c
-      delta = (1/2)*(h g)
-      v0' = (center c) + (-delta, -delta, delta)
-      v1' = (center c) + (delta, -delta, delta)
-      v2' = (center c) + (delta, -delta, -delta)
-      v3' = (center c) + (-delta, -delta, -delta)
-
-
--- | The right (in the direction of y) face of the cube.
-face5 :: Cube -> Face
-face5 c = Face c v0' v1' v2' v3'
-    where
-      g = grid c
-      delta = (1/2)*(h g)
-      v0' = (center c) + (-delta, delta, -delta)
-      v1' = (center c) + (delta, delta, -delta)
-      v2' = (center c) + (delta, delta, delta)
-      v3' = (center c) + (-delta, delta, delta)
-
-
-tetrahedron0 :: Face -> Tetrahedron
-tetrahedron0 f =
-    Tetrahedron c v0' v1' v2' v3'
-    where
-      c = cube f
-      v0' = v0 f
-      v1' = v1 f
-      v2' = center f
-      v3' = center c
-
-tetrahedron1 :: Face -> Tetrahedron
-tetrahedron1 f =
-    Tetrahedron c v0' v1' v2' v3'
-    where
-      c = cube f
-      v0' = v1 f
-      v1' = v2 f
-      v2' = center f
-      v3' = center c
-
-
-tetrahedron2 :: Face -> Tetrahedron
-tetrahedron2 f =
-    Tetrahedron c v0' v1' v2' v3'
-    where
-      c = cube f
-      v0' = v2 f
-      v1' = v3 f
-      v2' = center f
-      v3' = center c
-
-
-tetrahedron3 :: Face -> Tetrahedron
-tetrahedron3 f =
-    Tetrahedron c v0' v1' v2' v3'
-    where
-      c = cube f
-      v0' = v3 f
-      v1' = v0 f
-      v2' = center f
-      v3' = center c
-
-tetrahedrons :: Cube -> [Tetrahedron]
-tetrahedrons c =
-    concat [
- [tetrahedron0 f0, tetrahedron1 f0, tetrahedron2 f0, tetrahedron3 f0],
- [tetrahedron0 f1, tetrahedron1 f1, tetrahedron2 f1, tetrahedron3 f2],
- [tetrahedron0 f2, tetrahedron1 f2, tetrahedron2 f2, tetrahedron3 f2],
- [tetrahedron0 f3, tetrahedron1 f3, tetrahedron2 f3, tetrahedron3 f3],
- [tetrahedron0 f4, tetrahedron1 f4, tetrahedron2 f4, tetrahedron3 f4],
- [tetrahedron0 f5, tetrahedron1 f5, tetrahedron2 f5, tetrahedron3 f5] ]
-    where
-      f0 = face0 c
-      f1 = face1 c
-      f2 = face2 c
-      f3 = face3 c
-      f4 = face4 c
-      f5 = face5 c
+-- tetrahedron0 :: Face -> Tetrahedron
+-- tetrahedron0 f =
+--     Tetrahedron c v0' v1' v2' v3'
+--     where
+--       c = cube f
+--       v0' = v0 f
+--       v1' = v1 f
+--       v2' = center f
+--       v3' = center c
+
+-- tetrahedron1 :: Face -> Tetrahedron
+-- tetrahedron1 f =
+--     Tetrahedron c v0' v1' v2' v3'
+--     where
+--       c = cube f
+--       v0' = v1 f
+--       v1' = v2 f
+--       v2' = center f
+--       v3' = center c
+
+
+-- tetrahedron2 :: Face -> Tetrahedron
+-- tetrahedron2 f =
+--     Tetrahedron c v0' v1' v2' v3'
+--     where
+--       c = cube f
+--       v0' = v2 f
+--       v1' = v3 f
+--       v2' = center f
+--       v3' = center c
+
+
+-- tetrahedron3 :: Face -> Tetrahedron
+-- tetrahedron3 f =
+--     Tetrahedron c v0' v1' v2' v3'
+--     where
+--       c = cube f
+--       v0' = v3 f
+--       v1' = v0 f
+--       v2' = center f
+--       v3' = center c
+
+-- tetrahedrons :: Cube -> [Tetrahedron]
+-- tetrahedrons c =
+--     concat [
+--  [tetrahedron0 f0, tetrahedron1 f0, tetrahedron2 f0, tetrahedron3 f0],
+--  [tetrahedron0 f1, tetrahedron1 f1, tetrahedron2 f1, tetrahedron3 f2],
+--  [tetrahedron0 f2, tetrahedron1 f2, tetrahedron2 f2, tetrahedron3 f2],
+--  [tetrahedron0 f3, tetrahedron1 f3, tetrahedron2 f3, tetrahedron3 f3],
+--  [tetrahedron0 f4, tetrahedron1 f4, tetrahedron2 f4, tetrahedron3 f4],
+--  [tetrahedron0 f5, tetrahedron1 f5, tetrahedron2 f5, tetrahedron3 f5] ]
+--     where
+--       f0 = face0 c
+--       f1 = face1 c
+--       f2 = face2 c
+--       f3 = face3 c
+--       f4 = face4 c
+--       f5 = face5 c
diff --git a/src/FunctionValues.hs b/src/FunctionValues.hs
index d477a38..e9ff064 100644
--- a/src/FunctionValues.hs
+++ b/src/FunctionValues.hs
@@ -1,6 +1,8 @@
 module FunctionValues
 where
 
+import Prelude hiding (LT)
+
 import Cardinal
 
 data FunctionValues =
@@ -44,7 +46,63 @@ eval f L = left f
 eval f R = right f
 eval f T = top f
 eval f D = down f
+eval f FL = front_left f
+eval f FR = front_right f
+eval f FD = front_down f
+eval f FT = front_top f
+eval f BL = back_left f
+eval f BR = back_right f
+eval f BD = back_down f
+eval f BT = back_top f
+eval f LD = left_down f
+eval f LT = left_top f
+eval f RD = right_down f
+eval f RT = right_top f
+eval f FLD = front_left_down f
+eval f FLT = front_left_top f
+eval f FRD = front_right_down f
+eval f FRT = front_right_top f
+eval f BLD = back_left_down f
+eval f BLT = back_left_top f
+eval f BRD = back_right_down f
+eval f BRT = back_right_top f
+eval f I   = interior f
+eval _ (Scalar x) = x
 eval f (Sum x y) = (eval f x) + (eval f y)
 eval f (Difference x y) = (eval f x) - (eval f y)
 eval f (Product x y) = (eval f x) * (eval f y)
-eval f (ScalarProduct x y) = x * (eval f y)
+eval f (Quotient x y) = (eval f x) / (eval f y)
+
+value_at :: [[[Double]]] -> Int -> Int -> Int -> Double            
+value_at values i j k =
+    ((values !! k) !! j) !! i
+
+make_values :: [[[Double]]] -> Int -> Int -> Int -> FunctionValues
+make_values values i j k =
+    empty_values { front  = value_at values (i-1) j k,
+                   back   = value_at values (i+1) j k,
+                   left   = value_at values i (j-1) k,
+                   right  = value_at values i (j+1) k,
+                   down   = value_at values i j (k-1),
+                   top    = value_at values i j (k+1),
+                   front_left = value_at values (i-1) (j-1) k,
+                   front_right = value_at values (i-1) (j+1) k,
+                   front_down =value_at values (i-1) j (k-1),
+                   front_top = value_at values (i-1) j (k+1),
+                   back_left = value_at values (i+1) (j-1) k,
+                   back_right = value_at values (i+1) (j+1) k,
+                   back_down = value_at values (i+1) j (k-1),
+                   back_top = value_at values (i+1) j (k+1),
+                   left_down = value_at values i (j-1) (k-1),
+                   left_top = value_at values i (j-1) (k+1),
+                   right_top = value_at values i (j+1) (k+1),
+                   right_down = value_at values i (j+1) (k-1),
+                   front_left_down = value_at values (i-1) (j-1) (k-1),
+                   front_left_top = value_at values (i-1) (j-1) (k+1),
+                   front_right_down = value_at values (i-1) (j+1) (k-1),
+                   front_right_top = value_at values (i-1) (j+1) (k+1),
+                   back_left_down = value_at values (i-1) (j-1) (k-1),
+                   back_left_top = value_at values (i+1) (j-1) (k+1),
+                   back_right_down = value_at values (i+1) (j+1) (k-1),
+                   back_right_top = value_at values (i+1) (j+1) (k+1),
+                   interior = value_at values i j k }
diff --git a/src/Grid.hs b/src/Grid.hs
index e67ac76..6d360e2 100644
--- a/src/Grid.hs
+++ b/src/Grid.hs
@@ -1,13 +1,12 @@
 -- | The Grid module just contains the Grid type and two constructors
 --   for it. We hide the main Grid constructor because we don't want
 --   to allow instantiation of a grid with h <= 0.
-module Grid (empty_grid,
-             function_values,
-             Grid,
-             h,
-             make_grid)
+module Grid
 where
 
+import Cube (Cube(Cube))
+import FunctionValues
+
 -- | Our problem is defined on a Grid. The grid size is given by the
 --   positive number h. The function values are the values of the
 --   function at the grid points, which are distance h from one
@@ -28,3 +27,21 @@ make_grid grid_size values
 -- | Creates an empty grid with grid size 1.
 empty_grid :: Grid
 empty_grid = Grid 1 [[[]]]
+
+
+
+-- This is how we do a 'for' loop in Haskell.
+-- No, seriously.
+cubes :: Grid -> [[[Cube]]]
+cubes g
+    | fvs == [[[]]] = [[[]]]
+    | head fvs == [[]] = [[[]]]
+    | otherwise =
+        [[[ Cube (h g) i j k (make_values fvs i j k) | i <- [0..xsize]]
+              | j <- [0..ysize]]
+              | k <- [0..zsize]]
+    where
+      fvs = function_values g
+      zsize = (length fvs) - 1
+      ysize = (length $ head fvs) - 1
+      xsize = (length $ head $ head fvs) - 1
diff --git a/src/Main.hs b/src/Main.hs
index 4467d17..9ad6199 100644
--- a/src/Main.hs
+++ b/src/Main.hs
@@ -1,14 +1,14 @@
 module Main
 where
 
-import Cube
-import Face
-import Grid
-import Misc (flatten)
-import Point
-import RealFunction
-import Tetrahedron
-import ThreeDimensional
+--import Cube
+--import Face
+--import Grid
+--import Misc (flatten)
+--import Point
+--import RealFunction
+--import Tetrahedron
+--import ThreeDimensional
 
 trilinear :: [[[Double]]]
 trilinear = [ [ [ 1, 2, 3 ],
@@ -48,43 +48,44 @@ dummy = [ [ [ 0, 1, 2 ],
             [ 24, 25, 26 ]]]
 
 
-find_point_value :: RealFunction Point
-find_point_value p = poly p
-    where
-      g0 = make_grid 1 trilinear
-      the_cubes = flatten (cubes g0)
-      good_cubes = filter ((flip contains_point) p) the_cubes
-      target_cube = head good_cubes
-      good_tets = filter ((flip contains_point) p) (tetrahedrons target_cube)
-      target_tetrahedron = head good_tets
-      poly = polynomial target_tetrahedron
+--find_point_value :: RealFunction Point
+--find_point_value p = poly p
+--    where
+--      g0 = make_grid 1 trilinear
+--      the_cubes = flatten (cubes g0)
+--      good_cubes = filter ((flip contains_point) p) the_cubes
+--       target_cube = head good_cubes
+--       good_tets = filter ((flip contains_point) p) (tetrahedrons target_cube)
+--       target_tetrahedron = head good_tets
+--       poly = polynomial target_tetrahedron
 
 main :: IO ()
 main = do
-  print $ find_point_value (0,0,0)
-  print $ find_point_value (1,0,0)
-  print $ find_point_value (2,0,0)
-  print $ find_point_value (0,1,0)
-  print $ find_point_value (1,1,0)
-  print $ find_point_value (2,1,0)
-  print $ find_point_value (0,2,0)
-  print $ find_point_value (1,2,0)
-  print $ find_point_value (2,2,0)
-  print $ find_point_value (0,0,1)
-  print $ find_point_value (1,0,1)
-  print $ find_point_value (2,0,1)
-  print $ find_point_value (0,1,1)
-  print $ find_point_value (1,1,1)
-  print $ find_point_value (2,1,1)
-  print $ find_point_value (0,2,1)
-  print $ find_point_value (1,2,1)
-  print $ find_point_value (2,2,1)
-  print $ find_point_value (0,0,2)
-  print $ find_point_value (1,0,2)
-  print $ find_point_value (2,0,2)
-  print $ find_point_value (0,1,2)
-  print $ find_point_value (1,1,2)
-  print $ find_point_value (2,1,2)
-  print $ find_point_value (0,2,2)
-  print $ find_point_value (1,2,2)
-  print $ find_point_value (2,2,2)
+  putStrLn "Hello, World."
+  -- print $ find_point_value (0,0,0)
+  -- print $ find_point_value (1,0,0)
+  -- print $ find_point_value (2,0,0)
+  -- print $ find_point_value (0,1,0)
+  -- print $ find_point_value (1,1,0)
+  -- print $ find_point_value (2,1,0)
+  -- print $ find_point_value (0,2,0)
+  -- print $ find_point_value (1,2,0)
+  -- print $ find_point_value (2,2,0)
+  -- print $ find_point_value (0,0,1)
+  -- print $ find_point_value (1,0,1)
+  -- print $ find_point_value (2,0,1)
+  -- print $ find_point_value (0,1,1)
+  -- print $ find_point_value (1,1,1)
+  -- print $ find_point_value (2,1,1)
+  -- print $ find_point_value (0,2,1)
+  -- print $ find_point_value (1,2,1)
+  -- print $ find_point_value (2,2,1)
+  -- print $ find_point_value (0,0,2)
+  -- print $ find_point_value (1,0,2)
+  -- print $ find_point_value (2,0,2)
+  -- print $ find_point_value (0,1,2)
+  -- print $ find_point_value (1,1,2)
+  -- print $ find_point_value (2,1,2)
+  -- print $ find_point_value (0,2,2)
+  -- print $ find_point_value (1,2,2)
+  -- print $ find_point_value (2,2,2)
diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs
index ad164d2..81fe91c 100644
--- a/src/Tetrahedron.hs
+++ b/src/Tetrahedron.hs
@@ -2,40 +2,29 @@ module Tetrahedron
 where
 
 import Numeric.LinearAlgebra hiding (i, scale)
+import Prelude hiding (LT)
 
-import Cube (back,
-             Cube(datum),
-             down,
-             front,
-             Gridded,
-             left,
-             right,
-             top)
+import Cardinal
+import FunctionValues
 import Misc (factorial)
 import Point
 import RealFunction
 import ThreeDimensional
 
-data Tetrahedron = Tetrahedron { cube :: Cube,
-                                 v0   :: Point,
-                                 v1   :: Point,
-                                 v2   :: Point,
-                                 v3   :: Point }
+data Tetrahedron = Tetrahedron { fv :: FunctionValues,
+                                 v0 :: Point,
+                                 v1 :: Point,
+                                 v2 :: Point,
+                                 v3 :: Point }
                    deriving (Eq)
 
 instance Show Tetrahedron where
-    show t = "Tetrahedron (Cube: " ++ (show (cube t)) ++ ") " ++
-             "(v0: " ++ (show (v0 t)) ++ ") (v1: " ++ (show (v1 t)) ++
-             ") (v2: " ++ (show (v2 t)) ++ ") (v3: " ++ (show (v3 t)) ++
-             ")\n\n"
-
-instance Gridded Tetrahedron where
-    back t = back (cube t)
-    down t = down (cube t)
-    front t = front (cube t)
-    left t = left (cube t)
-    right t = right (cube t)
-    top t = top (cube t)
+    show t = "Tetrahedron:\n" ++
+             "  fv: " ++ (show (fv t)) ++ "\n" ++
+             "  v0: " ++ (show (v0 t)) ++ "\n" ++
+             "  v1: " ++ (show (v1 t)) ++ "\n" ++
+             "  v2: " ++ (show (v2 t)) ++ "\n" ++
+             "  v3: " ++ (show (v3 t)) ++ "\n"
 
 
 instance ThreeDimensional Tetrahedron where
@@ -88,394 +77,121 @@ beta t i j k l
 
 
 c :: Tetrahedron -> Int -> Int -> Int -> Int -> Double
-c x 0 0 3 0 = datum $ (1/8) * (i + f + l + t + lt + fl + ft + flt)
-    where
-      f = front x
-      flt = front (left (top x))
-      fl = front (left x)
-      ft =  front (top x)
-      i = cube x
-      l = left x
-      lt = left (top x)
-      t = top x
-
-
-c x 0 0 0 3 = datum $ (1/8) * (i + f + r + t + rt + fr + ft + frt)
-    where
-      f = front x
-      fr = front (right x)
-      frt = front (right (top x))
-      ft = front (top x)
-      i = cube x
-      r = right x
-      rt = right (top x)
-      t = top x
-
-
-c x 0 0 2 1 = datum $ (5/24)*(i + f + t + ft) + (1/24)*(l + fl + lt + flt)
-    where
-      f = front x
-      flt = front (left (top x))
-      fl = front (left x)
-      ft = front (top x)
-      i = cube x
-      l = left x
-      lt = left (top x)
-      t = top x
-
-
-c x 0 0 1 2 = datum $ (5/24)*(i + f + t + ft) + (1/24)*(r + fr + rt + frt)
-    where
-      f = front x
-      frt = front (right (top x))
-      fr = front (right x)
-      ft = front (top x)
-      i = cube x
-      r = right x
-      rt = right (top x)
-      t = top x
-
-
-c x 0 1 2 0 = datum $
-              (5/24)*(i + f) +
-              (1/8)*(l + t + fl + ft) +
-              (1/24)*(lt + flt)
-    where
-      f = front x
-      flt = front (left (top x))
-      fl = front (left x)
-      ft = front (top x)
-      i = cube x
-      l = left x
-      lt = left (top x)
-      t = top x
-
-
-c x 0 1 0 2 = datum $
-              (5/24)*(i + f) +
-              (1/8)*(r + t + fr + ft) +
-              (1/24)*(rt + frt)
-    where
-      f = front x
-      fr = front (right x)
-      frt = front (right (top x))
-      ft = front (top x)
-      i = cube x
-      r = right x
-      rt = right (top x)
-      t = top x
-
-
-c x 0 1 1 1 = datum $
-          (13/48)*(i + f) +
-          (7/48)*(t + ft) +
-          (1/32)*(l + r + fl + fr) +
-          (1/96)*(lt + rt + flt + frt)
-    where
-      f = front x
-      flt = front (left (top x))
-      fl = front (left x)
-      fr = front (right x)
-      frt = front (right (top x))
-      ft = front (top x)
-      i = cube x
-      l = left x
-      lt = left (top x)
-      r = right x
-      rt = right (top x)
-      t = top x
-
-
-c x 0 2 1 0 = datum $
-          (13/48)*(i + f) +
-          (17/192)*(l + t + fl + ft) +
-          (1/96)*(lt + flt) +
-          (1/64)*(r + d + fr + fd) +
-          (1/192)*(rt + ld + frt + fld)
-    where
-      d = down x
-      f = front x
-      fd = front (down x)
-      fld = front (left (down x))
-      flt = front (left (top x))
-      fl = front (left x)
-      fr = front (right x)
-      frt = front (right (top x))
-      ft = front (top x)
-      i = cube x
-      l = left x
-      ld = left (down x)
-      lt = left (top x)
-      r = right x
-      rt = right (top x)
-      t = top x
-
-
-c x 0 2 0 1 = datum $
-          (13/48)*(i + f) +
-          (17/192)*(r + t + fr + ft) +
-          (1/96)*(rt + frt) +
-          (1/64)*(l + d + fl + fd) +
-          (1/192)*(rd + lt + flt + frd)
-    where
-      d = down x
-      f = front x
-      fd = front (down x)
-      flt = front (left (top x))
-      fl = front (left x)
-      frd = front (right (down x))
-      fr = front (right x)
-      frt = front (right (top x))
-      ft = front (top x)
-      i = cube x
-      l = left x
-      lt = left (top x)
-      r = right x
-      rd = right (down x)
-      rt = right (top x)
-      t = top x
-
-
-c x 0 3 0 0 = datum $
-          (13/48)*(i + f) +
-          (5/96)*(l + r + t + d + fl + fr + ft + fd) +
-          (1/192)*(rt + rd + lt + ld + frt + frd + flt + fld)
-    where
-      d = down x
-      f = front x
-      fd = front (down x)
-      fld = front (left (down x))
-      flt = front (left (top x))
-      fl = front (left x)
-      frd = front (right (down x))
-      fr = front (right x)
-      frt = front (right (top x))
-      ft = front (top x)
-      i = cube x
-      l = left x
-      ld = left (down x)
-      lt = left (top x)
-      r = right x
-      rd = right (down x)
-      rt = right (top x)
-      t = top x
-
-c x 1 0 2 0 = datum $ (1/4)*i + (1/6)*(f + l + t) + (1/12)*(lt + fl + ft)
-    where
-      f = front x
-      fl = front (left x)
-      ft = front (top x)
-      i = cube x
-      l = left x
-      lt = left (top x)
-      t = top x
-
-
-c x 1 0 0 2 = datum $ (1/4)*i + (1/6)*(f + r + t) + (1/12)*(rt + fr + ft)
-    where
-      f = front x
-      fr = front (right x)
-      ft = front (top x)
-      i = cube x
-      r = right x
-      rt = right (top x)
-      t = top x
-
-
-c x 1 0 1 1 = datum $
-          (1/3)*i +
-          (5/24)*(f + t) +
-          (1/12)*ft +
-          (1/24)*(l + r) +
-          (1/48)*(lt + rt + fl + fr)
-    where
-      f = front x
-      fl = front (left x)
-      fr = front (right x)
-      ft = front (top x)
-      i = cube x
-      l = left x
-      lt = left (top x)
-      r = right x
-      rt = right (top x)
-      t = top x
-
-
-c x 1 1 1 0 = datum $
-          (1/3)*i +
-          (5/24)*f +
-          (1/8)*(l + t) +
-          (5/96)*(fl + ft) +
-          (1/48)*(d + r + lt) +
-          (1/96)*(fd + ld + rt + fr)
-    where
-      d = down x
-      f = front x
-      fd = front (down x)
-      fl = front (left x)
-      fr = front (right x)
-      ft = front (top x)
-      i = cube x
-      l = left x
-      ld = left (down x)
-      lt = left (top x)
-      r = right x
-      rt = right (top x)
-      t = top x
-
-
-c x 1 1 0 1 = datum $
-          (1/3)*i +
-          (5/24)*f +
-          (1/8)*(r + t) +
-          (5/96)*(fr + ft) +
-          (1/48)*(d + l + rt) +
-          (1/96)*(fd + lt + rd + fl)
-    where
-      d = down x
-      f = front x
-      fd = front (down x)
-      fl = front (left x)
-      fr = front (right x)
-      ft = front (top x)
-      i = cube x
-      l = left x
-      lt = left (top x)
-      r = right x
-      rd = right (down x)
-      rt = right (top x)
-      t = top x
-
-
-
-c x 1 2 0 0 = datum $
-          (1/3)*i +
-          (5/24)*f +
-          (7/96)*(l + r + t + d) +
-          (1/32)*(fl + fr + ft + fd) +
-          (1/96)*(rt + rd + lt + ld)
-    where
-      d = down x
-      f = front x
-      fd = front (down x)
-      fl = front (left x)
-      fr = front (right x)
-      ft = front (top x)
-      i = cube x
-      l = left x
-      ld = left (down x)
-      lt = left (top x)
-      r = right x
-      rd = right (down x)
-      rt = right (top x)
-      t = top x
-
-
-c x 2 0 1 0 = datum $
-          (3/8)*i +
-          (7/48)*(f + t + l) +
-          (1/48)*(r + d + b + lt + fl + ft) +
-          (1/96)*(rt + bt + fr + fd + ld + bl)
-    where
-      b = back x
-      bl = back (left x)
-      bt = back (top x)
-      d = down x
-      f = front x
-      fd = front (down x)
-      fl = front (left x)
-      fr = front (right x)
-      ft = front (top x)
-      i = cube x
-      l = left x
-      ld = left (down x)
-      lt = left (top x)
-      r = right x
-      rt = right (top x)
-      t = top x
-
-
-c x 2 0 0 1 = datum $
-          (3/8)*i +
-          (7/48)*(f + t + r) +
-          (1/48)*(l + d + b + rt + fr + ft) +
-          (1/96)*(lt + bt + fl + fd + rd + br)
-    where
-      b = back x
-      br = back (right x)
-      bt = back (top x)
-      d = down x
-      f = front x
-      fd = front (down x)
-      fl = front (left x)
-      fr = front (right x)
-      ft = front (top x)
-      i = cube x
-      l = left x
-      lt = left (top x)
-      r = right x
-      rd = right (down x)
-      rt = right (top x)
-      t = top x
-
-
-
-c x 2 1 0 0 = datum $
-          (3/8)*i +
-          (1/12)*(t + r + l + d) +
-          (1/64)*(ft + fr + fl + fd) +
-          (7/48)*f +
-          (1/48)*b +
-          (1/96)*(rt + ld + lt + rd) +
-          (1/192)*(bt + br + bl + bd)
-    where
-      b = back x
-      bd = back (down x)
-      bl = back (left x)
-      br = back (right x)
-      bt = back (top x)
-      d = down x
-      f = front x
-      fd = front (down x)
-      fl = front (left x)
-      fr = front (right x)
-      ft = front (top x)
-      i = cube x
-      l = left x
-      ld = left (down x)
-      lt = left (top x)
-      r = right x
-      rd = right (down x)
-      rt = right (top x)
-      t = top x
-
-
-c x 3 0 0 0 = datum $
-          (3/8)*i +
-          (1/12)*(t + f + l + r + d + b) +
-          (1/96)*(lt + fl + ft + rt + bt + fr) +
-          (1/96)*(fd + ld + bd + br + rd + bl)
-    where
-      b = back x
-      bd = back (down x)
-      bl = back (left x)
-      br = back (right x)
-      bt = back (top x)
-      d = down x
-      f = front x
-      fd = front (down x)
-      fl = front (left x)
-      fr = front (right x)
-      ft = front (top x)
-      i = cube x
-      l = left x
-      ld = left (down x)
-      lt = left (top x)
-      r = right x
-      rd = right (down x)
-      rt = right (top x)
-      t = top x
-
+c t 0 0 3 0 = eval (fv t) $
+              (1/8) * (I + F + L + T + LT + FL + FT + FLT)
+
+c t 0 0 0 3 = eval (fv t) $
+              (1/8) * (I + F + R + T + RT + FR + FT + FRT)
+
+c t 0 0 2 1 = eval (fv t) $
+              (5/24)*(I + F + T + FT) +
+              (1/24)*(L + FL + LT + FLT)
+
+c t 0 0 1 2 = eval (fv t) $
+              (5/24)*(I + F + T + FT) +
+              (1/24)*(R + FR + RT + FRT)
+
+c t 0 1 2 0 = eval (fv t) $
+              (5/24)*(I + F) +
+              (1/8)*(L + T + FL + FT) +
+              (1/24)*(LT + FLT)
+
+c t 0 1 0 2 = eval (fv t) $
+              (5/24)*(I + F) +
+              (1/8)*(R + T + FR + FT) +
+              (1/24)*(RT + FRT)
+
+c t 0 1 1 1 = eval (fv t) $
+          (13/48)*(I + F) +
+          (7/48)*(T + FT) +
+          (1/32)*(L + R + FL + FR) +
+          (1/96)*(LT + RT + FLT + FRT)
+
+c t 0 2 1 0 = eval (fv t) $
+          (13/48)*(I + F) +
+          (17/192)*(L + T + FL + FT) +
+          (1/96)*(LT + FLT) +
+          (1/64)*(R + D + FR + FD) +
+          (1/192)*(RT + LD + FRT + FLD)
+
+c t 0 2 0 1 = eval (fv t) $
+          (13/48)*(I + F) +
+          (17/192)*(R + T + FR + FT) +
+          (1/96)*(RT + FRT) +
+          (1/64)*(L + D + FL + FD) +
+          (1/192)*(RD + LT + FLT + FRD)
+
+c t 0 3 0 0 = eval (fv t) $
+          (13/48)*(I + F) +
+          (5/96)*(L + R + T + D + FL + FR + FT + FD) +
+          (1/192)*(RT + RD + LT + LD + FRT + FRD + FLT + FLD)
+
+c t 1 0 2 0 = eval (fv t) $
+              (1/4)*I +
+              (1/6)*(F + L + T) +
+              (1/12)*(LT + FL + FT)
+
+c t 1 0 0 2 = eval (fv t) $
+              (1/4)*I +
+              (1/6)*(F + R + T) +
+              (1/12)*(RT + FR + FT)
+
+c t 1 0 1 1 = eval (fv t) $
+          (1/3)*I +
+          (5/24)*(F + T) +
+          (1/12)*FT +
+          (1/24)*(L + R) +
+          (1/48)*(LT + RT + FL + FR)
+
+c t 1 1 1 0 = eval (fv t) $
+          (1/3)*I +
+          (5/24)*F +
+          (1/8)*(L + T) +
+          (5/96)*(FL + FT) +
+          (1/48)*(D + R + LT) +
+          (1/96)*(FD + LD + RT + FR)
+
+c t 1 1 0 1 = eval (fv t) $
+          (1/3)*I +
+          (5/24)*F +
+          (1/8)*(R + T) +
+          (5/96)*(FR + FT) +
+          (1/48)*(D + L + RT) +
+          (1/96)*(FD + LT + RD + FL)
+
+c t 1 2 0 0 = eval (fv t) $
+          (1/3)*I +
+          (5/24)*F +
+          (7/96)*(L + R + T + D) +
+          (1/32)*(FL + FR + FT + FD) +
+          (1/96)*(RT + RD + LT + LD)
+
+c t 2 0 1 0 = eval (fv t) $
+          (3/8)*I +
+          (7/48)*(F + T + L) +
+          (1/48)*(R + D + B + LT + FL + FT) +
+          (1/96)*(RT + BT + FR + FD + LD + BL)
+
+c t 2 0 0 1 = eval (fv t) $
+          (3/8)*I +
+          (7/48)*(F + T + R) +
+          (1/48)*(L + D + B + RT + FR + FT) +
+          (1/96)*(LT + BT + FL + FD + RD + BR)
+
+c t 2 1 0 0 = eval (fv t) $
+          (3/8)*I +
+          (1/12)*(T + R + L + D) +
+          (1/64)*(FT + FR + FL + FD) +
+          (7/48)*F +
+          (1/48)*B +
+          (1/96)*(RT + LD + LT + RD) +
+          (1/192)*(BT + BR + BL + BD)
+
+c t 3 0 0 0 = eval (fv t) $
+          (3/8)*I +
+          (1/12)*(T + F + L + R + D + B) +
+          (1/96)*(LT + FL + FT + RT + BT + FR) +
+          (1/96)*(FD + LD + BD + BR + RD + BL)
 
 c _ _ _ _ _ = error "coefficient index out of bounds"