X-Git-Url: https://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FCube.hs;h=ef2d9ce67ab29e52deae2928605951a1c3904185;hb=b2e1c440b9b1bb99ae564d6600230bbd1f7d204c;hp=23cbe4da70383d80c8fbf1f43d1ccc16a0260620;hpb=6fb9ab6b6068870323e996da931fc04c7710e3e4;p=spline3.git

diff --git a/src/Cube.hs b/src/Cube.hs
index 23cbe4d..ef2d9ce 100644
--- a/src/Cube.hs
+++ b/src/Cube.hs
@@ -1,33 +1,59 @@
 module Cube
 where
 
-import Face
+import Data.List ( (\\) )
+import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
+
+import Cardinal
+import qualified Face (Face(Face, v0, v1, v2, v3))
 import FunctionValues
---import Grid
 import Point
+import Tetrahedron hiding (c)
 import ThreeDimensional
 
 data Cube = Cube { h :: Double,
                    i :: Int,
                    j :: Int,
                    k :: Int,
-                   fv :: FunctionValues }
+                   fv :: FunctionValues,
+                   tetrahedra_volume :: Double }
             deriving (Eq)
 
 
+instance Arbitrary Cube where
+    arbitrary = do
+      (Positive h') <- arbitrary :: Gen (Positive Double)
+      i' <- choose (coordmin, coordmax)
+      j' <- choose (coordmin, coordmax)
+      k' <- choose (coordmin, coordmax)
+      fv' <- arbitrary :: Gen FunctionValues
+      (Positive tet_vol) <- arbitrary :: Gen (Positive Double)
+      return (Cube h' i' j' k' fv' tet_vol)
+        where
+          coordmin = -268435456 -- -(2^29 / 2)
+          coordmax = 268435456  -- +(2^29 / 2)
+
+
 instance Show Cube where
     show c =
-        "Cube_" ++ (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k 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)) ++ ")"
-
+        "Cube_" ++ subscript ++ "\n" ++
+        " h: " ++ (show (h c)) ++ "\n" ++
+        " Center: " ++ (show (center c)) ++ "\n" ++
+        " xmin: " ++ (show (xmin c)) ++ "\n" ++
+        " xmax: " ++ (show (xmax c)) ++ "\n" ++
+        " ymin: " ++ (show (ymin c)) ++ "\n" ++
+        " ymax: " ++ (show (ymax c)) ++ "\n" ++
+        " zmin: " ++ (show (zmin c)) ++ "\n" ++
+        " zmax: " ++ (show (zmax c)) ++ "\n" ++
+        " fv: " ++ (show (Cube.fv c)) ++ "\n"
+        where
+          subscript =
+              (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k c))
+
+
+-- | Returns an empty 'Cube'.
 empty_cube :: Cube
-empty_cube = Cube 0 0 0 0 empty_values
+empty_cube = Cube 0 0 0 0 empty_values 0
 
 
 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
@@ -91,130 +117,548 @@ instance ThreeDimensional Cube where
              y = delta * j'
              z = delta * k'
 
-    contains_point c p
-        | (x_coord p) < (xmin c) = False
-        | (x_coord p) > (xmax c) = False
-        | (y_coord p) < (ymin c) = False
-        | (y_coord p) > (ymax c) = False
-        | (z_coord p) < (zmin c) = False
-        | (z_coord p) > (zmax c) = False
+    -- | 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 c (x, y, z)
+        | x < (xmin c) = False
+        | x > (xmax c) = False
+        | y < (ymin c) = False
+        | y > (ymax c) = False
+        | z < (zmin c) = False
+        | z > (zmax c) = False
         | 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)
-
---     (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)
-
---     signum (Cube g1 i1 j1 k1 d1) =
---         Cube g1 (signum i1) (signum j1) (signum k1) (signum d1)
-
---     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)
-
---     recip (Cube g1 i1 j1 k1 d1) =
---         Cube g1 i1 j1 k1 (recip d1)
-
---     fromRational q = empty_cube { datum = fromRational q }
-
-
-
--- | 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'
-
-
-
-
-
 
 -- 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'
+top_face :: Cube -> Face.Face
+top_face c = Face.Face v0' v1' v2' v3'
     where
       delta = (1/2)*(h c)
-      v0' = (center c) + (-delta, delta, delta)
+      v0' = (center c) + (delta, -delta, delta)
       v1' = (center c) + (delta, delta, delta)
-      v2' = (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'
+back_face :: Cube -> Face.Face
+back_face c = Face.Face v0' v1' v2' v3'
     where
       delta = (1/2)*(h c)
-      v0' = (center c) + (delta, delta, delta)
+      v0' = (center c) + (delta, -delta, -delta)
       v1' = (center c) + (delta, delta, -delta)
-      v2' = (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'
+down_face :: Cube -> Face.Face
+down_face c = Face.Face v0' v1' v2' v3'
     where
       delta = (1/2)*(h c)
-      v0' = (center c) + (delta, delta, -delta)
+      v0' = (center c) + (-delta, -delta, -delta)
       v1' = (center c) + (-delta, delta, -delta)
-      v2' = (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'
+front_face :: Cube -> Face.Face
+front_face c = Face.Face v0' v1' v2' v3'
     where
       delta = (1/2)*(h c)
-      v0' = (center c) + (-delta, delta, -delta)
+      v0' = (center c) + (-delta, -delta, delta)
       v1' = (center c) + (-delta, delta, delta)
-      v2' = (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'
+left_face :: Cube -> Face.Face
+left_face c = Face.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)
+      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'
+right_face :: Cube -> Face.Face
+right_face c = Face.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)
+      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 :: Cube -> Tetrahedron
+tetrahedron0 c =
+    Tetrahedron (Cube.fv c) v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (front_face c)
+      v2' = Face.v0 (front_face c)
+      v3' = Face.v1 (front_face c)
+      vol = tetrahedra_volume c
+
+tetrahedron1 :: Cube -> Tetrahedron
+tetrahedron1 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (front_face c)
+      v2' = Face.v1 (front_face c)
+      v3' = Face.v2 (front_face c)
+      fv' = rotate ccwx (Cube.fv c)
+      vol = tetrahedra_volume c
+
+tetrahedron2 :: Cube -> Tetrahedron
+tetrahedron2 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (front_face c)
+      v2' = Face.v2 (front_face c)
+      v3' = Face.v3 (front_face c)
+      fv' = rotate ccwx $ rotate ccwx $ Cube.fv c
+      vol = tetrahedra_volume c
+
+tetrahedron3 :: Cube -> Tetrahedron
+tetrahedron3 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (front_face c)
+      v2' = Face.v3 (front_face c)
+      v3' = Face.v0 (front_face c)
+      fv' = rotate cwx (Cube.fv c)
+      vol = tetrahedra_volume c
+
+tetrahedron4 :: Cube -> Tetrahedron
+tetrahedron4 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (top_face c)
+      v2' = Face.v0 (top_face c)
+      v3' = Face.v1 (top_face c)
+      fv' = rotate cwy (Cube.fv c)
+      vol = tetrahedra_volume c
+
+tetrahedron5 :: Cube -> Tetrahedron
+tetrahedron5 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (top_face c)
+      v2' = Face.v1 (top_face c)
+      v3' = Face.v2 (top_face c)
+      fv' = rotate cwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+tetrahedron6 :: Cube -> Tetrahedron
+tetrahedron6 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (top_face c)
+      v2' = Face.v2 (top_face c)
+      v3' = Face.v3 (top_face c)
+      fv' = rotate cwy $ rotate cwz
+                       $ rotate cwz
+                       $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+tetrahedron7 :: Cube -> Tetrahedron
+tetrahedron7 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (top_face c)
+      v2' = Face.v3 (top_face c)
+      v3' = Face.v0 (top_face c)
+      fv' = rotate cwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+tetrahedron8 :: Cube -> Tetrahedron
+tetrahedron8 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (back_face c)
+      v2' = Face.v0 (back_face c)
+      v3' = Face.v1 (back_face c)
+      fv' = rotate cwy $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+tetrahedron9 :: Cube -> Tetrahedron
+tetrahedron9 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (back_face c)
+      v2' = Face.v1 (back_face c)
+      v3' = Face.v2 (back_face c)
+      fv' = rotate cwy $ rotate cwy
+                       $ rotate cwx
+                       $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+tetrahedron10 :: Cube -> Tetrahedron
+tetrahedron10 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (back_face c)
+      v2' = Face.v2 (back_face c)
+      v3' = Face.v3 (back_face c)
+      fv' = rotate cwy $ rotate cwy
+                       $ rotate cwx
+                       $ rotate cwx
+                       $ Tetrahedron.fv (tetrahedron0 c)
+
+      vol = tetrahedra_volume c
+
+tetrahedron11 :: Cube -> Tetrahedron
+tetrahedron11 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (back_face c)
+      v2' = Face.v3 (back_face c)
+      v3' = Face.v0 (back_face c)
+      fv' = rotate cwy $ rotate cwy
+                       $ rotate ccwx
+                       $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+
+tetrahedron12 :: Cube -> Tetrahedron
+tetrahedron12 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (down_face c)
+      v2' = Face.v0 (down_face c)
+      v3' = Face.v1 (down_face c)
+      fv' = rotate ccwy (Tetrahedron.fv (tetrahedron0 c))
+      vol = tetrahedra_volume c
+
+
+tetrahedron13 :: Cube -> Tetrahedron
+tetrahedron13 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (down_face c)
+      v2' = Face.v1 (down_face c)
+      v3' = Face.v2 (down_face c)
+      fv' = rotate ccwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+
+tetrahedron14 :: Cube -> Tetrahedron
+tetrahedron14 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (down_face c)
+      v2' = Face.v2 (down_face c)
+      v3' = Face.v3 (down_face c)
+      fv' = rotate ccwy $ rotate ccwz
+                        $ rotate ccwz
+                        $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+
+tetrahedron15 :: Cube -> Tetrahedron
+tetrahedron15 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (down_face c)
+      v2' = Face.v3 (down_face c)
+      v3' = Face.v0 (down_face c)
+      fv' = rotate ccwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+
+tetrahedron16 :: Cube -> Tetrahedron
+tetrahedron16 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (right_face c)
+      v2' = Face.v0 (right_face c)
+      v3' = Face.v1 (right_face c)
+      fv' = rotate ccwz (Tetrahedron.fv (tetrahedron0 c))
+      vol = tetrahedra_volume c
+
+
+tetrahedron17 :: Cube -> Tetrahedron
+tetrahedron17 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (right_face c)
+      v2' = Face.v1 (right_face c)
+      v3' = Face.v2 (right_face c)
+      fv' = rotate ccwz $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+
+tetrahedron18 :: Cube -> Tetrahedron
+tetrahedron18 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (right_face c)
+      v2' = Face.v2 (right_face c)
+      v3' = Face.v3 (right_face c)
+      fv' = rotate ccwz $ rotate cwy
+                        $ rotate cwy
+                        $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+
+tetrahedron19 :: Cube -> Tetrahedron
+tetrahedron19 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (right_face c)
+      v2' = Face.v3 (right_face c)
+      v3' = Face.v0 (right_face c)
+      fv' = rotate ccwz $ rotate ccwy
+                        $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+
+tetrahedron20 :: Cube -> Tetrahedron
+tetrahedron20 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (left_face c)
+      v2' = Face.v0 (left_face c)
+      v3' = Face.v1 (left_face c)
+      fv' = rotate cwz (Tetrahedron.fv (tetrahedron0 c))
+      vol = tetrahedra_volume c
 
+
+tetrahedron21 :: Cube -> Tetrahedron
+tetrahedron21 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (left_face c)
+      v2' = Face.v1 (left_face c)
+      v3' = Face.v2 (left_face c)
+      fv' = rotate cwz $ rotate ccwy $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+
+tetrahedron22 :: Cube -> Tetrahedron
+tetrahedron22 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (left_face c)
+      v2' = Face.v2 (left_face c)
+      v3' = Face.v3 (left_face c)
+      fv' = rotate cwz $ rotate ccwy
+                       $ rotate ccwy
+                       $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+
+tetrahedron23 :: Cube -> Tetrahedron
+tetrahedron23 c =
+    Tetrahedron fv' v0' v1' v2' v3' vol
+    where
+      v0' = center c
+      v1' = center (left_face c)
+      v2' = Face.v3 (left_face c)
+      v3' = Face.v0 (left_face c)
+      fv' = rotate cwz $ rotate cwy
+                       $ Tetrahedron.fv (tetrahedron0 c)
+      vol = tetrahedra_volume c
+
+
+tetrahedra :: Cube -> [Tetrahedron]
+tetrahedra c =
+    [tetrahedron0 c,
+     tetrahedron1 c,
+     tetrahedron2 c,
+     tetrahedron3 c,
+     tetrahedron4 c,
+     tetrahedron5 c,
+     tetrahedron6 c,
+     tetrahedron7 c,
+     tetrahedron8 c,
+     tetrahedron9 c,
+     tetrahedron10 c,
+     tetrahedron11 c,
+     tetrahedron12 c,
+     tetrahedron13 c,
+     tetrahedron14 c,
+     tetrahedron15 c,
+     tetrahedron16 c,
+     tetrahedron17 c,
+     tetrahedron18 c,
+     tetrahedron19 c,
+     tetrahedron20 c,
+     tetrahedron21 c,
+     tetrahedron22 c,
+     tetrahedron23 c]
+
+-- | All completely contained in the front half of the cube.
+front_half_tetrahedra :: Cube -> [Tetrahedron]
+front_half_tetrahedra c =
+  [tetrahedron0 c,
+   tetrahedron1 c,
+   tetrahedron2 c,
+   tetrahedron3 c,
+   tetrahedron6 c,
+   tetrahedron12 c,
+   tetrahedron19 c,
+   tetrahedron21 c]
+
+
+-- | All tetrahedra completely contained in the top half of the cube.
+top_half_tetrahedra :: Cube -> [Tetrahedron]
+top_half_tetrahedra c =
+  [tetrahedron4 c,
+   tetrahedron5 c,
+   tetrahedron6 c,
+   tetrahedron7 c,
+   tetrahedron0 c,
+   tetrahedron10 c,
+   tetrahedron16 c,
+   tetrahedron20 c]
+
+
+-- | All tetrahedra completely contained in the back half of the cube.
+back_half_tetrahedra :: Cube -> [Tetrahedron]
+back_half_tetrahedra c =
+  [tetrahedron8 c,
+   tetrahedron9 c,
+   tetrahedron10 c,
+   tetrahedron11 c,
+   tetrahedron4 c,
+   tetrahedron14 c,
+   tetrahedron17 c,
+   tetrahedron23 c]
+
+
+-- | All tetrahedra completely contained in the down half of the cube.
+down_half_tetrahedra :: Cube -> [Tetrahedron]
+down_half_tetrahedra c =
+  [tetrahedron12 c,
+   tetrahedron13 c,
+   tetrahedron14 c,
+   tetrahedron15 c,
+   tetrahedron2 c,
+   tetrahedron8 c,
+   tetrahedron18 c,
+   tetrahedron22 c]
+
+
+-- | All tetrahedra completely contained in the right half of the cube.
+right_half_tetrahedra :: Cube -> [Tetrahedron]
+right_half_tetrahedra c =
+  [tetrahedron16 c,
+   tetrahedron17 c,
+   tetrahedron18 c,
+   tetrahedron19 c,
+   tetrahedron1 c,
+   tetrahedron5 c,
+   tetrahedron9 c,
+   tetrahedron13 c]
+
+
+-- | All tetrahedra completely contained in the left half of the cube.
+left_half_tetrahedra :: Cube -> [Tetrahedron]
+left_half_tetrahedra c =
+  [tetrahedron20 c,
+   tetrahedron21 c,
+   tetrahedron22 c,
+   tetrahedron23 c,
+   tetrahedron3 c,
+   tetrahedron7 c,
+   tetrahedron11 c,
+   tetrahedron15 c]
+
+
+in_top_half :: Cube -> Point -> Bool
+in_top_half c (_,_,z) =
+  distance_from_top <= distance_from_bottom
+  where
+    distance_from_top = abs $ (zmax c) - z
+    distance_from_bottom = abs $ (zmin c) - z
+
+in_front_half :: Cube -> Point -> Bool
+in_front_half c (x,_,_) =
+    distance_from_front <= distance_from_back
+  where
+    distance_from_front = abs $ (xmin c) - x
+    distance_from_back = abs $ (xmax c) - x
+
+
+in_left_half :: Cube -> Point -> Bool
+in_left_half c (_,y,_) =
+    distance_from_left <= distance_from_right
+  where
+    distance_from_left = abs $ (ymin c) - y
+    distance_from_right = abs $ (ymax c) - y
+
+
+-- | Takes a 'Cube', and returns the Tetrahedra belonging to it that
+--   contain the given 'Point'. This should be faster than checking
+--   every tetrahedron individually, since we determine which half
+--   (hemisphere?) of the cube the point lies in three times: once in
+--   each dimension. This allows us to eliminate non-candidates
+--   quickly.
+--
+--   This can throw an exception, but the use of 'head' might
+--   save us some unnecessary computations.
+--
+find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
+find_containing_tetrahedron c p =
+  head containing_tetrahedra
+  where
+    candidates = tetrahedra c
+    non_candidates_x =
+        if (in_front_half c p) then
+          back_half_tetrahedra c
+        else
+          front_half_tetrahedra c
+
+    candidates_x = candidates \\ non_candidates_x
+
+    non_candidates_y =
+      if (in_left_half c p) then
+        right_half_tetrahedra c
+      else
+        left_half_tetrahedra c
+
+    candidates_xy = candidates_x \\ non_candidates_y
+
+    non_candidates_z =
+      if (in_top_half c p) then
+        down_half_tetrahedra c
+      else
+        top_half_tetrahedra c
+
+    candidates_xyz = candidates_xy \\ non_candidates_z
+
+    contains_our_point = flip contains_point p
+    containing_tetrahedra = filter contains_our_point candidates_xyz