Move the Tetrahedron tests into the Tetrehedron module.

[spline3.git] / src / Cube.hs

diff --git a/src/Cube.hs b/src/Cube.hs
index 0122aea298157326439dafebbabb0109ed758476..3c202a7f08b23ab188a2e62166a553bf01a96d8e 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 (Tetrahedron(Tetrahedron))

 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,68 +208,69 @@ 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.fv c) v0' v1' v2' v3'
+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)
     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'
+tetrahedron c 1 =
+    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)
     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'
+tetrahedron c 2 =
+    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
     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'
+tetrahedron c 3 =
+    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)
     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'
+tetrahedron c 4 =
+    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)
     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'
+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)
     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)
+      fv' = rotate cwy $ rotate cwz $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron6 :: Cube -> Tetrahedron
-tetrahedron6 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+tetrahedron c 6 =
+    Tetrahedron fv' v0' v1' v2' v3' vol

     where
       v0' = center c
       v1' = center (top_face c)
     where
       v0' = center c
       v1' = center (top_face c)


@@ -266,31 +278,31 @@ 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)
+                       $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron7 :: Cube -> Tetrahedron
-tetrahedron7 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+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)
     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)
+      fv' = rotate cwy $ rotate ccwz $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron8 :: Cube -> Tetrahedron
-tetrahedron8 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+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)
     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)
+      fv' = rotate cwy $ rotate cwy $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron9 :: Cube -> Tetrahedron
-tetrahedron9 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+tetrahedron c 9 =
+    Tetrahedron fv' v0' v1' v2' v3' vol

     where
       v0' = center c
       v1' = center (back_face c)
     where
       v0' = center c
       v1' = center (back_face c)


@@ -298,11 +310,11 @@ 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)
+                       $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron10 :: Cube -> Tetrahedron
-tetrahedron10 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+tetrahedron c 10 =
+    Tetrahedron fv' v0' v1' v2' v3' vol

     where
       v0' = center c
       v1' = center (back_face c)
     where
       v0' = center c
       v1' = center (back_face c)


@@ -311,12 +323,12 @@ 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 = tetrahedra_volume c

 
 

-tetrahedron11 :: Cube -> Tetrahedron
-tetrahedron11 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+tetrahedron c 11 =
+    Tetrahedron fv' v0' v1' v2' v3' vol

     where
       v0' = center c
       v1' = center (back_face c)
     where
       v0' = center c
       v1' = center (back_face c)


@@ -324,34 +336,31 @@ 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)
-
+                       $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron12 :: Cube -> Tetrahedron
-tetrahedron12 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+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)
     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))
-
+      fv' = rotate ccwy $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron13 :: Cube -> Tetrahedron
-tetrahedron13 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+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)
     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)
+      fv' = rotate ccwy $ rotate ccwz $ fv c
+      vol = tetrahedra_volume c

 
 

-
-tetrahedron14 :: Cube -> Tetrahedron
-tetrahedron14 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+tetrahedron c 14 =
+    Tetrahedron fv' v0' v1' v2' v3' vol

     where
       v0' = center c
       v1' = center (down_face c)
     where
       v0' = center c
       v1' = center (down_face c)


@@ -359,45 +368,41 @@ 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)
-
+                        $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron15 :: Cube -> Tetrahedron
-tetrahedron15 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+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)
     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)
-
+      fv' = rotate ccwy $ rotate cwz $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron16 :: Cube -> Tetrahedron
-tetrahedron16 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+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)
     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))
+      fv' = rotate ccwz $ fv c
+      vol = tetrahedra_volume c

 
 

-
-tetrahedron17 :: Cube -> Tetrahedron
-tetrahedron17 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+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)
     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)
-
+      fv' = rotate ccwz $ rotate cwy $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron18 :: Cube -> Tetrahedron
-tetrahedron18 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+tetrahedron c 18 =
+    Tetrahedron fv' v0' v1' v2' v3' vol

     where
       v0' = center c
       v1' = center (right_face c)
     where
       v0' = center c
       v1' = center (right_face c)


@@ -405,46 +410,42 @@ 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)
-
+                        $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron19 :: Cube -> Tetrahedron
-tetrahedron19 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+tetrahedron c 19 =
+    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
     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)
+                        $ fv c
+      vol = tetrahedra_volume c

 
 

-
-tetrahedron20 :: Cube -> Tetrahedron
-tetrahedron20 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+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)
     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))
-
+      fv' = rotate cwz $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron21 :: Cube -> Tetrahedron
-tetrahedron21 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+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)
     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)
+      fv' = rotate cwz $ rotate ccwy $ fv c
+      vol = tetrahedra_volume c

 
 

-
-tetrahedron22 :: Cube -> Tetrahedron
-tetrahedron22 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+tetrahedron c 22 =
+    Tetrahedron fv' v0' v1' v2' v3' vol

     where
       v0' = center c
       v1' = center (left_face c)
     where
       v0' = center c
       v1' = center (left_face c)


@@ -452,125 +453,100 @@ 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)
-
+                       $ fv c
+      vol = tetrahedra_volume c

 
 

-tetrahedron23 :: Cube -> Tetrahedron
-tetrahedron23 c =
-    Tetrahedron fv' v0' v1' v2' v3'
+tetrahedron c 23 =
+    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
     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)
+                       $ 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) =


@@ -607,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