4 import Numeric.LinearAlgebra hiding (i, scale)
14 import Misc (factorial)
17 import ThreeDimensional
19 data Tetrahedron = Tetrahedron { cube :: Cube,
26 instance Show Tetrahedron where
27 show t = "Tetrahedron (Cube: " ++ (show (cube t)) ++ ") " ++
28 "(v0: " ++ (show (v0 t)) ++ ") (v1: " ++ (show (v1 t)) ++
29 ") (v2: " ++ (show (v2 t)) ++ ") (v3: " ++ (show (v3 t)) ++
32 instance Gridded Tetrahedron where
33 back t = back (cube t)
34 down t = down (cube t)
35 front t = front (cube t)
36 left t = left (cube t)
37 right t = right (cube t)
41 instance ThreeDimensional Tetrahedron where
42 center t = ((v0 t) + (v1 t) + (v2 t) + (v3 t)) `scale` (1/4)
44 (b0 t p) >= 0 && (b1 t p) >= 0 && (b2 t p) >= 0 && (b3 t p) >= 0
47 polynomial :: Tetrahedron -> (RealFunction Point)
49 sum [ (c t i j k l) `cmult` (beta t i j k l) | i <- [0..3],
56 -- | Returns the domain point of t with indices i,j,k,l.
57 -- Simply an alias for the domain_point function.
58 xi :: Tetrahedron -> Int -> Int -> Int -> Int -> Point
61 -- | Returns the domain point of t with indices i,j,k,l.
62 domain_point :: Tetrahedron -> Int -> Int -> Int -> Int -> Point
63 domain_point t i j k l
64 | i + j + k + l == 3 = weighted_sum `scale` (1/3)
65 | otherwise = error "domain point index out of bounds"
67 v0' = (v0 t) `scale` (fromIntegral i)
68 v1' = (v1 t) `scale` (fromIntegral j)
69 v2' = (v2 t) `scale` (fromIntegral k)
70 v3' = (v3 t) `scale` (fromIntegral l)
71 weighted_sum = v0' + v1' + v2' + v3'
74 -- | The Bernstein polynomial on t with indices i,j,k,l. Denoted by a
75 -- capital 'B' in the Sorokina/Zeilfelder paper.
76 beta :: Tetrahedron -> Int -> Int -> Int -> Int -> (RealFunction Point)
78 | (i + j + k + l == 3) =
79 coefficient `cmult` (b0_term * b1_term * b2_term * b3_term)
80 | otherwise = error "basis function index out of bounds"
82 denominator = (factorial i)*(factorial j)*(factorial k)*(factorial l)
83 coefficient = 6 / (fromIntegral denominator)
84 b0_term = (b0 t) `fexp` i
85 b1_term = (b1 t) `fexp` j
86 b2_term = (b2 t) `fexp` k
87 b3_term = (b3 t) `fexp` l
90 c :: Tetrahedron -> Int -> Int -> Int -> Int -> Double
91 c x 0 0 3 0 = datum $ (1/8) * (i + f + l + t + lt + fl + ft + flt)
94 flt = front (left (top x))
103 c x 0 0 0 3 = datum $ (1/8) * (i + f + r + t + rt + fr + ft + frt)
107 frt = front (right (top x))
115 c x 0 0 2 1 = datum $ (5/24)*(i + f + t + ft) + (1/24)*(l + fl + lt + flt)
118 flt = front (left (top x))
127 c x 0 0 1 2 = datum $ (5/24)*(i + f + t + ft) + (1/24)*(r + fr + rt + frt)
130 frt = front (right (top x))
139 c x 0 1 2 0 = datum $
141 (1/8)*(l + t + fl + ft) +
145 flt = front (left (top x))
154 c x 0 1 0 2 = datum $
156 (1/8)*(r + t + fr + ft) +
161 frt = front (right (top x))
169 c x 0 1 1 1 = datum $
172 (1/32)*(l + r + fl + fr) +
173 (1/96)*(lt + rt + flt + frt)
176 flt = front (left (top x))
179 frt = front (right (top x))
189 c x 0 2 1 0 = datum $
191 (17/192)*(l + t + fl + ft) +
193 (1/64)*(r + d + fr + fd) +
194 (1/192)*(rt + ld + frt + fld)
199 fld = front (left (down x))
200 flt = front (left (top x))
203 frt = front (right (top x))
214 c x 0 2 0 1 = datum $
216 (17/192)*(r + t + fr + ft) +
218 (1/64)*(l + d + fl + fd) +
219 (1/192)*(rd + lt + flt + frd)
224 flt = front (left (top x))
226 frd = front (right (down x))
228 frt = front (right (top x))
239 c x 0 3 0 0 = datum $
241 (5/96)*(l + r + t + d + fl + fr + ft + fd) +
242 (1/192)*(rt + rd + lt + ld + frt + frd + flt + fld)
247 fld = front (left (down x))
248 flt = front (left (top x))
250 frd = front (right (down x))
252 frt = front (right (top x))
263 c x 1 0 2 0 = datum $ (1/4)*i + (1/6)*(f + l + t) + (1/12)*(lt + fl + ft)
274 c x 1 0 0 2 = datum $ (1/4)*i + (1/6)*(f + r + t) + (1/12)*(rt + fr + ft)
285 c x 1 0 1 1 = datum $
290 (1/48)*(lt + rt + fl + fr)
304 c x 1 1 1 0 = datum $
309 (1/48)*(d + r + lt) +
310 (1/96)*(fd + ld + rt + fr)
327 c x 1 1 0 1 = datum $
332 (1/48)*(d + l + rt) +
333 (1/96)*(fd + lt + rd + fl)
351 c x 1 2 0 0 = datum $
354 (7/96)*(l + r + t + d) +
355 (1/32)*(fl + fr + ft + fd) +
356 (1/96)*(rt + rd + lt + ld)
374 c x 2 0 1 0 = datum $
377 (1/48)*(r + d + b + lt + fl + ft) +
378 (1/96)*(rt + bt + fr + fd + ld + bl)
398 c x 2 0 0 1 = datum $
401 (1/48)*(l + d + b + rt + fr + ft) +
402 (1/96)*(lt + bt + fl + fd + rd + br)
423 c x 2 1 0 0 = datum $
425 (1/12)*(t + r + l + d) +
426 (1/64)*(ft + fr + fl + fd) +
429 (1/96)*(rt + ld + lt + rd) +
430 (1/192)*(bt + br + bl + bd)
453 c x 3 0 0 0 = datum $
455 (1/12)*(t + f + l + r + d + b) +
456 (1/96)*(lt + fl + ft + rt + bt + fr) +
457 (1/96)*(fd + ld + bd + br + rd + bl)
480 c _ _ _ _ _ = error "coefficient index out of bounds"
484 vol_matrix :: Tetrahedron -> Matrix Double
485 vol_matrix t = (4><4) $
504 -- Computed using the formula from Lai & Schumaker, Definition 15.4,
506 volume :: Tetrahedron -> Double
508 | (v0 t) == (v1 t) = 0
509 | (v0 t) == (v2 t) = 0
510 | (v0 t) == (v3 t) = 0
511 | (v1 t) == (v2 t) = 0
512 | (v1 t) == (v3 t) = 0
513 | (v2 t) == (v3 t) = 0
514 | otherwise = (1/6)*(det (vol_matrix t))
517 b0 :: Tetrahedron -> (RealFunction Point)
518 b0 t point = (volume inner_tetrahedron) / (volume t)
520 inner_tetrahedron = t { v0 = point }
522 b1 :: Tetrahedron -> (RealFunction Point)
523 b1 t point = (volume inner_tetrahedron) / (volume t)
525 inner_tetrahedron = t { v1 = point }
527 b2 :: Tetrahedron -> (RealFunction Point)
528 b2 t point = (volume inner_tetrahedron) / (volume t)
530 inner_tetrahedron = t { v2 = point }
532 b3 :: Tetrahedron -> (RealFunction Point)
533 b3 t point = (volume inner_tetrahedron) / (volume t)
535 inner_tetrahedron = t { v3 = point }