Remove the Tetrahedron 'number' field.

[spline3.git] / src / Cube.hs

diff --git a/src/Cube.hs b/src/Cube.hs
index 1e83a3e690d4a6b0f39cdc0b0025526650aa81db..b0b153e782e3fd6ab98b4856d7f0cf589ee60c40 100644 (file)
--- a/src/Cube.hs
+++ b/src/Cube.hs
@@ -1,21 +1,31 @@
 module Cube
 where
 
 module Cube
 where
 

-import Data.List ( (\\) )
+import Data.Maybe (fromJust)
+import qualified Data.Vector as V (
+  Vector,
+  findIndex,
+  map,
+  minimum,
+  singleton,
+  snoc,
+  unsafeIndex
+  )

 import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
 
 import Cardinal
 import qualified Face (Face(Face, v0, v1, v2, v3))
 import FunctionValues
 import Point
 import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
 
 import Cardinal
 import qualified Face (Face(Face, v0, v1, v2, v3))
 import FunctionValues
 import Point

-import Tetrahedron hiding (c)
+import Tetrahedron hiding (c, fv)

 import ThreeDimensional
 
 data Cube = Cube { h :: Double,
                    i :: Int,
                    j :: Int,
                    k :: Int,
 import ThreeDimensional
 
 data Cube = Cube { h :: Double,
                    i :: Int,
                    j :: Int,
                    k :: Int,

-                   fv :: FunctionValues }
+                   fv :: FunctionValues,
+                   tetrahedra_volume :: Double }

             deriving (Eq)
 
 
             deriving (Eq)
 
 


@@ -26,7 +36,8 @@ instance Arbitrary Cube where
       j' <- choose (coordmin, coordmax)
       k' <- choose (coordmin, coordmax)
       fv' <- arbitrary :: Gen FunctionValues
       j' <- choose (coordmin, coordmax)
       k' <- choose (coordmin, coordmax)
       fv' <- arbitrary :: Gen FunctionValues

-      return (Cube h' i' j' k' fv')
+      (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)
         where
           coordmin = -268435456 -- -(2^29 / 2)
           coordmax = 268435456  -- +(2^29 / 2)


@@ -51,7 +62,7 @@ instance Show Cube where
 
 -- | Returns an empty 'Cube'.
 empty_cube :: Cube
 
 -- | 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,
 
 
 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,


@@ -197,18 +208,18 @@ right_face c = Face.Face v0' v1' v2' v3'
       v3' = (center c) + (-delta, delta, -delta)
 
 
       v3' = (center c) + (-delta, delta, -delta)
 
 

-tetrahedron0 :: Cube -> Tetrahedron
-tetrahedron0 c =
+tetrahedron :: Cube -> Int -> Tetrahedron
+
+tetrahedron c 0 =

     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)
     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 = 0
+      vol = tetrahedra_volume c

 
 

-tetrahedron1 :: Cube -> Tetrahedron
-tetrahedron1 c =
+tetrahedron c 1 =

     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c


@@ -216,10 +227,9 @@ tetrahedron1 c =
       v2' = Face.v1 (front_face c)
       v3' = Face.v2 (front_face c)
       fv' = rotate ccwx (Cube.fv c)
       v2' = Face.v1 (front_face c)
       v3' = Face.v2 (front_face c)
       fv' = rotate ccwx (Cube.fv c)

-      vol = 0
+      vol = tetrahedra_volume c

 
 

-tetrahedron2 :: Cube -> Tetrahedron
-tetrahedron2 c =
+tetrahedron c 2 =

     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c


@@ -227,10 +237,9 @@ tetrahedron2 c =
       v2' = Face.v2 (front_face c)
       v3' = Face.v3 (front_face c)
       fv' = rotate ccwx $ rotate ccwx $ Cube.fv c
       v2' = Face.v2 (front_face c)
       v3' = Face.v3 (front_face c)
       fv' = rotate ccwx $ rotate ccwx $ Cube.fv c

-      vol = 0
+      vol = tetrahedra_volume c

 
 

-tetrahedron3 :: Cube -> Tetrahedron
-tetrahedron3 c =
+tetrahedron c 3 =

     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c


@@ -238,10 +247,9 @@ tetrahedron3 c =
       v2' = Face.v3 (front_face c)
       v3' = Face.v0 (front_face c)
       fv' = rotate cwx (Cube.fv c)
       v2' = Face.v3 (front_face c)
       v3' = Face.v0 (front_face c)
       fv' = rotate cwx (Cube.fv c)

-      vol = 0
+      vol = tetrahedra_volume c

 
 

-tetrahedron4 :: Cube -> Tetrahedron
-tetrahedron4 c =
+tetrahedron c 4 =

     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c


@@ -249,21 +257,19 @@ tetrahedron4 c =
       v2' = Face.v0 (top_face c)
       v3' = Face.v1 (top_face c)
       fv' = rotate cwy (Cube.fv c)
       v2' = Face.v0 (top_face c)
       v3' = Face.v1 (top_face c)
       fv' = rotate cwy (Cube.fv c)

-      vol = 0
+      vol = tetrahedra_volume c

 
 

-tetrahedron5 :: Cube -> Tetrahedron
-tetrahedron5 c =
+tetrahedron c 5 =

     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)
     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 = 0
+      fv' = rotate cwy $ rotate cwz $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron6 :: Cube -> Tetrahedron
-tetrahedron6 c =
+tetrahedron c 6 =

     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c


@@ -272,33 +278,30 @@ tetrahedron6 c =
       v3' = Face.v3 (top_face c)
       fv' = rotate cwy $ rotate cwz
                        $ rotate cwz
       v3' = Face.v3 (top_face c)
       fv' = rotate cwy $ rotate cwz
                        $ rotate cwz

-                       $ Tetrahedron.fv (tetrahedron0 c)
-      vol = 0
+                       $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron7 :: Cube -> Tetrahedron
-tetrahedron7 c =
+tetrahedron c 7 =

     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)
     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 = 0
+      fv' = rotate cwy $ rotate ccwz $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron8 :: Cube -> Tetrahedron
-tetrahedron8 c =
+tetrahedron c 8 =

     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)
     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 = 0
+      fv' = rotate cwy $ rotate cwy $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron9 :: Cube -> Tetrahedron
-tetrahedron9 c =
+tetrahedron c 9 =

     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c


@@ -307,11 +310,10 @@ tetrahedron9 c =
       v3' = Face.v2 (back_face c)
       fv' = rotate cwy $ rotate cwy
                        $ rotate cwx
       v3' = Face.v2 (back_face c)
       fv' = rotate cwy $ rotate cwy
                        $ rotate cwx

-                       $ Tetrahedron.fv (tetrahedron0 c)
-      vol = 0
+                       $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron10 :: Cube -> Tetrahedron
-tetrahedron10 c =
+tetrahedron c 10 =

     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c


@@ -321,12 +323,11 @@ tetrahedron10 c =
       fv' = rotate cwy $ rotate cwy
                        $ rotate cwx
                        $ rotate cwx
       fv' = rotate cwy $ rotate cwy
                        $ rotate cwx
                        $ rotate cwx

-                       $ Tetrahedron.fv (tetrahedron0 c)
+                       $ fv c

 
 

-      vol = 0
+      vol = tetrahedra_volume c

 
 

-tetrahedron11 :: Cube -> Tetrahedron
-tetrahedron11 c =
+tetrahedron c 11 =

     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c


@@ -335,36 +336,30 @@ tetrahedron11 c =
       v3' = Face.v0 (back_face c)
       fv' = rotate cwy $ rotate cwy
                        $ rotate ccwx
       v3' = Face.v0 (back_face c)
       fv' = rotate cwy $ rotate cwy
                        $ rotate ccwx

-                       $ Tetrahedron.fv (tetrahedron0 c)
-      vol = 0
-
+                       $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron12 :: Cube -> Tetrahedron
-tetrahedron12 c =
+tetrahedron c 12 =

     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)
     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 = 0
-
+      fv' = rotate ccwy $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron13 :: Cube -> Tetrahedron
-tetrahedron13 c =
+tetrahedron c 13 =

     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)
     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 = 0
+      fv' = rotate ccwy $ rotate ccwz $ fv c
+      vol = tetrahedra_volume c

 
 

-
-tetrahedron14 :: Cube -> Tetrahedron
-tetrahedron14 c =
+tetrahedron c 14 =

     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c


@@ -373,48 +368,40 @@ tetrahedron14 c =
       v3' = Face.v3 (down_face c)
       fv' = rotate ccwy $ rotate ccwz
                         $ rotate ccwz
       v3' = Face.v3 (down_face c)
       fv' = rotate ccwy $ rotate ccwz
                         $ rotate ccwz

-                        $ Tetrahedron.fv (tetrahedron0 c)
-      vol = 0
-
+                        $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron15 :: Cube -> Tetrahedron
-tetrahedron15 c =
+tetrahedron c 15 =

     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)
     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 = 0
-
+      fv' = rotate ccwy $ rotate cwz $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron16 :: Cube -> Tetrahedron
-tetrahedron16 c =
+tetrahedron c 16 =

     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)
     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 = 0
+      fv' = rotate ccwz $ fv c
+      vol = tetrahedra_volume c

 
 

-
-tetrahedron17 :: Cube -> Tetrahedron
-tetrahedron17 c =
+tetrahedron c 17 =

     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)
     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 = 0
-
+      fv' = rotate ccwz $ rotate cwy $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron18 :: Cube -> Tetrahedron
-tetrahedron18 c =
+tetrahedron c 18 =

     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c


@@ -423,12 +410,10 @@ tetrahedron18 c =
       v3' = Face.v3 (right_face c)
       fv' = rotate ccwz $ rotate cwy
                         $ rotate cwy
       v3' = Face.v3 (right_face c)
       fv' = rotate ccwz $ rotate cwy
                         $ rotate cwy

-                        $ Tetrahedron.fv (tetrahedron0 c)
-      vol = 0
-
+                        $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron19 :: Cube -> Tetrahedron
-tetrahedron19 c =
+tetrahedron c 19 =

     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c


@@ -436,36 +421,30 @@ tetrahedron19 c =
       v2' = Face.v3 (right_face c)
       v3' = Face.v0 (right_face c)
       fv' = rotate ccwz $ rotate ccwy
       v2' = Face.v3 (right_face c)
       v3' = Face.v0 (right_face c)
       fv' = rotate ccwz $ rotate ccwy

-                        $ Tetrahedron.fv (tetrahedron0 c)
-      vol = 0
+                        $ fv c
+      vol = tetrahedra_volume c

 
 

-
-tetrahedron20 :: Cube -> Tetrahedron
-tetrahedron20 c =
+tetrahedron c 20 =

     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)
     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 = 0
-
+      fv' = rotate cwz $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron21 :: Cube -> Tetrahedron
-tetrahedron21 c =
+tetrahedron c 21 =

     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)
     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 = 0
+      fv' = rotate cwz $ rotate ccwy $ fv c
+      vol = tetrahedra_volume c

 
 

-
-tetrahedron22 :: Cube -> Tetrahedron
-tetrahedron22 c =
+tetrahedron c 22 =

     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c


@@ -474,12 +453,10 @@ tetrahedron22 c =
       v3' = Face.v3 (left_face c)
       fv' = rotate cwz $ rotate ccwy
                        $ rotate ccwy
       v3' = Face.v3 (left_face c)
       fv' = rotate cwz $ rotate ccwy
                        $ rotate ccwy

-                       $ Tetrahedron.fv (tetrahedron0 c)
-      vol = 0
-
+                       $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron23 :: Cube -> Tetrahedron
-tetrahedron23 c =
+tetrahedron c 23 =

     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c
     Tetrahedron fv' v0' v1' v2' v3' vol
     where
       v0' = center c


@@ -487,114 +464,89 @@ tetrahedron23 c =
       v2' = Face.v3 (left_face c)
       v3' = Face.v0 (left_face c)
       fv' = rotate cwz $ rotate cwy
       v2' = Face.v3 (left_face c)
       v3' = Face.v0 (left_face c)
       fv' = rotate cwz $ rotate cwy

-                       $ Tetrahedron.fv (tetrahedron0 c)
-      vol = 0
+                       $ fv c
+      vol = tetrahedra_volume c

 
 

+-- Feels dirty, but whatever.
+tetrahedron _ _ = error "asked for a nonexistent tetrahedron"

 
 

-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]

 
 

+-- Only used in tests, so we don't need the added speed
+-- of Data.Vector.
+tetrahedra :: Cube -> [Tetrahedron]
+tetrahedra c = [ tetrahedron c n | n <- [0..23] ]
+
+front_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
+front_left_top_tetrahedra  c =
+  V.singleton (tetrahedron c 0) `V.snoc`
+    (tetrahedron c 3) `V.snoc`
+    (tetrahedron c 6) `V.snoc`
+    (tetrahedron c 7) `V.snoc`
+    (tetrahedron c 20) `V.snoc`
+    (tetrahedron c 21)
+
+front_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
+front_left_down_tetrahedra  c =
+  V.singleton (tetrahedron c 0) `V.snoc`
+    (tetrahedron c 2) `V.snoc`
+    (tetrahedron c 3) `V.snoc`
+    (tetrahedron c 12) `V.snoc`
+    (tetrahedron c 15) `V.snoc`
+    (tetrahedron c 21)
+
+front_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
+front_right_top_tetrahedra  c =
+  V.singleton (tetrahedron c 0) `V.snoc`
+    (tetrahedron c 1) `V.snoc`
+    (tetrahedron c 5) `V.snoc`
+    (tetrahedron c 6) `V.snoc`
+    (tetrahedron c 16) `V.snoc`
+    (tetrahedron c 19)
+
+front_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
+front_right_down_tetrahedra  c =
+  V.singleton (tetrahedron c 1) `V.snoc`
+    (tetrahedron c 2) `V.snoc`
+    (tetrahedron c 12) `V.snoc`
+    (tetrahedron c 13) `V.snoc`
+    (tetrahedron c 18) `V.snoc`
+    (tetrahedron c 19)
+
+back_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
+back_left_top_tetrahedra  c =
+  V.singleton (tetrahedron c 0) `V.snoc`
+    (tetrahedron c 3) `V.snoc`
+    (tetrahedron c 6) `V.snoc`
+    (tetrahedron c 7) `V.snoc`
+    (tetrahedron c 20) `V.snoc`
+    (tetrahedron c 21)
+
+back_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
+back_left_down_tetrahedra  c =
+  V.singleton (tetrahedron c 8) `V.snoc`
+    (tetrahedron c 11) `V.snoc`
+    (tetrahedron c 14) `V.snoc`
+    (tetrahedron c 15) `V.snoc`
+    (tetrahedron c 22) `V.snoc`
+    (tetrahedron c 23)
+
+back_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
+back_right_top_tetrahedra  c =
+  V.singleton (tetrahedron c 4) `V.snoc`
+    (tetrahedron c 5) `V.snoc`
+    (tetrahedron c 9) `V.snoc`
+    (tetrahedron c 10) `V.snoc`
+    (tetrahedron c 16) `V.snoc`
+    (tetrahedron c 17)
+
+back_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
+back_right_down_tetrahedra  c =
+  V.singleton (tetrahedron c 8) `V.snoc`
+    (tetrahedron c 9) `V.snoc`
+    (tetrahedron c 13) `V.snoc`
+    (tetrahedron c 14) `V.snoc`
+    (tetrahedron c 17) `V.snoc`
+    (tetrahedron c 18)

 
 in_top_half :: Cube -> Point -> Bool
 in_top_half c (_,_,z) =
 
 in_top_half :: Cube -> Point -> Bool
 in_top_half c (_,_,z) =


@@ -631,32 +583,43 @@ in_left_half c (_,y,_) =
 --
 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
 find_containing_tetrahedron c p =
 --
 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
 find_containing_tetrahedron c p =

-  head containing_tetrahedra
+  candidates `V.unsafeIndex` (fromJust lucky_idx)

   where
   where

-    candidates = tetrahedra c
-    non_candidates_x =
-        if (in_front_half c p) then
-          back_half_tetrahedra c
+    front_half = in_front_half c p
+    top_half = in_top_half c p
+    left_half = in_left_half c p
+
+    candidates =
+      if front_half then
+
+        if left_half then
+          if top_half then
+            front_left_top_tetrahedra c
+          else
+            front_left_down_tetrahedra c

         else
         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
+          if top_half then
+            front_right_top_tetrahedra c
+          else
+            front_right_down_tetrahedra c
+
+      else -- bottom half
+
+        if left_half then
+          if top_half then
+            back_left_top_tetrahedra c
+          else
+            back_left_down_tetrahedra c
+        else
+          if top_half then
+            back_right_top_tetrahedra c
+          else
+            back_right_down_tetrahedra c
+
+    -- Use the dot product instead of 'distance' here to save a
+    -- sqrt(). So, "distances" below really means "distances squared."
+    distances = V.map ((dot p) . center) candidates
+    shortest_distance = V.minimum distances
+    lucky_idx = V.findIndex
+                  (\t -> (center t) `dot` p == shortest_distance)
+                  candidates