module Tetrahedron
where

import Numeric.LinearAlgebra hiding (i, scale)
import Prelude hiding (LT)
import Test.QuickCheck (Arbitrary(..), Gen)

import Cardinal
import FunctionValues
import Misc (factorial)
import Point
import RealFunction
import ThreeDimensional

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


instance Arbitrary Tetrahedron where
    arbitrary = do
      rnd_v0 <- arbitrary :: Gen Point
      rnd_v1 <- arbitrary :: Gen Point
      rnd_v2 <- arbitrary :: Gen Point
      rnd_v3 <- arbitrary :: Gen Point
      rnd_fv <- arbitrary :: Gen FunctionValues
      return (Tetrahedron rnd_fv rnd_v0 rnd_v1 rnd_v2 rnd_v3)


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


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.
--   Simply an alias for the domain_point function.
xi :: Tetrahedron -> Int -> Int -> Int -> Int -> Point
xi = domain_point

-- | Returns the domain point of t with indices i,j,k,l.
domain_point :: Tetrahedron -> Int -> Int -> Int -> Int -> Point
domain_point t i j k l
   | i + j + k + l == 3 = weighted_sum `scale` (1/3)
   | otherwise = error "domain point 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


-- | The coefficient function. c t i j k l returns the coefficient
--   c_ijkl with respect to the tetrahedron t. The definition uses
--   pattern matching to mimic the definitions given in Sorokina and
--   Zeilfelder, pp. 84-86. If incorrect indices are supplied, the
--   function will simply error.
c :: Tetrahedron -> Int -> Int -> Int -> Int -> Double
c t 0 0 3 0 = eval (fv t) $
              (1/8) * (I + F + L + T + LT + FL + FT + FLT)

c t 0 0 0 3 = eval (fv t) $
              (1/8) * (I + F + R + T + RT + FR + FT + FRT)

c t 0 0 2 1 = eval (fv t) $
              (5/24)*(I + F + T + FT) +
              (1/24)*(L + FL + LT + FLT)

c t 0 0 1 2 = eval (fv t) $
              (5/24)*(I + F + T + FT) +
              (1/24)*(R + FR + RT + FRT)

c t 0 1 2 0 = eval (fv t) $
              (5/24)*(I + F) +
              (1/8)*(L + T + FL + FT) +
              (1/24)*(LT + FLT)

c t 0 1 0 2 = eval (fv t) $
              (5/24)*(I + F) +
              (1/8)*(R + T + FR + FT) +
              (1/24)*(RT + FRT)

c t 0 1 1 1 = eval (fv t) $
          (13/48)*(I + F) +
          (7/48)*(T + FT) +
          (1/32)*(L + R + FL + FR) +
          (1/96)*(LT + RT + FLT + FRT)

c t 0 2 1 0 = eval (fv t) $
          (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)

c t 0 2 0 1 = eval (fv t) $
          (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)

c t 0 3 0 0 = eval (fv t) $
          (13/48)*(I + F) +
          (5/96)*(L + R + T + D + FL + FR + FT + FD) +
          (1/192)*(RT + RD + LT + LD + FRT + FRD + FLT + FLD)

c t 1 0 2 0 = eval (fv t) $
              (1/4)*I +
              (1/6)*(F + L + T) +
              (1/12)*(LT + FL + FT)

c t 1 0 0 2 = eval (fv t) $
              (1/4)*I +
              (1/6)*(F + R + T) +
              (1/12)*(RT + FR + FT)

c t 1 0 1 1 = eval (fv t) $
          (1/3)*I +
          (5/24)*(F + T) +
          (1/12)*FT +
          (1/24)*(L + R) +
          (1/48)*(LT + RT + FL + FR)

c t 1 1 1 0 = eval (fv t) $
          (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)

c t 1 1 0 1 = eval (fv t) $
          (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)

c t 1 2 0 0 = eval (fv t) $
          (1/3)*I +
          (5/24)*F +
          (7/96)*(L + R + T + D) +
          (1/32)*(FL + FR + FT + FD) +
          (1/96)*(RT + RD + LT + LD)

c t 2 0 1 0 = eval (fv t) $
          (3/8)*I +
          (7/48)*(F + T + L) +
          (1/48)*(R + D + B + LT + FL + FT) +
          (1/96)*(RT + BT + FR + FD + LD + BL)

c t 2 0 0 1 = eval (fv t) $
          (3/8)*I +
          (7/48)*(F + T + R) +
          (1/48)*(L + D + B + RT + FR + FT) +
          (1/96)*(LT + BT + FL + FD + RD + BR)

c t 2 1 0 0 = eval (fv t) $
          (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)

c t 3 0 0 0 = eval (fv t) $
          (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)

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



-- | The matrix used in the tetrahedron volume calculation as given in
--   Lai & Schumaker, Definition 15.4, page 436.
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))


-- | The barycentric coordinates of a point with respect to v0.
b0 :: Tetrahedron -> (RealFunction Point)
b0 t point = (volume inner_tetrahedron) / (volume t)
             where
               inner_tetrahedron = t { v0 = point }


-- | The barycentric coordinates of a point with respect to v1.
b1 :: Tetrahedron -> (RealFunction Point)
b1 t point = (volume inner_tetrahedron) / (volume t)
             where
               inner_tetrahedron = t { v1 = point }


-- | The barycentric coordinates of a point with respect to v2.
b2 :: Tetrahedron -> (RealFunction Point)
b2 t point = (volume inner_tetrahedron) / (volume t)
             where
               inner_tetrahedron = t { v2 = point }


-- | The barycentric coordinates of a point with respect to v3.
b3 :: Tetrahedron -> (RealFunction Point)
b3 t point = (volume inner_tetrahedron) / (volume t)
             where
               inner_tetrahedron = t { v3 = point }