src/Tetrahedron.hs

   1 module Tetrahedron
   2 where
   3
   4 import Numeric.LinearAlgebra hiding (i, scale)
   5
   6 import Cube (back,
   7              Cube(datum),
   8              down,
   9              front,
  10              Gridded,
  11              left,
  12              right,
  13              top)
  14 import Misc (factorial)
  15 import Point
  16 import RealFunction
  17 import ThreeDimensional
  18
  19 data Tetrahedron = Tetrahedron { cube :: Cube,
  20                                  v0   :: Point,
  21                                  v1   :: Point,
  22                                  v2   :: Point,
  23                                  v3   :: Point }
  24                    deriving (Eq)
  25
  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)) ++
  30              ")\n\n"
  31
  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)
  38     top t = top (cube t)
  39
  40
  41 instance ThreeDimensional Tetrahedron where
  42     center t = ((v0 t) + (v1 t) + (v2 t) + (v3 t)) `scale` (1/4)
  43     contains_point t p =
  44         (b0 t p) >= 0 && (b1 t p) >= 0 && (b2 t p) >= 0 && (b3 t p) >= 0
  45
  46
  47 polynomial :: Tetrahedron -> (RealFunction Point)
  48 polynomial t =
  49     sum [ (c t i j k l) `cmult` (beta t i j k l) | i <- [0..3],
  50                                                    j <- [0..3],
  51                                                    k <- [0..3],
  52                                                    l <- [0..3],
  53                                                    i + j + k + l == 3]
  54
  55
  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
  59 xi = domain_point
  60
  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"
  66    where
  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'
  72
  73
  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)
  77 beta t i j k l
  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"
  81   where
  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
  88
  89
  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)
  92     where
  93       f = front x
  94       flt = front (left (top x))
  95       fl = front (left x)
  96       ft =  front (top x)
  97       i = cube x
  98       l = left x
  99       lt = left (top x)
 100       t = top x
 101
 102
 103 c x 0 0 0 3 = datum $ (1/8) * (i + f + r + t + rt + fr + ft + frt)
 104     where
 105       f = front x
 106       fr = front (right x)
 107       frt = front (right (top x))
 108       ft = front (top x)
 109       i = cube x
 110       r = right x
 111       rt = right (top x)
 112       t = top x
 113
 114
 115 c x 0 0 2 1 = datum $ (5/24)*(i + f + t + ft) + (1/24)*(l + fl + lt + flt)
 116     where
 117       f = front x
 118       flt = front (left (top x))
 119       fl = front (left x)
 120       ft = front (top x)
 121       i = cube x
 122       l = left x
 123       lt = left (top x)
 124       t = top x
 125
 126
 127 c x 0 0 1 2 = datum $ (5/24)*(i + f + t + ft) + (1/24)*(r + fr + rt + frt)
 128     where
 129       f = front x
 130       frt = front (right (top x))
 131       fr = front (right x)
 132       ft = front (top x)
 133       i = cube x
 134       r = right x
 135       rt = right (top x)
 136       t = top x
 137
 138
 139 c x 0 1 2 0 = datum $
 140               (5/24)*(i + f) +
 141               (1/8)*(l + t + fl + ft) +
 142               (1/24)*(lt + flt)
 143     where
 144       f = front x
 145       flt = front (left (top x))
 146       fl = front (left x)
 147       ft = front (top x)
 148       i = cube x
 149       l = left x
 150       lt = left (top x)
 151       t = top x
 152
 153
 154 c x 0 1 0 2 = datum $
 155               (5/24)*(i + f) +
 156               (1/8)*(r + t + fr + ft) +
 157               (1/24)*(rt + frt)
 158     where
 159       f = front x
 160       fr = front (right x)
 161       frt = front (right (top x))
 162       ft = front (top x)
 163       i = cube x
 164       r = right x
 165       rt = right (top x)
 166       t = top x
 167
 168
 169 c x 0 1 1 1 = datum $
 170           (13/48)*(i + f) +
 171           (7/48)*(t + ft) +
 172           (1/32)*(l + r + fl + fr) +
 173           (1/96)*(lt + rt + flt + frt)
 174     where
 175       f = front x
 176       flt = front (left (top x))
 177       fl = front (left x)
 178       fr = front (right x)
 179       frt = front (right (top x))
 180       ft = front (top x)
 181       i = cube x
 182       l = left x
 183       lt = left (top x)
 184       r = right x
 185       rt = right (top x)
 186       t = top x
 187
 188
 189 c x 0 2 1 0 = datum $
 190           (13/48)*(i + f) +
 191           (17/192)*(l + t + fl + ft) +
 192           (1/96)*(lt + flt) +
 193           (1/64)*(r + d + fr + fd) +
 194           (1/192)*(rt + ld + frt + fld)
 195     where
 196       d = down x
 197       f = front x
 198       fd = front (down x)
 199       fld = front (left (down x))
 200       flt = front (left (top x))
 201       fl = front (left x)
 202       fr = front (right x)
 203       frt = front (right (top x))
 204       ft = front (top x)
 205       i = cube x
 206       l = left x
 207       ld = left (down x)
 208       lt = left (top x)
 209       r = right x
 210       rt = right (top x)
 211       t = top x
 212
 213
 214 c x 0 2 0 1 = datum $
 215           (13/48)*(i + f) +
 216           (17/192)*(r + t + fr + ft) +
 217           (1/96)*(rt + frt) +
 218           (1/64)*(l + d + fl + fd) +
 219           (1/192)*(rd + lt + flt + frd)
 220     where
 221       d = down x
 222       f = front x
 223       fd = front (down x)
 224       flt = front (left (top x))
 225       fl = front (left x)
 226       frd = front (right (down x))
 227       fr = front (right x)
 228       frt = front (right (top x))
 229       ft = front (top x)
 230       i = cube x
 231       l = left x
 232       lt = left (top x)
 233       r = right x
 234       rd = right (down x)
 235       rt = right (top x)
 236       t = top x
 237
 238
 239 c x 0 3 0 0 = datum $
 240           (13/48)*(i + f) +
 241           (5/96)*(l + r + t + d + fl + fr + ft + fd) +
 242           (1/192)*(rt + rd + lt + ld + frt + frd + flt + fld)
 243     where
 244       d = down x
 245       f = front x
 246       fd = front (down x)
 247       fld = front (left (down x))
 248       flt = front (left (top x))
 249       fl = front (left x)
 250       frd = front (right (down x))
 251       fr = front (right x)
 252       frt = front (right (top x))
 253       ft = front (top x)
 254       i = cube x
 255       l = left x
 256       ld = left (down x)
 257       lt = left (top x)
 258       r = right x
 259       rd = right (down x)
 260       rt = right (top x)
 261       t = top x
 262
 263 c x 1 0 2 0 = datum $ (1/4)*i + (1/6)*(f + l + t) + (1/12)*(lt + fl + ft)
 264     where
 265       f = front x
 266       fl = front (left x)
 267       ft = front (top x)
 268       i = cube x
 269       l = left x
 270       lt = left (top x)
 271       t = top x
 272
 273
 274 c x 1 0 0 2 = datum $ (1/4)*i + (1/6)*(f + r + t) + (1/12)*(rt + fr + ft)
 275     where
 276       f = front x
 277       fr = front (right x)
 278       ft = front (top x)
 279       i = cube x
 280       r = right x
 281       rt = right (top x)
 282       t = top x
 283
 284
 285 c x 1 0 1 1 = datum $
 286           (1/3)*i +
 287           (5/24)*(f + t) +
 288           (1/12)*ft +
 289           (1/24)*(l + r) +
 290           (1/48)*(lt + rt + fl + fr)
 291     where
 292       f = front x
 293       fl = front (left x)
 294       fr = front (right x)
 295       ft = front (top x)
 296       i = cube x
 297       l = left x
 298       lt = left (top x)
 299       r = right x
 300       rt = right (top x)
 301       t = top x
 302
 303
 304 c x 1 1 1 0 = datum $
 305           (1/3)*i +
 306           (5/24)*f +
 307           (1/8)*(l + t) +
 308           (5/96)*(fl + ft) +
 309           (1/48)*(d + r + lt) +
 310           (1/96)*(fd + ld + rt + fr)
 311     where
 312       d = down x
 313       f = front x
 314       fd = front (down x)
 315       fl = front (left x)
 316       fr = front (right x)
 317       ft = front (top x)
 318       i = cube x
 319       l = left x
 320       ld = left (down x)
 321       lt = left (top x)
 322       r = right x
 323       rt = right (top x)
 324       t = top x
 325
 326
 327 c x 1 1 0 1 = datum $
 328           (1/3)*i +
 329           (5/24)*f +
 330           (1/8)*(r + t) +
 331           (5/96)*(fr + ft) +
 332           (1/48)*(d + l + rt) +
 333           (1/96)*(fd + lt + rd + fl)
 334     where
 335       d = down x
 336       f = front x
 337       fd = front (down x)
 338       fl = front (left x)
 339       fr = front (right x)
 340       ft = front (top x)
 341       i = cube x
 342       l = left x
 343       lt = left (top x)
 344       r = right x
 345       rd = right (down x)
 346       rt = right (top x)
 347       t = top x
 348
 349
 350
 351 c x 1 2 0 0 = datum $
 352           (1/3)*i +
 353           (5/24)*f +
 354           (7/96)*(l + r + t + d) +
 355           (1/32)*(fl + fr + ft + fd) +
 356           (1/96)*(rt + rd + lt + ld)
 357     where
 358       d = down x
 359       f = front x
 360       fd = front (down x)
 361       fl = front (left x)
 362       fr = front (right x)
 363       ft = front (top x)
 364       i = cube x
 365       l = left x
 366       ld = left (down x)
 367       lt = left (top x)
 368       r = right x
 369       rd = right (down x)
 370       rt = right (top x)
 371       t = top x
 372
 373
 374 c x 2 0 1 0 = datum $
 375           (3/8)*i +
 376           (7/48)*(f + t + l) +
 377           (1/48)*(r + d + b + lt + fl + ft) +
 378           (1/96)*(rt + bt + fr + fd + ld + bl)
 379     where
 380       b = back x
 381       bl = back (left x)
 382       bt = back (top x)
 383       d = down x
 384       f = front x
 385       fd = front (down x)
 386       fl = front (left x)
 387       fr = front (right x)
 388       ft = front (top x)
 389       i = cube x
 390       l = left x
 391       ld = left (down x)
 392       lt = left (top x)
 393       r = right x
 394       rt = right (top x)
 395       t = top x
 396
 397
 398 c x 2 0 0 1 = datum $
 399           (3/8)*i +
 400           (7/48)*(f + t + r) +
 401           (1/48)*(l + d + b + rt + fr + ft) +
 402           (1/96)*(lt + bt + fl + fd + rd + br)
 403     where
 404       b = back x
 405       br = back (right x)
 406       bt = back (top x)
 407       d = down x
 408       f = front x
 409       fd = front (down x)
 410       fl = front (left x)
 411       fr = front (right x)
 412       ft = front (top x)
 413       i = cube x
 414       l = left x
 415       lt = left (top x)
 416       r = right x
 417       rd = right (down x)
 418       rt = right (top x)
 419       t = top x
 420
 421
 422
 423 c x 2 1 0 0 = datum $
 424           (3/8)*i +
 425           (1/12)*(t + r + l + d) +
 426           (1/64)*(ft + fr + fl + fd) +
 427           (7/48)*f +
 428           (1/48)*b +
 429           (1/96)*(rt + ld + lt + rd) +
 430           (1/192)*(bt + br + bl + bd)
 431     where
 432       b = back x
 433       bd = back (down x)
 434       bl = back (left x)
 435       br = back (right x)
 436       bt = back (top x)
 437       d = down x
 438       f = front x
 439       fd = front (down x)
 440       fl = front (left x)
 441       fr = front (right x)
 442       ft = front (top x)
 443       i = cube x
 444       l = left x
 445       ld = left (down x)
 446       lt = left (top x)
 447       r = right x
 448       rd = right (down x)
 449       rt = right (top x)
 450       t = top x
 451
 452
 453 c x 3 0 0 0 = datum $
 454           (3/8)*i +
 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)
 458     where
 459       b = back x
 460       bd = back (down x)
 461       bl = back (left x)
 462       br = back (right x)
 463       bt = back (top x)
 464       d = down x
 465       f = front x
 466       fd = front (down x)
 467       fl = front (left x)
 468       fr = front (right x)
 469       ft = front (top x)
 470       i = cube x
 471       l = left x
 472       ld = left (down x)
 473       lt = left (top x)
 474       r = right x
 475       rd = right (down x)
 476       rt = right (top x)
 477       t = top x
 478
 479
 480 c _ _ _ _ _ = error "coefficient index out of bounds"
 481
 482
 483
 484 vol_matrix :: Tetrahedron -> Matrix Double
 485 vol_matrix t = (4><4) $
 486                [1,  1,  1,  1,
 487                 x1, x2, x3, x4,
 488                 y1, y2, y3, y4,
 489                 z1, z2, z3, z4 ]
 490     where
 491       x1 = x_coord (v0 t)
 492       x2 = x_coord (v1 t)
 493       x3 = x_coord (v2 t)
 494       x4 = x_coord (v3 t)
 495       y1 = y_coord (v0 t)
 496       y2 = y_coord (v1 t)
 497       y3 = y_coord (v2 t)
 498       y4 = y_coord (v3 t)
 499       z1 = z_coord (v0 t)
 500       z2 = z_coord (v1 t)
 501       z3 = z_coord (v2 t)
 502       z4 = z_coord (v3 t)
 503
 504 -- Computed using the formula from Lai & Schumaker, Definition 15.4,
 505 -- page 436.
 506 volume :: Tetrahedron -> Double
 507 volume t
 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))
 515
 516
 517 b0 :: Tetrahedron -> (RealFunction Point)
 518 b0 t point = (volume inner_tetrahedron) / (volume t)
 519              where
 520                inner_tetrahedron = t { v0 = point }
 521
 522 b1 :: Tetrahedron -> (RealFunction Point)
 523 b1 t point = (volume inner_tetrahedron) / (volume t)
 524              where
 525                inner_tetrahedron = t { v1 = point }
 526
 527 b2 :: Tetrahedron -> (RealFunction Point)
 528 b2 t point = (volume inner_tetrahedron) / (volume t)
 529              where
 530                inner_tetrahedron = t { v2 = point }
 531
 532 b3 :: Tetrahedron -> (RealFunction Point)
 533 b3 t point = (volume inner_tetrahedron) / (volume t)
 534              where
 535                inner_tetrahedron = t { v3 = point }