]> gitweb.michael.orlitzky.com - spline3.git/blob - Tetrahedron.hs
d468f0c570a920fe3967fdd7e20973ef2cd9945e
[spline3.git] / Tetrahedron.hs
1 module Tetrahedron
2 where
3
4 import Numeric.LinearAlgebra hiding (i, scale)
5 import Prelude hiding (LT)
6 import Test.QuickCheck (Arbitrary(..), Gen, choose)
7
8 import Cardinal
9 import Comparisons (nearly_ge)
10 import FunctionValues
11 import Misc (factorial)
12 import Point
13 import RealFunction
14 import ThreeDimensional
15
16 data Tetrahedron =
17 Tetrahedron { fv :: FunctionValues,
18 v0 :: Point,
19 v1 :: Point,
20 v2 :: Point,
21 v3 :: Point,
22 precomputed_volume :: Double,
23
24 -- | Between 0 and 23; used to quickly determine which
25 -- tetrahedron I am in the parent 'Cube' without
26 -- having to compare them all.
27 number :: Int
28 }
29 deriving (Eq)
30
31
32 instance Arbitrary Tetrahedron where
33 arbitrary = do
34 rnd_v0 <- arbitrary :: Gen Point
35 rnd_v1 <- arbitrary :: Gen Point
36 rnd_v2 <- arbitrary :: Gen Point
37 rnd_v3 <- arbitrary :: Gen Point
38 rnd_fv <- arbitrary :: Gen FunctionValues
39 rnd_no <- choose (0,23)
40
41 -- We can't assign an incorrect precomputed volume,
42 -- so we have to calculate the correct one here.
43 let t' = Tetrahedron rnd_fv rnd_v0 rnd_v1 rnd_v2 rnd_v3 0 rnd_no
44 let vol = volume t'
45 return (Tetrahedron rnd_fv rnd_v0 rnd_v1 rnd_v2 rnd_v3 vol rnd_no)
46
47
48 instance Show Tetrahedron where
49 show t = "Tetrahedron:\n" ++
50 " fv: " ++ (show (fv t)) ++ "\n" ++
51 " v0: " ++ (show (v0 t)) ++ "\n" ++
52 " v1: " ++ (show (v1 t)) ++ "\n" ++
53 " v2: " ++ (show (v2 t)) ++ "\n" ++
54 " v3: " ++ (show (v3 t)) ++ "\n"
55
56
57 instance ThreeDimensional Tetrahedron where
58 center t = ((v0 t) + (v1 t) + (v2 t) + (v3 t)) `scale` (1/4)
59 contains_point t p =
60 b0_unscaled `nearly_ge` 0 &&
61 b1_unscaled `nearly_ge` 0 &&
62 b2_unscaled `nearly_ge` 0 &&
63 b3_unscaled `nearly_ge` 0
64 where
65 -- Drop the useless division and volume calculation that we
66 -- would do if we used the regular b0,..b3 functions.
67 b0_unscaled :: Double
68 b0_unscaled = volume inner_tetrahedron
69 where inner_tetrahedron = t { v0 = p }
70
71 b1_unscaled :: Double
72 b1_unscaled = volume inner_tetrahedron
73 where inner_tetrahedron = t { v1 = p }
74
75 b2_unscaled :: Double
76 b2_unscaled = volume inner_tetrahedron
77 where inner_tetrahedron = t { v2 = p }
78
79 b3_unscaled :: Double
80 b3_unscaled = volume inner_tetrahedron
81 where inner_tetrahedron = t { v3 = p }
82
83
84 polynomial :: Tetrahedron -> (RealFunction Point)
85 polynomial t =
86 sum [ (c t i j k l) `cmult` (beta t i j k l) | i <- [0..3],
87 j <- [0..3],
88 k <- [0..3],
89 l <- [0..3],
90 i + j + k + l == 3]
91
92
93 -- | Returns the domain point of t with indices i,j,k,l.
94 -- Simply an alias for the domain_point function.
95 xi :: Tetrahedron -> Int -> Int -> Int -> Int -> Point
96 xi = domain_point
97
98 -- | Returns the domain point of t with indices i,j,k,l.
99 domain_point :: Tetrahedron -> Int -> Int -> Int -> Int -> Point
100 domain_point t i j k l
101 | i + j + k + l == 3 = weighted_sum `scale` (1/3)
102 | otherwise = error "domain point index out of bounds"
103 where
104 v0' = (v0 t) `scale` (fromIntegral i)
105 v1' = (v1 t) `scale` (fromIntegral j)
106 v2' = (v2 t) `scale` (fromIntegral k)
107 v3' = (v3 t) `scale` (fromIntegral l)
108 weighted_sum = v0' + v1' + v2' + v3'
109
110
111 -- | The Bernstein polynomial on t with indices i,j,k,l. Denoted by a
112 -- capital 'B' in the Sorokina/Zeilfelder paper.
113 beta :: Tetrahedron -> Int -> Int -> Int -> Int -> (RealFunction Point)
114 beta t i j k l
115 | (i + j + k + l == 3) =
116 coefficient `cmult` (b0_term * b1_term * b2_term * b3_term)
117 | otherwise = error "basis function index out of bounds"
118 where
119 denominator = (factorial i)*(factorial j)*(factorial k)*(factorial l)
120 coefficient = 6 / (fromIntegral denominator)
121 b0_term = (b0 t) `fexp` i
122 b1_term = (b1 t) `fexp` j
123 b2_term = (b2 t) `fexp` k
124 b3_term = (b3 t) `fexp` l
125
126
127 -- | The coefficient function. c t i j k l returns the coefficient
128 -- c_ijkl with respect to the tetrahedron t. The definition uses
129 -- pattern matching to mimic the definitions given in Sorokina and
130 -- Zeilfelder, pp. 84-86. If incorrect indices are supplied, the
131 -- function will simply error.
132 c :: Tetrahedron -> Int -> Int -> Int -> Int -> Double
133 c t 0 0 3 0 = eval (fv t) $
134 (1/8) * (I + F + L + T + LT + FL + FT + FLT)
135
136 c t 0 0 0 3 = eval (fv t) $
137 (1/8) * (I + F + R + T + RT + FR + FT + FRT)
138
139 c t 0 0 2 1 = eval (fv t) $
140 (5/24)*(I + F + T + FT) +
141 (1/24)*(L + FL + LT + FLT)
142
143 c t 0 0 1 2 = eval (fv t) $
144 (5/24)*(I + F + T + FT) +
145 (1/24)*(R + FR + RT + FRT)
146
147 c t 0 1 2 0 = eval (fv t) $
148 (5/24)*(I + F) +
149 (1/8)*(L + T + FL + FT) +
150 (1/24)*(LT + FLT)
151
152 c t 0 1 0 2 = eval (fv t) $
153 (5/24)*(I + F) +
154 (1/8)*(R + T + FR + FT) +
155 (1/24)*(RT + FRT)
156
157 c t 0 1 1 1 = eval (fv t) $
158 (13/48)*(I + F) +
159 (7/48)*(T + FT) +
160 (1/32)*(L + R + FL + FR) +
161 (1/96)*(LT + RT + FLT + FRT)
162
163 c t 0 2 1 0 = eval (fv t) $
164 (13/48)*(I + F) +
165 (17/192)*(L + T + FL + FT) +
166 (1/96)*(LT + FLT) +
167 (1/64)*(R + D + FR + FD) +
168 (1/192)*(RT + LD + FRT + FLD)
169
170 c t 0 2 0 1 = eval (fv t) $
171 (13/48)*(I + F) +
172 (17/192)*(R + T + FR + FT) +
173 (1/96)*(RT + FRT) +
174 (1/64)*(L + D + FL + FD) +
175 (1/192)*(RD + LT + FLT + FRD)
176
177 c t 0 3 0 0 = eval (fv t) $
178 (13/48)*(I + F) +
179 (5/96)*(L + R + T + D + FL + FR + FT + FD) +
180 (1/192)*(RT + RD + LT + LD + FRT + FRD + FLT + FLD)
181
182 c t 1 0 2 0 = eval (fv t) $
183 (1/4)*I +
184 (1/6)*(F + L + T) +
185 (1/12)*(LT + FL + FT)
186
187 c t 1 0 0 2 = eval (fv t) $
188 (1/4)*I +
189 (1/6)*(F + R + T) +
190 (1/12)*(RT + FR + FT)
191
192 c t 1 0 1 1 = eval (fv t) $
193 (1/3)*I +
194 (5/24)*(F + T) +
195 (1/12)*FT +
196 (1/24)*(L + R) +
197 (1/48)*(LT + RT + FL + FR)
198
199 c t 1 1 1 0 = eval (fv t) $
200 (1/3)*I +
201 (5/24)*F +
202 (1/8)*(L + T) +
203 (5/96)*(FL + FT) +
204 (1/48)*(D + R + LT) +
205 (1/96)*(FD + LD + RT + FR)
206
207 c t 1 1 0 1 = eval (fv t) $
208 (1/3)*I +
209 (5/24)*F +
210 (1/8)*(R + T) +
211 (5/96)*(FR + FT) +
212 (1/48)*(D + L + RT) +
213 (1/96)*(FD + LT + RD + FL)
214
215 c t 1 2 0 0 = eval (fv t) $
216 (1/3)*I +
217 (5/24)*F +
218 (7/96)*(L + R + T + D) +
219 (1/32)*(FL + FR + FT + FD) +
220 (1/96)*(RT + RD + LT + LD)
221
222 c t 2 0 1 0 = eval (fv t) $
223 (3/8)*I +
224 (7/48)*(F + T + L) +
225 (1/48)*(R + D + B + LT + FL + FT) +
226 (1/96)*(RT + BT + FR + FD + LD + BL)
227
228 c t 2 0 0 1 = eval (fv t) $
229 (3/8)*I +
230 (7/48)*(F + T + R) +
231 (1/48)*(L + D + B + RT + FR + FT) +
232 (1/96)*(LT + BT + FL + FD + RD + BR)
233
234 c t 2 1 0 0 = eval (fv t) $
235 (3/8)*I +
236 (1/12)*(T + R + L + D) +
237 (1/64)*(FT + FR + FL + FD) +
238 (7/48)*F +
239 (1/48)*B +
240 (1/96)*(RT + LD + LT + RD) +
241 (1/192)*(BT + BR + BL + BD)
242
243 c t 3 0 0 0 = eval (fv t) $
244 (3/8)*I +
245 (1/12)*(T + F + L + R + D + B) +
246 (1/96)*(LT + FL + FT + RT + BT + FR) +
247 (1/96)*(FD + LD + BD + BR + RD + BL)
248
249 c _ _ _ _ _ = error "coefficient index out of bounds"
250
251
252
253 -- | The matrix used in the tetrahedron volume calculation as given in
254 -- Lai & Schumaker, Definition 15.4, page 436.
255 vol_matrix :: Tetrahedron -> Matrix Double
256 vol_matrix t = (4><4)
257 [1, 1, 1, 1,
258 x1, x2, x3, x4,
259 y1, y2, y3, y4,
260 z1, z2, z3, z4 ]
261 where
262 (x1, y1, z1) = v0 t
263 (x2, y2, z2) = v1 t
264 (x3, y3, z3) = v2 t
265 (x4, y4, z4) = v3 t
266
267 -- | Computed using the formula from Lai & Schumaker, Definition 15.4,
268 -- page 436.
269 volume :: Tetrahedron -> Double
270 volume t
271 | (v0 t) == (v1 t) = 0
272 | (v0 t) == (v2 t) = 0
273 | (v0 t) == (v3 t) = 0
274 | (v1 t) == (v2 t) = 0
275 | (v1 t) == (v3 t) = 0
276 | (v2 t) == (v3 t) = 0
277 | otherwise = (1/6)*(det (vol_matrix t))
278
279
280 -- | The barycentric coordinates of a point with respect to v0.
281 b0 :: Tetrahedron -> (RealFunction Point)
282 b0 t point = (volume inner_tetrahedron) / (precomputed_volume t)
283 where
284 inner_tetrahedron = t { v0 = point }
285
286
287 -- | The barycentric coordinates of a point with respect to v1.
288 b1 :: Tetrahedron -> (RealFunction Point)
289 b1 t point = (volume inner_tetrahedron) / (precomputed_volume t)
290 where
291 inner_tetrahedron = t { v1 = point }
292
293
294 -- | The barycentric coordinates of a point with respect to v2.
295 b2 :: Tetrahedron -> (RealFunction Point)
296 b2 t point = (volume inner_tetrahedron) / (precomputed_volume t)
297 where
298 inner_tetrahedron = t { v2 = point }
299
300
301 -- | The barycentric coordinates of a point with respect to v3.
302 b3 :: Tetrahedron -> (RealFunction Point)
303 b3 t point = (volume inner_tetrahedron) / (precomputed_volume t)
304 where
305 inner_tetrahedron = t { v3 = point }