[spline3.git] / Cube.hs

module Cube
where

import Data.List ( (\\) )
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 ThreeDimensional

data Cube = Cube { h :: Double,
                   i :: Int,
                   j :: Int,
                   k :: Int,
                   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_" ++ 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 0


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

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

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

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

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

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

instance ThreeDimensional Cube where
    -- | The center of Cube_ijk coincides with v_ijk at
    --   (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
    center c = (x, y, z)
           where
             delta = h c
             i' = fromIntegral (i c) :: Double
             j' = fromIntegral (j c) :: Double
             k' = fromIntegral (k c) :: Double
             x = delta * i'
             y = delta * j'
             z = delta * k'

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



-- Face stuff.

-- | The top (in the direction of z) face of the cube.
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)
      v1' = (center c) + (delta, delta, delta)
      v2' = (center c) + (-delta, delta, delta)
      v3' = (center c) + (-delta, -delta, delta)



-- | The back (in the direction of x) face of the cube.
back_face :: Cube -> Face.Face
back_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)


-- The bottom face (in the direction of -z) of the cube.
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)
      v1' = (center c) + (-delta, delta, -delta)
      v2' = (center c) + (delta, delta, -delta)
      v3' = (center c) + (delta, -delta, -delta)



-- | The front (in the direction of -x) face of the cube.
front_face :: Cube -> Face.Face
front_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)

-- | The left (in the direction of -y) face of the cube.
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)


-- | The right (in the direction of y) face of the cube.
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)


tetrahedron :: Cube -> Int -> Tetrahedron

tetrahedron c 0 =
    Tetrahedron (Cube.fv c) v0' v1' v2' v3' vol 0
    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

tetrahedron c 1 =
    Tetrahedron fv' v0' v1' v2' v3' vol 1
    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

tetrahedron c 2 =
    Tetrahedron fv' v0' v1' v2' v3' vol 2
    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

tetrahedron c 3 =
    Tetrahedron fv' v0' v1' v2' v3' vol 3
    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

tetrahedron c 4 =
    Tetrahedron fv' v0' v1' v2' v3' vol 4
    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

tetrahedron c 5 =
    Tetrahedron fv' v0' v1' v2' v3' vol 5
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 6 =
    Tetrahedron fv' v0' v1' v2' v3' vol 6
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 7 =
    Tetrahedron fv' v0' v1' v2' v3' vol 7
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 8 =
    Tetrahedron fv' v0' v1' v2' v3' vol 8
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 9 =
    Tetrahedron fv' v0' v1' v2' v3' vol 9
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 10 =
    Tetrahedron fv' v0' v1' v2' v3' vol 10
    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 (tetrahedron c 0)

      vol = tetrahedra_volume c

tetrahedron c 11 =
    Tetrahedron fv' v0' v1' v2' v3' vol 11
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 12 =
    Tetrahedron fv' v0' v1' v2' v3' vol 12
    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 (tetrahedron c 0))
      vol = tetrahedra_volume c

tetrahedron c 13 =
    Tetrahedron fv' v0' v1' v2' v3' vol 13
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 14 =
    Tetrahedron fv' v0' v1' v2' v3' vol 14
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 15 =
    Tetrahedron fv' v0' v1' v2' v3' vol 15
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 16 =
    Tetrahedron fv' v0' v1' v2' v3' vol 16
    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 (tetrahedron c 0))
      vol = tetrahedra_volume c

tetrahedron c 17 =
    Tetrahedron fv' v0' v1' v2' v3' vol 17
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 18 =
    Tetrahedron fv' v0' v1' v2' v3' vol 18
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 19 =
    Tetrahedron fv' v0' v1' v2' v3' vol 19
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 20 =
    Tetrahedron fv' v0' v1' v2' v3' vol 20
    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 (tetrahedron c 0))
      vol = tetrahedra_volume c

tetrahedron c 21 =
    Tetrahedron fv' v0' v1' v2' v3' vol 21
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 22 =
    Tetrahedron fv' v0' v1' v2' v3' vol 22
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

tetrahedron c 23 =
    Tetrahedron fv' v0' v1' v2' v3' vol 23
    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 (tetrahedron c 0)
      vol = tetrahedra_volume c

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


tetrahedra :: Cube -> [Tetrahedron]
tetrahedra c =
    [ tetrahedron c n | n <- [0..23] ]

-- | All completely contained in the front half of the cube.
front_half_tetrahedra :: Cube -> [Tetrahedron]
front_half_tetrahedra c =
  [ tetrahedron c n | n <- [0,1,2,3,6,12,19,21] ]

-- | All tetrahedra completely contained in the top half of the cube.
top_half_tetrahedra :: Cube -> [Tetrahedron]
top_half_tetrahedra c =
  [ tetrahedron c n | n <- [4,5,6,7,0,10,16,20] ]

-- | All tetrahedra completely contained in the back half of the cube.
back_half_tetrahedra :: Cube -> [Tetrahedron]
back_half_tetrahedra c =
  [ tetrahedron c n | n <- [8,9,10,11,4,14,17,23] ]

-- | All tetrahedra completely contained in the down half of the cube.
down_half_tetrahedra :: Cube -> [Tetrahedron]
down_half_tetrahedra c =
  [ tetrahedron c n | n <- [12,13,14,15,2,8,18,22] ]

-- | All tetrahedra completely contained in the right half of the cube.
right_half_tetrahedra :: Cube -> [Tetrahedron]
right_half_tetrahedra c =
  [ tetrahedron c n | n <- [16,17,18,19,1,5,9,13] ]

-- | All tetrahedra completely contained in the left half of the cube.
left_half_tetrahedra :: Cube -> [Tetrahedron]
left_half_tetrahedra c =
  [ tetrahedron c n | n <- [20,21,22,23,3,7,11,15] ]

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