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