[spline3.git] / src / Cube.hs

-- The "tetrahedron" function pattern matches on the integers zero
-- through twenty-three, but doesn't handle the "otherwise" case, for
-- performance reasons.
{-# OPTIONS_GHC -Wno-incomplete-patterns #-}

module Cube (
  Cube(..),
  cube_properties,
  find_containing_tetrahedron,
  tetrahedra,
  tetrahedron )
where

import Data.Maybe ( fromJust )
import qualified Data.Vector as V (
  Vector,
  findIndex,
  map,
  minimum,
  singleton,
  snoc,
  unsafeIndex)
import Prelude hiding ( LT )
import Test.Tasty ( TestTree, testGroup )
import Test.Tasty.QuickCheck (
  Arbitrary( arbitrary ),
  Gen,
  Positive( Positive ),
  choose,
  testProperty )
import Cardinal (
  Cardinal(F, B, L, R, D, T, FL, FR, FD, FT,
           BL, BR, BD, BT, LD, LT, RD, RT, I),
  ccwx,
  ccwy,
  ccwz,
  cwx,
  cwy,
  cwz )
import Comparisons ( (~=), (~~=) )
import qualified Face ( Face(..), center )
import FunctionValues ( FunctionValues, eval, rotate )
import Misc ( all_equal, disjoint )
import Point ( Point( Point ), dot )
import Tetrahedron (
  Tetrahedron(Tetrahedron, function_values, v0, v1, v2, v3),
  barycenter,
  c,
  volume )

data Cube = Cube { i  :: !Int,
                   j  :: !Int,
                   k  :: !Int,
                   fv :: !FunctionValues,
                   tetrahedra_volume :: !Double }
            deriving (Eq)


instance Arbitrary Cube where
    arbitrary = do
      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 i' j' k' fv' tet_vol)
      where
        -- The idea here is that, when cubed in the volume formula,
        -- these numbers don't overflow 64 bits. This number is not
        -- magic in any other sense than that it does not cause test
        -- failures, while 2^23 does.
        coordmax = 4194304 :: Int -- 2^22
        coordmin = -coordmax


instance Show Cube where
    show cube =
        "Cube_" ++ subscript ++ "\n" ++
        " Center: " ++ (show (center cube)) ++ "\n" ++
        " xmin: " ++ (show (xmin cube)) ++ "\n" ++
        " xmax: " ++ (show (xmax cube)) ++ "\n" ++
        " ymin: " ++ (show (ymin cube)) ++ "\n" ++
        " ymax: " ++ (show (ymax cube)) ++ "\n" ++
        " zmin: " ++ (show (zmin cube)) ++ "\n" ++
        " zmax: " ++ (show (zmax cube)) ++ "\n"
        where
          subscript =
              (show (i cube)) ++ "," ++ (show (j cube)) ++ "," ++ (show (k cube))


-- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
--   p. 76.
xmin :: Cube -> Double
xmin cube = (i' - 1/2)
    where
      i' = fromIntegral (i cube) :: Double

-- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
--   p. 76.
xmax :: Cube -> Double
xmax cube = (i' + 1/2)
    where
      i' = fromIntegral (i cube) :: Double

-- | The front boundary of the cube. See Sorokina and Zeilfelder,
--   p. 76.
ymin :: Cube -> Double
ymin cube = (j' - 1/2)
    where
      j' = fromIntegral (j cube) :: Double

-- | The back boundary of the cube. See Sorokina and Zeilfelder,
--   p. 76.
ymax :: Cube -> Double
ymax cube = (j' + 1/2)
    where
      j' = fromIntegral (j cube) :: Double

-- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
--   p. 76.
zmin :: Cube -> Double
zmin cube = (k' - 1/2)
    where
      k' = fromIntegral (k cube) :: Double

-- | The top boundary of the cube. See Sorokina and Zeilfelder,
--   p. 76.
zmax :: Cube -> Double
zmax cube = (k' + 1/2)
    where
      k' = fromIntegral (k cube) :: Double


-- | The center of Cube_ijk coincides with v_ijk at
--   (i, j, k). See Sorokina and Zeilfelder, p. 76.
center :: Cube -> Point
center cube =
  Point x y z
  where
    x = fromIntegral (i cube) :: Double
    y = fromIntegral (j cube) :: Double
    z = fromIntegral (k cube) :: Double


-- Face stuff.

-- | The top (in the direction of z) face of the cube.
top_face :: Cube -> Face.Face
top_face cube = Face.Face v0' v1' v2' v3'
    where
      delta = (1/2) :: Double
      cc  = center cube
      v0' = cc + ( Point delta (-delta) delta )
      v1' = cc + ( Point delta delta delta )
      v2' = cc + ( Point (-delta) delta delta )
      v3' = cc + ( Point (-delta) (-delta) delta )



-- | The back (in the direction of x) face of the cube.
back_face :: Cube -> Face.Face
back_face cube = Face.Face v0' v1' v2' v3'
    where
      delta = (1/2) :: Double
      cc  = center cube
      v0' = cc + ( Point delta (-delta) (-delta) )
      v1' = cc + ( Point delta delta (-delta) )
      v2' = cc + ( Point delta delta delta )
      v3' = cc + ( Point delta (-delta) delta )


-- The bottom face (in the direction of -z) of the cube.
down_face :: Cube -> Face.Face
down_face cube = Face.Face v0' v1' v2' v3'
    where
      delta = (1/2) :: Double
      cc  = center cube
      v0' = cc + ( Point (-delta) (-delta) (-delta) )
      v1' = cc + ( Point (-delta) delta (-delta) )
      v2' = cc + ( Point delta delta (-delta) )
      v3' = cc + ( Point delta (-delta) (-delta) )



-- | The front (in the direction of -x) face of the cube.
front_face :: Cube -> Face.Face
front_face cube = Face.Face v0' v1' v2' v3'
    where
      delta = (1/2) :: Double
      cc  = center cube
      v0' = cc + ( Point (-delta) (-delta) delta )
      v1' = cc + ( Point (-delta) delta delta )
      v2' = cc + ( Point (-delta) delta (-delta) )
      v3' = cc + ( Point (-delta) (-delta) (-delta) )

-- | The left (in the direction of -y) face of the cube.
left_face :: Cube -> Face.Face
left_face cube = Face.Face v0' v1' v2' v3'
    where
      delta = (1/2) :: Double
      cc  = center cube
      v0' = cc + ( Point delta (-delta) delta )
      v1' = cc + ( Point (-delta) (-delta) delta )
      v2' = cc + ( Point (-delta) (-delta) (-delta) )
      v3' = cc + ( Point delta (-delta) (-delta) )


-- | The right (in the direction of y) face of the cube.
right_face :: Cube -> Face.Face
right_face cube = Face.Face v0' v1' v2' v3'
    where
      delta = (1/2) :: Double
      cc  = center cube
      v0' = cc + ( Point (-delta) delta delta)
      v1' = cc + ( Point delta  delta delta )
      v2' = cc + ( Point delta delta (-delta) )
      v3' = cc + ( Point (-delta) delta (-delta) )


tetrahedron :: Cube -> Int -> Tetrahedron

tetrahedron cube 0 =
    Tetrahedron (fv cube) v0' v1' v2' v3' vol
    where
      v0' = center cube
      ff  = front_face cube
      v1' = Face.center ff
      v2' = Face.v0 ff
      v3' = Face.v1 ff
      vol = tetrahedra_volume cube

tetrahedron cube 1 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      ff  = front_face cube
      v1' = Face.center ff
      v2' = Face.v1 ff
      v3' = Face.v2 ff
      fv' = rotate ccwx (fv cube)
      vol = tetrahedra_volume cube

tetrahedron cube 2 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      ff  = front_face cube
      v1' = Face.center ff
      v2' = Face.v2 ff
      v3' = Face.v3 ff
      fv' = rotate ccwx $ rotate ccwx $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 3 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      ff  = front_face cube
      v1' = Face.center ff
      v2' = Face.v3 ff
      v3' = Face.v0 ff
      fv' = rotate cwx (fv cube)
      vol = tetrahedra_volume cube

tetrahedron cube 4 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      tf  = top_face cube
      v1' = Face.center tf
      v2' = Face.v0 tf
      v3' = Face.v1 tf
      fv' = rotate cwy (fv cube)
      vol = tetrahedra_volume cube

tetrahedron cube 5 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      tf  = top_face cube
      v1' = Face.center tf
      v2' = Face.v1 tf
      v3' = Face.v2 tf
      fv' = rotate cwy $ rotate cwz $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 6 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      tf  = top_face cube
      v1' = Face.center tf
      v2' = Face.v2 tf
      v3' = Face.v3 tf
      fv' = rotate cwy $ rotate cwz
                       $ rotate cwz
                       $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 7 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      tf  = top_face cube
      v1' = Face.center tf
      v2' = Face.v3 tf
      v3' = Face.v0 tf
      fv' = rotate cwy $ rotate ccwz $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 8 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      bf  = back_face cube
      v1' = Face.center bf
      v2' = Face.v0 bf
      v3' = Face.v1 bf
      fv' = rotate cwy $ rotate cwy $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 9 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      bf  = back_face cube
      v1' = Face.center bf
      v2' = Face.v1 bf
      v3' = Face.v2 bf
      fv' = rotate cwy $ rotate cwy
                       $ rotate cwx
                       $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 10 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      bf  = back_face cube
      v1' = Face.center bf
      v2' = Face.v2 bf
      v3' = Face.v3 bf
      fv' = rotate cwy $ rotate cwy
                       $ rotate cwx
                       $ rotate cwx
                       $ fv cube

      vol = tetrahedra_volume cube

tetrahedron cube 11 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      bf  = back_face cube
      v1' = Face.center bf
      v2' = Face.v3 bf
      v3' = Face.v0 bf
      fv' = rotate cwy $ rotate cwy
                       $ rotate ccwx
                       $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 12 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      df  = down_face cube
      v1' = Face.center df
      v2' = Face.v0 df
      v3' = Face.v1 df
      fv' = rotate ccwy $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 13 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      df  = down_face cube
      v1' = Face.center df
      v2' = Face.v1 df
      v3' = Face.v2 df
      fv' = rotate ccwy $ rotate ccwz $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 14 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      df  = down_face cube
      v1' = Face.center df
      v2' = Face.v2 df
      v3' = Face.v3 df
      fv' = rotate ccwy $ rotate ccwz
                        $ rotate ccwz
                        $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 15 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      df  = down_face cube
      v1' = Face.center df
      v2' = Face.v3 df
      v3' = Face.v0 df
      fv' = rotate ccwy $ rotate cwz $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 16 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      rf  = right_face cube
      v1' = Face.center rf
      v2' = Face.v0 rf
      v3' = Face.v1 rf
      fv' = rotate ccwz $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 17 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      rf  = right_face cube
      v1' = Face.center rf
      v2' = Face.v1 rf
      v3' = Face.v2 rf
      fv' = rotate ccwz $ rotate cwy $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 18 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      rf  = right_face cube
      v1' = Face.center rf
      v2' = Face.v2 rf
      v3' = Face.v3 rf
      fv' = rotate ccwz $ rotate cwy
                        $ rotate cwy
                        $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 19 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      rf  = right_face cube
      v1' = Face.center rf
      v2' = Face.v3 rf
      v3' = Face.v0 rf
      fv' = rotate ccwz $ rotate ccwy
                        $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 20 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      lf  = left_face cube
      v1' = Face.center lf
      v2' = Face.v0 lf
      v3' = Face.v1 lf
      fv' = rotate cwz $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 21 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      lf  = left_face cube
      v1' = Face.center lf
      v2' = Face.v1 lf
      v3' = Face.v2 lf
      fv' = rotate cwz $ rotate ccwy $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 22 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      lf  = left_face cube
      v1' = Face.center lf
      v2' = Face.v2 lf
      v3' = Face.v3 lf
      fv' = rotate cwz $ rotate ccwy
                       $ rotate ccwy
                       $ fv cube
      vol = tetrahedra_volume cube

tetrahedron cube 23 =
    Tetrahedron fv' v0' v1' v2' v3' vol
    where
      v0' = center cube
      lf  = left_face cube
      v1' = Face.center lf
      v2' = Face.v3 lf
      v3' = Face.v0 lf
      fv' = rotate cwz $ rotate cwy
                       $ fv cube
      vol = tetrahedra_volume cube


-- Only used in tests, so we don't need the added speed
-- of Data.Vector.
tetrahedra :: Cube -> [Tetrahedron]
tetrahedra cube = [ tetrahedron cube n | n <- [0..23] ]

front_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
front_left_top_tetrahedra  cube =
  V.singleton (tetrahedron cube 0) `V.snoc`
    (tetrahedron cube 3) `V.snoc`
    (tetrahedron cube 6) `V.snoc`
    (tetrahedron cube 7) `V.snoc`
    (tetrahedron cube 20) `V.snoc`
    (tetrahedron cube 21)

front_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
front_left_down_tetrahedra  cube =
  V.singleton (tetrahedron cube 0) `V.snoc`
    (tetrahedron cube 2) `V.snoc`
    (tetrahedron cube 3) `V.snoc`
    (tetrahedron cube 12) `V.snoc`
    (tetrahedron cube 15) `V.snoc`
    (tetrahedron cube 21)

front_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
front_right_top_tetrahedra  cube =
  V.singleton (tetrahedron cube 0) `V.snoc`
    (tetrahedron cube 1) `V.snoc`
    (tetrahedron cube 5) `V.snoc`
    (tetrahedron cube 6) `V.snoc`
    (tetrahedron cube 16) `V.snoc`
    (tetrahedron cube 19)

front_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
front_right_down_tetrahedra  cube =
  V.singleton (tetrahedron cube 1) `V.snoc`
    (tetrahedron cube 2) `V.snoc`
    (tetrahedron cube 12) `V.snoc`
    (tetrahedron cube 13) `V.snoc`
    (tetrahedron cube 18) `V.snoc`
    (tetrahedron cube 19)

back_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
back_left_top_tetrahedra  cube =
  V.singleton (tetrahedron cube 0) `V.snoc`
    (tetrahedron cube 3) `V.snoc`
    (tetrahedron cube 6) `V.snoc`
    (tetrahedron cube 7) `V.snoc`
    (tetrahedron cube 20) `V.snoc`
    (tetrahedron cube 21)

back_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
back_left_down_tetrahedra  cube =
  V.singleton (tetrahedron cube 8) `V.snoc`
    (tetrahedron cube 11) `V.snoc`
    (tetrahedron cube 14) `V.snoc`
    (tetrahedron cube 15) `V.snoc`
    (tetrahedron cube 22) `V.snoc`
    (tetrahedron cube 23)

back_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
back_right_top_tetrahedra  cube =
  V.singleton (tetrahedron cube 4) `V.snoc`
    (tetrahedron cube 5) `V.snoc`
    (tetrahedron cube 9) `V.snoc`
    (tetrahedron cube 10) `V.snoc`
    (tetrahedron cube 16) `V.snoc`
    (tetrahedron cube 17)

back_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
back_right_down_tetrahedra  cube =
  V.singleton (tetrahedron cube 8) `V.snoc`
    (tetrahedron cube 9) `V.snoc`
    (tetrahedron cube 13) `V.snoc`
    (tetrahedron cube 14) `V.snoc`
    (tetrahedron cube 17) `V.snoc`
    (tetrahedron cube 18)

in_top_half :: Cube -> Point -> Bool
in_top_half cube (Point _ _ z) =
  distance_from_top <= distance_from_bottom
  where
    distance_from_top = abs $ (zmax cube) - z
    distance_from_bottom = abs $ (zmin cube) - z

in_front_half :: Cube -> Point -> Bool
in_front_half cube (Point x _ _) =
    distance_from_front <= distance_from_back
  where
    distance_from_front = abs $ (xmin cube) - x
    distance_from_back = abs $ (xmax cube) - x


in_left_half :: Cube -> Point -> Bool
in_left_half cube (Point _ y _) =
    distance_from_left <= distance_from_right
  where
    distance_from_left = abs $ (ymin cube) - y
    distance_from_right = abs $ (ymax cube) - 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.
--
{-# INLINE find_containing_tetrahedron #-}
find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
find_containing_tetrahedron cube p =
  candidates `V.unsafeIndex` (fromJust lucky_idx)
  where
    front_half = in_front_half cube p
    top_half = in_top_half cube p
    left_half = in_left_half cube p

    candidates :: V.Vector Tetrahedron
    candidates
      | front_half =
        if left_half then
          if top_half then
            front_left_top_tetrahedra cube
          else
            front_left_down_tetrahedra cube
        else
          if top_half then
            front_right_top_tetrahedra cube
          else
            front_right_down_tetrahedra cube

      | otherwise = -- back half
        if left_half then
          if top_half then
            back_left_top_tetrahedra cube
          else
            back_left_down_tetrahedra cube
        else
          if top_half then
            back_right_top_tetrahedra cube
          else
            back_right_down_tetrahedra cube

    -- Use the dot product instead of Euclidean distance here to save
    -- a sqrt(). So, "distances" below really means "distances
    -- squared."
    distances :: V.Vector Double
    distances = V.map ((dot p) . barycenter) candidates

    shortest_distance :: Double
    shortest_distance = V.minimum distances

    -- Compute the index of the tetrahedron with the center closest to
    -- p. This is a bad algorithm, but don't change it! If you make it
    -- smarter by finding the index of shortest_distance in distances
    -- (this should give the same answer and avoids recomputing the
    -- dot product), the program gets slower. Seriously!
    lucky_idx :: Maybe Int
    lucky_idx = V.findIndex
                  (\t -> (barycenter t) `dot` p == shortest_distance)
                  candidates






-- * Tests

prop_opposite_octant_tetrahedra_disjoint1 :: Cube -> Bool
prop_opposite_octant_tetrahedra_disjoint1 cube =
  disjoint (front_left_top_tetrahedra cube) (front_right_down_tetrahedra cube)

prop_opposite_octant_tetrahedra_disjoint2 :: Cube -> Bool
prop_opposite_octant_tetrahedra_disjoint2 cube =
  disjoint (back_left_top_tetrahedra cube) (back_right_down_tetrahedra cube)

prop_opposite_octant_tetrahedra_disjoint3 :: Cube -> Bool
prop_opposite_octant_tetrahedra_disjoint3 cube =
  disjoint (front_left_top_tetrahedra cube) (back_right_top_tetrahedra cube)

prop_opposite_octant_tetrahedra_disjoint4 :: Cube -> Bool
prop_opposite_octant_tetrahedra_disjoint4 cube =
  disjoint (front_left_down_tetrahedra cube) (back_right_down_tetrahedra cube)

prop_opposite_octant_tetrahedra_disjoint5 :: Cube -> Bool
prop_opposite_octant_tetrahedra_disjoint5 cube =
  disjoint (front_left_top_tetrahedra cube) (back_left_down_tetrahedra cube)

prop_opposite_octant_tetrahedra_disjoint6 :: Cube -> Bool
prop_opposite_octant_tetrahedra_disjoint6 cube =
  disjoint (front_right_top_tetrahedra cube) (back_right_down_tetrahedra cube)


-- | Since the grid size is necessarily positive, all tetrahedra
--   (which comprise cubes of positive volume) must have positive
--   volume as well.
prop_all_volumes_positive :: Cube -> Bool
prop_all_volumes_positive cube =
    all (>= 0) volumes
    where
      ts = tetrahedra cube
      volumes = map volume ts


-- | In fact, since all of the tetrahedra are identical, we should
--   already know their volumes. There's 24 tetrahedra to a cube, so
--   we'd expect the volume of each one to be 1/24.
prop_all_volumes_exact :: Cube -> Bool
prop_all_volumes_exact cube =
    and [volume t ~~= 1/24 | t <- tetrahedra cube]

-- | All tetrahedron should have their v0 located at the center of the
--   cube.
prop_v0_all_equal :: Cube -> Bool
prop_v0_all_equal cube = (v0 t0) == (v0 t1)
    where
      t0 = head (tetrahedra cube) -- Doesn't matter which two we choose.
      t1 = head $ tail (tetrahedra cube)


-- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Note that the
--   third and fourth indices of c-t3 have been switched. This is
--   because we store the triangles oriented such that their volume is
--   positive. If T and T-tilde share \<v0,v1,v2\> and v3,v3-tilde point
--   in opposite directions, one of them has to have negative volume!
prop_c0120_identity1 :: Cube -> Bool
prop_c0120_identity1 cube =
   c t0 0 1 2 0 ~= (c t0 0 0 2 1 + c t3 0 0 1 2) / 2
     where
       t0 = tetrahedron cube 0
       t3 = tetrahedron cube 3


-- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
--   'prop_c0120_identity1' with tetrahedrons 1 and 2.
prop_c0120_identity2 :: Cube -> Bool
prop_c0120_identity2 cube =
   c t1 0 1 2 0 ~= (c t1 0 0 2 1 + c t0 0 0 1 2) / 2
     where
       t0 = tetrahedron cube 0
       t1 = tetrahedron cube 1

-- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
--   'prop_c0120_identity1' with tetrahedrons 1 and 2.
prop_c0120_identity3 :: Cube -> Bool
prop_c0120_identity3 cube =
   c t2 0 1 2 0 ~= (c t2 0 0 2 1 + c t1 0 0 1 2) / 2
     where
       t1 = tetrahedron cube 1
       t2 = tetrahedron cube 2

-- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
--   'prop_c0120_identity1' with tetrahedrons 2 and 3.
prop_c0120_identity4 :: Cube -> Bool
prop_c0120_identity4 cube =
   c t3 0 1 2 0 ~= (c t3 0 0 2 1 + c t2 0 0 1 2) / 2
     where
       t2 = tetrahedron cube 2
       t3 = tetrahedron cube 3


-- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
--   'prop_c0120_identity1' with tetrahedrons 4 and 5.
prop_c0120_identity5 :: Cube -> Bool
prop_c0120_identity5 cube =
    c t5 0 1 2 0 ~= (c t5 0 0 2 1 + c t4 0 0 1 2) / 2
    where
      t4 = tetrahedron cube 4
      t5 = tetrahedron cube 5

-- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
--   'prop_c0120_identity1' with tetrahedrons 5 and 6.
prop_c0120_identity6 :: Cube -> Bool
prop_c0120_identity6 cube =
   c t6 0 1 2 0 ~= (c t6 0 0 2 1 + c t5 0 0 1 2) / 2
     where
       t5 = tetrahedron cube 5
       t6 = tetrahedron cube 6


-- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
--   'prop_c0120_identity1' with tetrahedrons 6 and 7.
prop_c0120_identity7 :: Cube -> Bool
prop_c0120_identity7 cube =
   c t7 0 1 2 0 ~= (c t7 0 0 2 1 + c t6 0 0 1 2) / 2
     where
       t6 = tetrahedron cube 6
       t7 = tetrahedron cube 7


-- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
--   'prop_c0120_identity1'.
prop_c0210_identity1 :: Cube -> Bool
prop_c0210_identity1 cube =
    c t0 0 2 1 0 ~= (c t0 0 1 1 1 + c t3 0 1 1 1) / 2
      where
        t0 = tetrahedron cube 0
        t3 = tetrahedron cube 3


-- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
--   'prop_c0120_identity1'.
prop_c0300_identity1 :: Cube -> Bool
prop_c0300_identity1 cube =
    c t0 0 3 0 0 ~= (c t0 0 2 0 1 + c t3 0 2 1 0) / 2
      where
        t0 = tetrahedron cube 0
        t3 = tetrahedron cube 3


-- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
--   'prop_c0120_identity1'.
prop_c1110_identity :: Cube -> Bool
prop_c1110_identity cube =
    c t0 1 1 1 0 ~= (c t0 1 0 1 1 + c t3 1 0 1 1) / 2
      where
        t0 = tetrahedron cube 0
        t3 = tetrahedron cube 3


-- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
--   'prop_c0120_identity1'.
prop_c1200_identity1 :: Cube -> Bool
prop_c1200_identity1 cube =
    c t0 1 2 0 0 ~= (c t0 1 1 0 1 + c t3 1 1 1 0) / 2
      where
        t0 = tetrahedron cube 0
        t3 = tetrahedron cube 3


-- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
--   'prop_c0120_identity1'.
prop_c2100_identity1 :: Cube -> Bool
prop_c2100_identity1 cube =
    c t0 2 1 0 0 ~= (c t0 2 0 0 1 + c t3 2 0 1 0) / 2
      where
        t0 = tetrahedron cube 0
        t3 = tetrahedron cube 3



-- | Given in Sorokina and Zeilfelder, p. 79, (2.7). Note that the
--   third and fourth indices of c-t3 have been switched. This is
--   because we store the triangles oriented such that their volume is
--   positive. If T and T-tilde share \<v0,v1,v2\> and v3,v3-tilde
--   point in opposite directions, one of them has to have negative
--   volume!
prop_c0102_identity1 :: Cube -> Bool
prop_c0102_identity1 cube =
    c t0 0 1 0 2 ~= (c t0 0 0 1 2 + c t1 0 0 2 1) / 2
      where
        t0 = tetrahedron cube 0
        t1 = tetrahedron cube 1


-- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
--   'prop_c0102_identity1'.
prop_c0201_identity1 :: Cube -> Bool
prop_c0201_identity1 cube =
    c t0 0 2 0 1 ~= (c t0 0 1 1 1 + c t1 0 1 1 1) / 2
      where
        t0 = tetrahedron cube 0
        t1 = tetrahedron cube 1


-- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
--   'prop_c0102_identity1'.
prop_c0300_identity2 :: Cube -> Bool
prop_c0300_identity2 cube =
    c t0 0 3 0 0 ~= (c t0 0 2 1 0 + c t1 0 2 0 1) / 2
      where
        t0 = tetrahedron cube 0
        t1 = tetrahedron cube 1


-- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
--   'prop_c0102_identity1'.
prop_c1101_identity :: Cube -> Bool
prop_c1101_identity cube =
    c t0 1 1 0 1 ~= (c t0 1 0 1 1 + c t1 1 0 1 1) / 2
      where
        t0 = tetrahedron cube 0
        t1 = tetrahedron cube 1


-- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
--   'prop_c0102_identity1'.
prop_c1200_identity2 :: Cube -> Bool
prop_c1200_identity2 cube =
    c t0 1 2 0 0 ~= (c t0 1 1 1 0 + c t1 1 1 0 1) / 2
      where
        t0 = tetrahedron cube 0
        t1 = tetrahedron cube 1


-- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
--   'prop_c0102_identity1'.
prop_c2100_identity2 :: Cube -> Bool
prop_c2100_identity2 cube =
    c t0 2 1 0 0 ~= (c t0 2 0 1 0 + c t1 2 0 0 1) / 2
      where
        t0 = tetrahedron cube 0
        t1 = tetrahedron cube 1


-- | Given in Sorokina and Zeilfelder, p. 79, (2.8). The third and
--   fourth indices of c-t6 have been switched. This is because we
--   store the triangles oriented such that their volume is
--   positive. If T and T-tilde share \<v0,v1,v2\> and v3,v3-tilde
--   point in opposite directions, one of them has to have negative
--   volume!
prop_c3000_identity :: Cube -> Bool
prop_c3000_identity cube =
    c t0 3 0 0 0 ~= c t0 2 1 0 0 + c t6 2 1 0 0
                    - ((c t0 2 0 1 0 + c t0 2 0 0 1)/ 2)
      where
        t0 = tetrahedron cube 0
        t6 = tetrahedron cube 6


-- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
--   'prop_c3000_identity'.
prop_c2010_identity :: Cube -> Bool
prop_c2010_identity cube =
    c t0 2 0 1 0 ~= c t0 1 1 1 0 + c t6 1 1 0 1
                    - ((c t0 1 0 2 0 + c t0 1 0 1 1)/ 2)
      where
        t0 = tetrahedron cube 0
        t6 = tetrahedron cube 6


-- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
--   'prop_c3000_identity'.
prop_c2001_identity :: Cube -> Bool
prop_c2001_identity cube =
    c t0 2 0 0 1 ~= c t0 1 1 0 1 + c t6 1 1 1 0
                    - ((c t0 1 0 0 2 + c t0 1 0 1 1)/ 2)
      where
        t0 = tetrahedron cube 0
        t6 = tetrahedron cube 6


-- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
--   'prop_c3000_identity'.
prop_c1020_identity :: Cube -> Bool
prop_c1020_identity cube =
    c t0 1 0 2 0 ~= c t0 0 1 2 0 + c t6 0 1 0 2
                    - ((c t0 0 0 3 0 + c t0 0 0 2 1)/ 2)
      where
        t0 = tetrahedron cube 0
        t6 = tetrahedron cube 6


-- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
--   'prop_c3000_identity'.
prop_c1002_identity :: Cube -> Bool
prop_c1002_identity cube =
    c t0 1 0 0 2 ~= c t0 0 1 0 2 + c t6 0 1 2 0
                    - ((c t0 0 0 0 3 + c t0 0 0 1 2)/ 2)
      where
        t0 = tetrahedron cube 0
        t6 = tetrahedron cube 6


-- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
--   'prop_c3000_identity'.
prop_c1011_identity :: Cube -> Bool
prop_c1011_identity cube =
    c t0 1 0 1 1 ~= c t0 0 1 1 1 + c t6 0 1 1 1 -
                    ((c t0 0 0 1 2 + c t0 0 0 2 1)/ 2)
      where
        t0 = tetrahedron cube 0
        t6 = tetrahedron cube 6


-- | The function values at the interior should be the same for all
--   tetrahedra.
prop_interior_values_all_identical :: Cube -> Bool
prop_interior_values_all_identical cube =
    all_equal [ eval (function_values tet) I | tet <- tetrahedra cube ]


-- | We know what (c t6 2 1 0 0) should be from Sorokina and Zeilfelder, p. 87.
--   This test checks the rotation works as expected.
prop_c_tilde_2100_rotation_correct :: Cube -> Bool
prop_c_tilde_2100_rotation_correct cube =
    expr1 ~= expr2
    where
      t0 = tetrahedron cube 0
      t6 = tetrahedron cube 6

      -- What gets computed for c2100 of t6.
      expr1 = eval (function_values t6) $
          (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)

      -- What should be computed for c2100 of t6.
      expr2 = eval (function_values t0) $
              (3/8)*I +
              (1/12)*(F + R + L + B) +
              (1/64)*(FT + RT + LT + BT) +
              (7/48)*T +
              (1/48)*D +
              (1/96)*(FR + FL + BR + BL) +
              (1/192)*(FD + RD + LD + BD)


-- | We know what (c t6 2 1 0 0) should be from Sorokina and
--   Zeilfelder, p. 87.  This test checks the actual value based on
--   the FunctionValues of the cube.
--
--   If 'prop_c_tilde_2100_rotation_correct' passes, then this test is
--   even meaningful!
prop_c_tilde_2100_correct :: Cube -> Bool
prop_c_tilde_2100_correct cube =
    c t6 2 1 0 0 ~= expected
    where
      t0 = tetrahedron cube 0
      t6 = tetrahedron cube 6
      fvs = function_values t0
      expected = eval fvs $
                  (3/8)*I +
                  (1/12)*(F + R + L + B) +
                  (1/64)*(FT + RT + LT + BT) +
                  (7/48)*T +
                  (1/48)*D +
                  (1/96)*(FR + FL + BR + BL) +
                  (1/192)*(FD + RD + LD + BD)


-- Tests to check that the correct edges are incidental.
prop_t0_shares_edge_with_t1 :: Cube -> Bool
prop_t0_shares_edge_with_t1 cube =
    (v1 t0) == (v1 t1) && (v3 t0) == (v2 t1)
      where
        t0 = tetrahedron cube 0
        t1 = tetrahedron cube 1

prop_t0_shares_edge_with_t3 :: Cube -> Bool
prop_t0_shares_edge_with_t3 cube =
    (v1 t0) == (v1 t3) && (v2 t0) == (v3 t3)
      where
        t0 = tetrahedron cube 0
        t3 = tetrahedron cube 3

prop_t0_shares_edge_with_t6 :: Cube -> Bool
prop_t0_shares_edge_with_t6 cube =
    (v2 t0) == (v3 t6) && (v3 t0) == (v2 t6)
      where
        t0 = tetrahedron cube 0
        t6 = tetrahedron cube 6

prop_t1_shares_edge_with_t2 :: Cube -> Bool
prop_t1_shares_edge_with_t2 cube =
    (v1 t1) == (v1 t2) && (v3 t1) == (v2 t2)
      where
        t1 = tetrahedron cube 1
        t2 = tetrahedron cube 2

prop_t1_shares_edge_with_t19 :: Cube -> Bool
prop_t1_shares_edge_with_t19 cube =
    (v2 t1) == (v3 t19) && (v3 t1) == (v2 t19)
      where
        t1 = tetrahedron cube 1
        t19 = tetrahedron cube 19

prop_t2_shares_edge_with_t3 :: Cube -> Bool
prop_t2_shares_edge_with_t3 cube =
    (v1 t1) == (v1 t2) && (v3 t1) == (v2 t2)
      where
        t1 = tetrahedron cube 1
        t2 = tetrahedron cube 2

prop_t2_shares_edge_with_t12 :: Cube -> Bool
prop_t2_shares_edge_with_t12 cube =
    (v2 t2) == (v3 t12) && (v3 t2) == (v2 t12)
      where
        t2 = tetrahedron cube 2
        t12 = tetrahedron cube 12

prop_t3_shares_edge_with_t21 :: Cube -> Bool
prop_t3_shares_edge_with_t21 cube =
    (v2 t3) == (v3 t21) && (v3 t3) == (v2 t21)
      where
        t3 = tetrahedron cube 3
        t21 = tetrahedron cube 21

prop_t4_shares_edge_with_t5 :: Cube -> Bool
prop_t4_shares_edge_with_t5 cube =
    (v1 t4) == (v1 t5) && (v3 t4) == (v2 t5)
      where
        t4 = tetrahedron cube 4
        t5 = tetrahedron cube 5

prop_t4_shares_edge_with_t7 :: Cube -> Bool
prop_t4_shares_edge_with_t7 cube =
    (v1 t4) == (v1 t7) && (v2 t4) == (v3 t7)
      where
        t4 = tetrahedron cube 4
        t7 = tetrahedron cube 7

prop_t4_shares_edge_with_t10 :: Cube -> Bool
prop_t4_shares_edge_with_t10 cube =
    (v2 t4) == (v3 t10) && (v3 t4) == (v2 t10)
      where
        t4 = tetrahedron cube 4
        t10 = tetrahedron cube 10

prop_t5_shares_edge_with_t6 :: Cube -> Bool
prop_t5_shares_edge_with_t6 cube =
    (v1 t5) == (v1 t6) && (v3 t5) == (v2 t6)
      where
        t5 = tetrahedron cube 5
        t6 = tetrahedron cube 6

prop_t5_shares_edge_with_t16 :: Cube -> Bool
prop_t5_shares_edge_with_t16 cube =
    (v2 t5) == (v3 t16) && (v3 t5) == (v2 t16)
      where
        t5 = tetrahedron cube 5
        t16 = tetrahedron cube 16

prop_t6_shares_edge_with_t7 :: Cube -> Bool
prop_t6_shares_edge_with_t7 cube =
    (v1 t6) == (v1 t7) && (v3 t6) == (v2 t7)
      where
        t6 = tetrahedron cube 6
        t7 = tetrahedron cube 7

prop_t7_shares_edge_with_t20 :: Cube -> Bool
prop_t7_shares_edge_with_t20 cube =
    (v2 t7) == (v3 t20) && (v2 t7) == (v3 t20)
      where
        t7 = tetrahedron cube 7
        t20 = tetrahedron cube 20


p79_26_properties :: TestTree
p79_26_properties =
    testGroup "p. 79, Section (2.6) properties" [
      testProperty "c0120 identity1" prop_c0120_identity1,
      testProperty "c0120 identity2" prop_c0120_identity2,
      testProperty "c0120 identity3" prop_c0120_identity3,
      testProperty "c0120 identity4" prop_c0120_identity4,
      testProperty "c0120 identity5" prop_c0120_identity5,
      testProperty "c0120 identity6" prop_c0120_identity6,
      testProperty "c0120 identity7" prop_c0120_identity7,
      testProperty "c0210 identity1" prop_c0210_identity1,
      testProperty "c0300 identity1" prop_c0300_identity1,
      testProperty "c1110 identity" prop_c1110_identity,
      testProperty "c1200 identity1" prop_c1200_identity1,
      testProperty "c2100 identity1" prop_c2100_identity1]

p79_27_properties :: TestTree
p79_27_properties =
    testGroup "p. 79, Section (2.7) properties" [
      testProperty "c0102 identity1" prop_c0102_identity1,
      testProperty "c0201 identity1" prop_c0201_identity1,
      testProperty "c0300 identity2" prop_c0300_identity2,
      testProperty "c1101 identity" prop_c1101_identity,
      testProperty "c1200 identity2" prop_c1200_identity2,
      testProperty "c2100 identity2" prop_c2100_identity2 ]


p79_28_properties :: TestTree
p79_28_properties =
  testGroup "p. 79, Section (2.8) properties" [
    testProperty "c3000 identity" prop_c3000_identity,
    testProperty "c2010 identity" prop_c2010_identity,
    testProperty "c2001 identity" prop_c2001_identity,
    testProperty "c1020 identity" prop_c1020_identity,
    testProperty "c1002 identity" prop_c1002_identity,
    testProperty "c1011 identity" prop_c1011_identity ]


edge_incidence_tests :: TestTree
edge_incidence_tests =
    testGroup "Edge incidence tests" [
      testProperty "t0 shares edge with t6" prop_t0_shares_edge_with_t6,
      testProperty "t0 shares edge with t1" prop_t0_shares_edge_with_t1,
      testProperty "t0 shares edge with t3" prop_t0_shares_edge_with_t3,
      testProperty "t1 shares edge with t2" prop_t1_shares_edge_with_t2,
      testProperty "t1 shares edge with t19" prop_t1_shares_edge_with_t19,
      testProperty "t2 shares edge with t3" prop_t2_shares_edge_with_t3,
      testProperty "t2 shares edge with t12" prop_t2_shares_edge_with_t12,
      testProperty "t3 shares edge with t21" prop_t3_shares_edge_with_t21,
      testProperty "t4 shares edge with t5" prop_t4_shares_edge_with_t5,
      testProperty "t4 shares edge with t7" prop_t4_shares_edge_with_t7,
      testProperty "t4 shares edge with t10" prop_t4_shares_edge_with_t10,
      testProperty "t5 shares edge with t6" prop_t5_shares_edge_with_t6,
      testProperty "t5 shares edge with t16" prop_t5_shares_edge_with_t16,
      testProperty "t6 shares edge with t7" prop_t6_shares_edge_with_t7,
      testProperty "t7 shares edge with t20" prop_t7_shares_edge_with_t20 ]

cube_properties :: TestTree
cube_properties =
  testGroup "Cube properties" [
    p79_26_properties,
    p79_27_properties,
    p79_28_properties,
    edge_incidence_tests,
    testProperty "opposite octant tetrahedra are disjoint (1)"
      prop_opposite_octant_tetrahedra_disjoint1,
    testProperty "opposite octant tetrahedra are disjoint (2)"
      prop_opposite_octant_tetrahedra_disjoint2,
    testProperty "opposite octant tetrahedra are disjoint (3)"
     prop_opposite_octant_tetrahedra_disjoint3,
    testProperty "opposite octant tetrahedra are disjoint (4)"
      prop_opposite_octant_tetrahedra_disjoint4,
    testProperty "opposite octant tetrahedra are disjoint (5)"
      prop_opposite_octant_tetrahedra_disjoint5,
    testProperty "opposite octant tetrahedra are disjoint (6)"
      prop_opposite_octant_tetrahedra_disjoint6,
    testProperty "all volumes positive" prop_all_volumes_positive,
    testProperty "all volumes exact" prop_all_volumes_exact,
    testProperty "v0 all equal" prop_v0_all_equal,
    testProperty "interior values all identical"
      prop_interior_values_all_identical,
    testProperty "c-tilde_2100 rotation correct"
      prop_c_tilde_2100_rotation_correct,
    testProperty "c-tilde_2100 correct"
      prop_c_tilde_2100_correct ]