module Tetrahedron
where

import Numeric.LinearAlgebra hiding (i, scale)

import Cube (back,
             Cube(datum),
             down,
             front,
             Gridded,
             left,
             right,
             top)
import Misc (factorial)
import Point
import RealFunction
import ThreeDimensional

data Tetrahedron = Tetrahedron { cube :: Cube,
                                 v0   :: Point,
                                 v1   :: Point,
                                 v2   :: Point,
                                 v3   :: Point }
                   deriving (Eq)

instance Show Tetrahedron where
    show t = "Tetrahedron (Cube: " ++ (show (cube t)) ++ ") " ++
             "(v0: " ++ (show (v0 t)) ++ ") (v1: " ++ (show (v1 t)) ++
             ") (v2: " ++ (show (v2 t)) ++ ") (v3: " ++ (show (v3 t)) ++
             ")\n\n"

instance Gridded Tetrahedron where
    back t = back (cube t)
    down t = down (cube t)
    front t = front (cube t)
    left t = left (cube t)
    right t = right (cube t)
    top t = top (cube t)


instance ThreeDimensional Tetrahedron where
    center t = ((v0 t) + (v1 t) + (v2 t) + (v3 t)) `scale` (1/4)
    contains_point t p =
        (b0 t p) >= 0 && (b1 t p) >= 0 && (b2 t p) >= 0 && (b3 t p) >= 0


polynomial :: Tetrahedron -> (RealFunction Point)
polynomial t =
    sum [ (c t i j k l) `cmult` (beta t i j k l) | i <- [0..3],
                                                   j <- [0..3],
                                                   k <- [0..3],
                                                   l <- [0..3],
                                                   i + j + k + l == 3]


-- | Returns the domain point of t with indices i,j,k,l.
xi :: Tetrahedron -> Int -> Int -> Int -> Int -> Point
xi t i j k l
   | i + j + k + l == 3 = weighted_sum `scale` (1/3)
   | otherwise = error "xi index out of bounds"
   where
     v0' = (v0 t) `scale` (fromIntegral i)
     v1' = (v1 t) `scale` (fromIntegral j)
     v2' = (v2 t) `scale` (fromIntegral k)
     v3' = (v3 t) `scale` (fromIntegral l)
     weighted_sum = v0' + v1' + v2' + v3'



-- | The Bernstein polynomial on t with indices i,j,k,l. Denoted by a
--   capital 'B' in the Sorokina/Zeilfelder paper.
beta :: Tetrahedron -> Int -> Int -> Int -> Int -> (RealFunction Point)
beta t i j k l
  | (i + j + k + l == 3) =
      coefficient `cmult` (b0_term * b1_term * b2_term * b3_term)
  | otherwise = error "basis function index out of bounds"
  where
    denominator = (factorial i)*(factorial j)*(factorial k)*(factorial l)
    coefficient = 6 / (fromIntegral denominator)
    b0_term = (b0 t) `fexp` i
    b1_term = (b1 t) `fexp` j
    b2_term = (b2 t) `fexp` k
    b3_term = (b3 t) `fexp` l


c :: Tetrahedron -> Int -> Int -> Int -> Int -> Double
c x 0 0 3 0 = datum $ (1/8) * (i + f + l + t + lt + fl + ft + flt)
    where
      f = front x
      flt = front (left (top x))
      fl = front (left x)
      ft =  front (top x)
      i = cube x
      l = left x
      lt = left (top x)
      t = top x


c x 0 0 0 3 = datum $ (1/8) * (i + f + r + t + rt + fr + ft + frt)
    where
      f = front x
      fr = front (right x)
      frt = front (right (top x))
      ft = front (top x)
      i = cube x
      r = right x
      rt = right (top x)
      t = top x


c x 0 0 2 1 = datum $ (5/24)*(i + f + t + ft) + (1/24)*(l + fl + lt + flt)
    where
      f = front x
      flt = front (left (top x))
      fl = front (left x)
      ft = front (top x)
      i = cube x
      l = left x
      lt = left (top x)
      t = top x


c x 0 0 1 2 = datum $ (5/24)*(i + f + t + ft) + (1/24)*(r + fr + rt + frt)
    where
      f = front x
      frt = front (right (top x))
      fr = front (right x)
      ft = front (top x)
      i = cube x
      r = right x
      rt = right (top x)
      t = top x


c x 0 1 2 0 = datum $
              (5/24)*(i + f) +
              (1/8)*(l + t + fl + ft) +
              (1/24)*(lt + flt)
    where
      f = front x
      flt = front (left (top x))
      fl = front (left x)
      ft = front (top x)
      i = cube x
      l = left x
      lt = left (top x)
      t = top x


c x 0 1 0 2 = datum $
              (5/24)*(i + f) +
              (1/8)*(r + t + fr + ft) +
              (1/24)*(rt + frt)
    where
      f = front x
      fr = front (right x)
      frt = front (right (top x))
      ft = front (top x)
      i = cube x
      r = right x
      rt = right (top x)
      t = top x


c x 0 1 1 1 = datum $
          (13/48)*(i + f) +
          (7/48)*(t + ft) +
          (1/32)*(l + r + fl + fr) +
          (1/96)*(lt + rt + flt + frt)
    where
      f = front x
      flt = front (left (top x))
      fl = front (left x)
      fr = front (right x)
      frt = front (right (top x))
      ft = front (top x)
      i = cube x
      l = left x
      lt = left (top x)
      r = right x
      rt = right (top x)
      t = top x


c x 0 2 1 0 = datum $
          (13/48)*(i + f) +
          (17/192)*(l + t + fl + ft) +
          (1/96)*(lt + flt) +
          (1/64)*(r + d + fr + fd) +
          (1/192)*(rt + ld + frt + fld)
    where
      d = down x
      f = front x
      fd = front (down x)
      fld = front (left (down x))
      flt = front (left (top x))
      fl = front (left x)
      fr = front (right x)
      frt = front (right (top x))
      ft = front (top x)
      i = cube x
      l = left x
      ld = left (down x)
      lt = left (top x)
      r = right x
      rt = right (top x)
      t = top x


c x 0 2 0 1 = datum $
          (13/48)*(i + f) +
          (17/192)*(r + t + fr + ft) +
          (1/96)*(rt + frt) +
          (1/64)*(l + d + fl + fd) +
          (1/192)*(rd + lt + flt + frd)
    where
      d = down x
      f = front x
      fd = front (down x)
      flt = front (left (top x))
      fl = front (left x)
      frd = front (right (down x))
      fr = front (right x)
      frt = front (right (top x))
      ft = front (top x)
      i = cube x
      l = left x
      lt = left (top x)
      r = right x
      rd = right (down x)
      rt = right (top x)
      t = top x


c x 0 3 0 0 = datum $
          (13/48)*(i + f) +
          (5/96)*(l + r + t + d + fl + fr + ft + fd) +
          (1/192)*(rt + rd + lt + ld + frt + frd + flt + fld)
    where
      d = down x
      f = front x
      fd = front (down x)
      fld = front (left (down x))
      flt = front (left (top x))
      fl = front (left x)
      frd = front (right (down x))
      fr = front (right x)
      frt = front (right (top x))
      ft = front (top x)
      i = cube x
      l = left x
      ld = left (down x)
      lt = left (top x)
      r = right x
      rd = right (down x)
      rt = right (top x)
      t = top x

c x 1 0 2 0 = datum $ (1/4)*i + (1/6)*(f + l + t) + (1/12)*(lt + fl + ft)
    where
      f = front x
      fl = front (left x)
      ft = front (top x)
      i = cube x
      l = left x
      lt = left (top x)
      t = top x


c x 1 0 0 2 = datum $ (1/4)*i + (1/6)*(f + r + t) + (1/12)*(rt + fr + ft)
    where
      f = front x
      fr = front (right x)
      ft = front (top x)
      i = cube x
      r = right x
      rt = right (top x)
      t = top x


c x 1 0 1 1 = datum $
          (1/3)*i +
          (5/24)*(f + t) +
          (1/12)*ft +
          (1/24)*(l + r) +
          (1/48)*(lt + rt + fl + fr)
    where
      f = front x
      fl = front (left x)
      fr = front (right x)
      ft = front (top x)
      i = cube x
      l = left x
      lt = left (top x)
      r = right x
      rt = right (top x)
      t = top x


c x 1 1 1 0 = datum $
          (1/3)*i +
          (5/24)*f +
          (1/8)*(l + t) +
          (5/96)*(fl + ft) +
          (1/48)*(d + r + lt) +
          (1/96)*(fd + ld + rt + fr)
    where
      d = down x
      f = front x
      fd = front (down x)
      fl = front (left x)
      fr = front (right x)
      ft = front (top x)
      i = cube x
      l = left x
      ld = left (down x)
      lt = left (top x)
      r = right x
      rt = right (top x)
      t = top x


c x 1 1 0 1 = datum $
          (1/3)*i +
          (5/24)*f +
          (1/8)*(r + t) +
          (5/96)*(fr + ft) +
          (1/48)*(d + l + rt) +
          (1/96)*(fd + lt + rd + fl)
    where
      d = down x
      f = front x
      fd = front (down x)
      fl = front (left x)
      fr = front (right x)
      ft = front (top x)
      i = cube x
      l = left x
      lt = left (top x)
      r = right x
      rd = right (down x)
      rt = right (top x)
      t = top x



c x 1 2 0 0 = datum $
          (1/3)*i +
          (5/24)*f +
          (7/96)*(l + r + t + d) +
          (1/32)*(fl + fr + ft + fd) +
          (1/96)*(rt + rd + lt + ld)
    where
      d = down x
      f = front x
      fd = front (down x)
      fl = front (left x)
      fr = front (right x)
      ft = front (top x)
      i = cube x
      l = left x
      ld = left (down x)
      lt = left (top x)
      r = right x
      rd = right (down x)
      rt = right (top x)
      t = top x


c x 2 0 1 0 = datum $
          (3/8)*i +
          (7/48)*(f + t + l) +
          (1/48)*(r + d + b + lt + fl + ft) +
          (1/96)*(rt + bt + fr + fd + ld + bl)
    where
      b = back x
      bl = back (left x)
      bt = back (top x)
      d = down x
      f = front x
      fd = front (down x)
      fl = front (left x)
      fr = front (right x)
      ft = front (top x)
      i = cube x
      l = left x
      ld = left (down x)
      lt = left (top x)
      r = right x
      rt = right (top x)
      t = top x


c x 2 0 0 1 = datum $
          (3/8)*i +
          (7/48)*(f + t + r) +
          (1/48)*(l + d + b + rt + fr + ft) +
          (1/96)*(lt + bt + fl + fd + rd + br)
    where
      b = back x
      br = back (right x)
      bt = back (top x)
      d = down x
      f = front x
      fd = front (down x)
      fl = front (left x)
      fr = front (right x)
      ft = front (top x)
      i = cube x
      l = left x
      lt = left (top x)
      r = right x
      rd = right (down x)
      rt = right (top x)
      t = top x



c x 2 1 0 0 = datum $
          (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)
    where
      b = back x
      bd = back (down x)
      bl = back (left x)
      br = back (right x)
      bt = back (top x)
      d = down x
      f = front x
      fd = front (down x)
      fl = front (left x)
      fr = front (right x)
      ft = front (top x)
      i = cube x
      l = left x
      ld = left (down x)
      lt = left (top x)
      r = right x
      rd = right (down x)
      rt = right (top x)
      t = top x


c x 3 0 0 0 = datum $
          (3/8)*i +
          (1/12)*(t + f + l + r + d + b) +
          (1/96)*(lt + fl + ft + rt + bt + fr) +
          (1/96)*(fd + ld + bd + br + rd + bl)
    where
      b = back x
      bd = back (down x)
      bl = back (left x)
      br = back (right x)
      bt = back (top x)
      d = down x
      f = front x
      fd = front (down x)
      fl = front (left x)
      fr = front (right x)
      ft = front (top x)
      i = cube x
      l = left x
      ld = left (down x)
      lt = left (top x)
      r = right x
      rd = right (down x)
      rt = right (top x)
      t = top x


c _ _ _ _ _ = error "coefficient index out of bounds"



vol_matrix :: Tetrahedron -> Matrix Double
vol_matrix t = (4><4) $
               [1,  1,  1,  1,
                x1, x2, x3, x4,
                y1, y2, y3, y4,
                z1, z2, z3, z4 ]
    where
      x1 = x_coord (v0 t)
      x2 = x_coord (v1 t)
      x3 = x_coord (v2 t)
      x4 = x_coord (v3 t)
      y1 = y_coord (v0 t)
      y2 = y_coord (v1 t)
      y3 = y_coord (v2 t)
      y4 = y_coord (v3 t)
      z1 = z_coord (v0 t)
      z2 = z_coord (v1 t)
      z3 = z_coord (v2 t)
      z4 = z_coord (v3 t)

-- Computed using the formula from Lai & Schumaker, Definition 15.4,
-- page 436.
volume :: Tetrahedron -> Double
volume t
       | (v0 t) == (v1 t) = 0
       | (v0 t) == (v2 t) = 0
       | (v0 t) == (v3 t) = 0
       | (v1 t) == (v2 t) = 0
       | (v1 t) == (v3 t) = 0
       | (v2 t) == (v3 t) = 0
       | otherwise = (1/6)*(det (vol_matrix t))


b0 :: Tetrahedron -> (RealFunction Point)
b0 t point = (volume inner_tetrahedron) / (volume t)
             where
               inner_tetrahedron = t { v0 = point }

b1 :: Tetrahedron -> (RealFunction Point)
b1 t point = (volume inner_tetrahedron) / (volume t)
             where
               inner_tetrahedron = t { v1 = point }

b2 :: Tetrahedron -> (RealFunction Point)
b2 t point = (volume inner_tetrahedron) / (volume t)
             where
               inner_tetrahedron = t { v2 = point }

b3 :: Tetrahedron -> (RealFunction Point)
b3 t point = (volume inner_tetrahedron) / (volume t)
             where
               inner_tetrahedron = t { v3 = point }