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 xi :: Tetrahedron -> Int -> Int -> Int -> Int -> Point
59 | i + j + k + l == 3 = weighted_sum `scale` (1/3)
60 | otherwise = error "xi index out of bounds"
62 v0' = (v0 t) `scale` (fromIntegral i)
63 v1' = (v1 t) `scale` (fromIntegral j)
64 v2' = (v2 t) `scale` (fromIntegral k)
65 v3' = (v3 t) `scale` (fromIntegral l)
66 weighted_sum = v0' + v1' + v2' + v3'
70 -- | The Bernstein polynomial on t with indices i,j,k,l. Denoted by a
71 -- capital 'B' in the Sorokina/Zeilfelder paper.
72 beta :: Tetrahedron -> Int -> Int -> Int -> Int -> (RealFunction Point)
74 | (i + j + k + l == 3) =
75 coefficient `cmult` (b0_term * b1_term * b2_term * b3_term)
76 | otherwise = error "basis function index out of bounds"
78 denominator = (factorial i)*(factorial j)*(factorial k)*(factorial l)
79 coefficient = 6 / (fromIntegral denominator)
80 b0_term = (b0 t) `fexp` i
81 b1_term = (b1 t) `fexp` j
82 b2_term = (b2 t) `fexp` k
83 b3_term = (b3 t) `fexp` l
86 c :: Tetrahedron -> Int -> Int -> Int -> Int -> Double
87 c x 0 0 3 0 = datum $ (1/8) * (i + f + l + t + lt + fl + ft + flt)
90 flt = front (left (top x))
99 c x 0 0 0 3 = datum $ (1/8) * (i + f + r + t + rt + fr + ft + frt)
103 frt = front (right (top x))
111 c x 0 0 2 1 = datum $ (5/24)*(i + f + t + ft) + (1/24)*(l + fl + lt + flt)
114 flt = front (left (top x))
123 c x 0 0 1 2 = datum $ (5/24)*(i + f + t + ft) + (1/24)*(r + fr + rt + frt)
126 frt = front (right (top x))
135 c x 0 1 2 0 = datum $
137 (1/8)*(l + t + fl + ft) +
141 flt = front (left (top x))
150 c x 0 1 0 2 = datum $
152 (1/8)*(r + t + fr + ft) +
157 frt = front (right (top x))
165 c x 0 1 1 1 = datum $
168 (1/32)*(l + r + fl + fr) +
169 (1/96)*(lt + rt + flt + frt)
172 flt = front (left (top x))
175 frt = front (right (top x))
185 c x 0 2 1 0 = datum $
187 (17/192)*(l + t + fl + ft) +
189 (1/64)*(r + d + fr + fd) +
190 (1/192)*(rt + ld + frt + fld)
195 fld = front (left (down x))
196 flt = front (left (top x))
199 frt = front (right (top x))
210 c x 0 2 0 1 = datum $
212 (17/192)*(r + t + fr + ft) +
214 (1/64)*(l + d + fl + fd) +
215 (1/192)*(rd + lt + flt + frd)
220 flt = front (left (top x))
222 frd = front (right (down x))
224 frt = front (right (top x))
235 c x 0 3 0 0 = datum $
237 (5/96)*(l + r + t + d + fl + fr + ft + fd) +
238 (1/192)*(rt + rd + lt + ld + frt + frd + flt + fld)
243 fld = front (left (down x))
244 flt = front (left (top x))
246 frd = front (right (down x))
248 frt = front (right (top x))
259 c x 1 0 2 0 = datum $ (1/4)*i + (1/6)*(f + l + t) + (1/12)*(lt + fl + ft)
270 c x 1 0 0 2 = datum $ (1/4)*i + (1/6)*(f + r + t) + (1/12)*(rt + fr + ft)
281 c x 1 0 1 1 = datum $
286 (1/48)*(lt + rt + fl + fr)
300 c x 1 1 1 0 = datum $
305 (1/48)*(d + r + lt) +
306 (1/96)*(fd + ld + rt + fr)
323 c x 1 1 0 1 = datum $
328 (1/48)*(d + l + rt) +
329 (1/96)*(fd + lt + rd + fl)
347 c x 1 2 0 0 = datum $
350 (7/96)*(l + r + t + d) +
351 (1/32)*(fl + fr + ft + fd) +
352 (1/96)*(rt + rd + lt + ld)
370 c x 2 0 1 0 = datum $
373 (1/48)*(r + d + b + lt + fl + ft) +
374 (1/96)*(rt + bt + fr + fd + ld + bl)
394 c x 2 0 0 1 = datum $
397 (1/48)*(l + d + b + rt + fr + ft) +
398 (1/96)*(lt + bt + fl + fd + rd + br)
419 c x 2 1 0 0 = datum $
421 (1/12)*(t + r + l + d) +
422 (1/64)*(ft + fr + fl + fd) +
425 (1/96)*(rt + ld + lt + rd) +
426 (1/192)*(bt + br + bl + bd)
449 c x 3 0 0 0 = datum $
451 (1/12)*(t + f + l + r + d + b) +
452 (1/96)*(lt + fl + ft + rt + bt + fr) +
453 (1/96)*(fd + ld + bd + br + rd + bl)
476 c _ _ _ _ _ = error "coefficient index out of bounds"
480 vol_matrix :: Tetrahedron -> Matrix Double
481 vol_matrix t = (4><4) $
500 -- Computed using the formula from Lai & Schumaker, Definition 15.4,
502 volume :: Tetrahedron -> Double
504 | (v0 t) == (v1 t) = 0
505 | (v0 t) == (v2 t) = 0
506 | (v0 t) == (v3 t) = 0
507 | (v1 t) == (v2 t) = 0
508 | (v1 t) == (v3 t) = 0
509 | (v2 t) == (v3 t) = 0
510 | otherwise = (1/6)*(det (vol_matrix t))
513 b0 :: Tetrahedron -> (RealFunction Point)
514 b0 t point = (volume inner_tetrahedron) / (volume t)
516 inner_tetrahedron = t { v0 = point }
518 b1 :: Tetrahedron -> (RealFunction Point)
519 b1 t point = (volume inner_tetrahedron) / (volume t)
521 inner_tetrahedron = t { v1 = point }
523 b2 :: Tetrahedron -> (RealFunction Point)
524 b2 t point = (volume inner_tetrahedron) / (volume t)
526 inner_tetrahedron = t { v2 = point }
528 b3 :: Tetrahedron -> (RealFunction Point)
529 b3 t point = (volume inner_tetrahedron) / (volume t)
531 inner_tetrahedron = t { v3 = point }