src/Tetrahedron.hs

   1 {-# LANGUAGE BangPatterns #-}
   2 module Tetrahedron (
   3   Tetrahedron(..),
   4   b0, -- Cube test
   5   b1, -- Cube test
   6   b2, -- Cube test
   7   b3, -- Cube test
   8   barycenter,
   9   c,
  10   contains_point,
  11   polynomial,
  12   tetrahedron_properties,
  13   tetrahedron_tests,
  14   volume -- Cube test
  15   )
  16 where
  17
  18 import qualified Data.Vector as V (
  19   singleton,
  20   snoc,
  21   sum
  22   )
  23
  24 import Test.Framework (Test, testGroup)
  25 import Test.Framework.Providers.HUnit (testCase)
  26 import Test.Framework.Providers.QuickCheck2 (testProperty)
  27 import Test.HUnit (Assertion, assertEqual)
  28 import Test.QuickCheck (Arbitrary(..), Gen, Property, (==>))
  29
  30 import Comparisons ((~=), nearly_ge)
  31 import FunctionValues (FunctionValues(..), empty_values)
  32 import Misc (factorial)
  33 import Point (Point(..), scale)
  34 import RealFunction (RealFunction, cmult, fexp)
  35
  36 data Tetrahedron =
  37   Tetrahedron { function_values :: FunctionValues,
  38                 v0 :: !Point,
  39                 v1 :: !Point,
  40                 v2 :: !Point,
  41                 v3 :: !Point,
  42                 precomputed_volume :: !Double
  43               }
  44     deriving (Eq)
  45
  46
  47 instance Arbitrary Tetrahedron where
  48     arbitrary = do
  49       rnd_v0 <- arbitrary :: Gen Point
  50       rnd_v1 <- arbitrary :: Gen Point
  51       rnd_v2 <- arbitrary :: Gen Point
  52       rnd_v3 <- arbitrary :: Gen Point
  53       rnd_fv <- arbitrary :: Gen FunctionValues
  54
  55       -- We can't assign an incorrect precomputed volume,
  56       -- so we have to calculate the correct one here.
  57       let t' = Tetrahedron rnd_fv rnd_v0 rnd_v1 rnd_v2 rnd_v3 0
  58       let vol = volume t'
  59       return (Tetrahedron rnd_fv rnd_v0 rnd_v1 rnd_v2 rnd_v3 vol)
  60
  61
  62 instance Show Tetrahedron where
  63     show t = "Tetrahedron:\n" ++
  64              "  function_values: " ++ (show (function_values t)) ++ "\n" ++
  65              "  v0: " ++ (show (v0 t)) ++ "\n" ++
  66              "  v1: " ++ (show (v1 t)) ++ "\n" ++
  67              "  v2: " ++ (show (v2 t)) ++ "\n" ++
  68              "  v3: " ++ (show (v3 t)) ++ "\n"
  69
  70
  71 -- | Find the barycenter of the given tetrahedron.
  72 --   We just average the four vertices.
  73 barycenter :: Tetrahedron -> Point
  74 barycenter (Tetrahedron _ v0' v1' v2' v3' _) =
  75   (v0' + v1' + v2' + v3') `scale` (1/4)
  76
  77 -- | A point is internal to a tetrahedron if all of its barycentric
  78 --   coordinates with respect to that tetrahedron are non-negative.
  79 contains_point :: Tetrahedron -> Point -> Bool
  80 contains_point t p0 =
  81   b0_unscaled `nearly_ge` 0 &&
  82   b1_unscaled `nearly_ge` 0 &&
  83   b2_unscaled `nearly_ge` 0 &&
  84   b3_unscaled `nearly_ge` 0
  85   where
  86     -- Drop the useless division and volume calculation that we
  87     -- would do if we used the regular b0,..b3 functions.
  88     b0_unscaled :: Double
  89     b0_unscaled = volume inner_tetrahedron
  90       where
  91         inner_tetrahedron = t { v0 = p0 }
  92
  93     b1_unscaled :: Double
  94     b1_unscaled = volume inner_tetrahedron
  95       where inner_tetrahedron = t { v1 = p0 }
  96
  97     b2_unscaled :: Double
  98     b2_unscaled = volume inner_tetrahedron
  99       where inner_tetrahedron = t { v2 = p0 }
 100
 101     b3_unscaled :: Double
 102     b3_unscaled = volume inner_tetrahedron
 103       where inner_tetrahedron = t { v3 = p0 }
 104
 105
 106 {-# INLINE polynomial #-}
 107 polynomial :: Tetrahedron -> (RealFunction Point)
 108 polynomial t =
 109     V.sum $ V.singleton ((c t 0 0 0 3) `cmult` (beta t 0 0 0 3)) `V.snoc`
 110             ((c t 0 0 1 2) `cmult` (beta t 0 0 1 2)) `V.snoc`
 111             ((c t 0 0 2 1) `cmult` (beta t 0 0 2 1)) `V.snoc`
 112             ((c t 0 0 3 0) `cmult` (beta t 0 0 3 0)) `V.snoc`
 113             ((c t 0 1 0 2) `cmult` (beta t 0 1 0 2)) `V.snoc`
 114             ((c t 0 1 1 1) `cmult` (beta t 0 1 1 1)) `V.snoc`
 115             ((c t 0 1 2 0) `cmult` (beta t 0 1 2 0)) `V.snoc`
 116             ((c t 0 2 0 1) `cmult` (beta t 0 2 0 1)) `V.snoc`
 117             ((c t 0 2 1 0) `cmult` (beta t 0 2 1 0)) `V.snoc`
 118             ((c t 0 3 0 0) `cmult` (beta t 0 3 0 0)) `V.snoc`
 119             ((c t 1 0 0 2) `cmult` (beta t 1 0 0 2)) `V.snoc`
 120             ((c t 1 0 1 1) `cmult` (beta t 1 0 1 1)) `V.snoc`
 121             ((c t 1 0 2 0) `cmult` (beta t 1 0 2 0)) `V.snoc`
 122             ((c t 1 1 0 1) `cmult` (beta t 1 1 0 1)) `V.snoc`
 123             ((c t 1 1 1 0) `cmult` (beta t 1 1 1 0)) `V.snoc`
 124             ((c t 1 2 0 0) `cmult` (beta t 1 2 0 0)) `V.snoc`
 125             ((c t 2 0 0 1) `cmult` (beta t 2 0 0 1)) `V.snoc`
 126             ((c t 2 0 1 0) `cmult` (beta t 2 0 1 0)) `V.snoc`
 127             ((c t 2 1 0 0) `cmult` (beta t 2 1 0 0)) `V.snoc`
 128             ((c t 3 0 0 0) `cmult` (beta t 3 0 0 0))
 129
 130
 131
 132 -- | The Bernstein polynomial on t with indices i,j,k,l. Denoted by a
 133 --   capital 'B' in the Sorokina/Zeilfelder paper.
 134 beta :: Tetrahedron -> Int -> Int -> Int -> Int -> (RealFunction Point)
 135 beta t i j k l =
 136   coefficient `cmult` (b0_term * b1_term * b2_term * b3_term)
 137   where
 138     denominator = (factorial i)*(factorial j)*(factorial k)*(factorial l)
 139     coefficient = 6 / (fromIntegral denominator)
 140     b0_term = (b0 t) `fexp` i
 141     b1_term = (b1 t) `fexp` j
 142     b2_term = (b2 t) `fexp` k
 143     b3_term = (b3 t) `fexp` l
 144
 145
 146 -- | The coefficient function. c t i j k l returns the coefficient
 147 --   c_ijkl with respect to the tetrahedron t. The definition uses
 148 --   pattern matching to mimic the definitions given in Sorokina and
 149 --   Zeilfelder, pp. 84-86. If incorrect indices are supplied, the world
 150 --   will end. This is for performance reasons.
 151 c :: Tetrahedron -> Int -> Int -> Int -> Int -> Double
 152 c !t !i !j !k !l =
 153   coefficient i j k l
 154   where
 155     fvs = function_values t
 156     f = front fvs
 157     b = back fvs
 158     r = right fvs
 159     l' = left fvs
 160     t' = top fvs
 161     d  = down fvs
 162     fl = front_left fvs
 163     fr = front_right fvs
 164     fd = front_down fvs
 165     ft = front_top fvs
 166     bl = back_left fvs
 167     br = back_right fvs
 168     bd = back_down fvs
 169     bt = back_top fvs
 170     ld = left_down fvs
 171     lt = left_top fvs
 172     rd = right_down fvs
 173     rt = right_top fvs
 174     fld = front_left_down fvs
 175     flt = front_left_top fvs
 176     frd = front_right_down fvs
 177     frt = front_right_top fvs
 178     i' = interior fvs
 179
 180     coefficient :: Int -> Int -> Int -> Int -> Double
 181     coefficient 0 0 3 0 =
 182       (1/8) * (i' + f + l' + t' + lt + fl + ft + flt)
 183
 184     coefficient 0 0 0 3 =
 185       (1/8) * (i' + f + r + t' + rt + fr + ft + frt)
 186
 187     coefficient 0 0 2 1 =
 188       (5/24)*(i' + f + t' + ft) + (1/24)*(l' + fl + lt + flt)
 189
 190     coefficient 0 0 1 2 =
 191       (5/24)*(i' + f + t' + ft) + (1/24)*(r + fr + rt + frt)
 192
 193     coefficient 0 1 2 0 =
 194       (5/24)*(i' + f) + (1/8)*(l' + t' + fl + ft)
 195                       + (1/24)*(lt + flt)
 196
 197     coefficient 0 1 0 2 =
 198       (5/24)*(i' + f) + (1/8)*(r + t' + fr + ft)
 199                       + (1/24)*(rt + frt)
 200
 201     coefficient 0 1 1 1 =
 202       (13/48)*(i' + f) + (7/48)*(t' + ft)
 203                        + (1/32)*(l' + r + fl + fr)
 204                        + (1/96)*(lt + rt + flt + frt)
 205
 206     coefficient 0 2 1 0 =
 207       (13/48)*(i' + f) + (17/192)*(l' + t' + fl + ft)
 208                        + (1/96)*(lt + flt)
 209                        + (1/64)*(r + d + fr + fd)
 210                        + (1/192)*(rt + ld + frt + fld)
 211
 212     coefficient 0 2 0 1 =
 213       (13/48)*(i' + f) + (17/192)*(r + t' + fr + ft)
 214                        + (1/96)*(rt + frt)
 215                        + (1/64)*(l' + d + fl + fd)
 216                        + (1/192)*(rd + lt + flt + frd)
 217
 218     coefficient 0 3 0 0 =
 219      (13/48)*(i' + f) + (5/96)*(l' + r + t' + d + fl + fr + ft + fd)
 220                       + (1/192)*(rt + rd + lt + ld + frt + frd + flt + fld)
 221
 222     coefficient 1 0 2 0 =
 223       (1/4)*i' + (1/6)*(f + l' + t')
 224                + (1/12)*(lt + fl + ft)
 225
 226     coefficient 1 0 0 2 =
 227       (1/4)*i' + (1/6)*(f + r + t')
 228                + (1/12)*(rt + fr + ft)
 229
 230     coefficient 1 0 1 1 =
 231       (1/3)*i' + (5/24)*(f + t')
 232                + (1/12)*ft
 233                + (1/24)*(l' + r)
 234                + (1/48)*(lt + rt + fl + fr)
 235
 236     coefficient 1 1 1 0 =
 237       (1/3)*i' + (5/24)*f
 238                + (1/8)*(l' + t')
 239                + (5/96)*(fl + ft)
 240                + (1/48)*(d + r + lt)
 241                + (1/96)*(fd + ld + rt + fr)
 242
 243     coefficient 1 1 0 1 =
 244       (1/3)*i' + (5/24)*f
 245                + (1/8)*(r + t')
 246                + (5/96)*(fr + ft)
 247                + (1/48)*(d + l' + rt)
 248                + (1/96)*(fd + lt + rd + fl)
 249
 250     coefficient 1 2 0 0 =
 251       (1/3)*i' + (5/24)*f
 252                + (7/96)*(l' + r + t' + d)
 253                + (1/32)*(fl + fr + ft + fd)
 254                + (1/96)*(rt + rd + lt + ld)
 255
 256     coefficient 2 0 1 0 =
 257       (3/8)*i' + (7/48)*(f + t' + l')
 258                + (1/48)*(r + d + b + lt + fl + ft)
 259                + (1/96)*(rt + bt + fr + fd + ld + bl)
 260
 261     coefficient 2 0 0 1 =
 262       (3/8)*i' + (7/48)*(f + t' + r)
 263                + (1/48)*(l' + d + b + rt + fr + ft)
 264                + (1/96)*(lt + bt + fl + fd + rd + br)
 265
 266     coefficient 2 1 0 0 =
 267       (3/8)*i' + (1/12)*(t' + r + l' + d)
 268                + (1/64)*(ft + fr + fl + fd)
 269                + (7/48)*f
 270                + (1/48)*b
 271                + (1/96)*(rt + ld + lt + rd)
 272                + (1/192)*(bt + br + bl + bd)
 273
 274     coefficient 3 0 0 0 =
 275       (3/8)*i' + (1/12)*(t' + f + l' + r + d + b)
 276                + (1/96)*(lt + fl + ft + rt + bt + fr)
 277                + (1/96)*(fd + ld + bd + br + rd + bl)
 278
 279
 280
 281 -- | Compute the determinant of the 4x4 matrix,
 282 --
 283 --   [1]
 284 --   [x]
 285 --   [y]
 286 --   [z]
 287 --
 288 --   where [1] = [1, 1, 1, 1],
 289 --         [x] = [x1,x2,x3,x4],
 290 --
 291 --   et cetera.
 292 --
 293 --   The termX nonsense is an attempt to prevent Double overflow.
 294 --   which has been observed to happen with large coordinates.
 295 --
 296 det :: Point -> Point -> Point -> Point -> Double
 297 det p0 p1 p2 p3 =
 298   term5 + term6
 299   where
 300     Point x1 y1 z1 = p0
 301     Point x2 y2 z2 = p1
 302     Point x3 y3 z3 = p2
 303     Point x4 y4 z4 = p3
 304     term1 = ((x2 - x4)*y1 - (x1 - x4)*y2 + (x1 - x2)*y4)*z3
 305     term2 = ((x2 - x3)*y1 - (x1 - x3)*y2 + (x1 - x2)*y3)*z4
 306     term3 = ((x3 - x4)*y2 - (x2 - x4)*y3 + (x2 - x3)*y4)*z1
 307     term4 = ((x3 - x4)*y1 - (x1 - x4)*y3 + (x1 - x3)*y4)*z2
 308     term5 = term1 - term2
 309     term6 = term3 - term4
 310
 311
 312 -- | Computed using the formula from Lai & Schumaker, Definition 15.4,
 313 --   page 436.
 314 {-# INLINE volume #-}
 315 volume :: Tetrahedron -> Double
 316 volume (Tetrahedron _ v0' v1' v2' v3' _) =
 317   (1/6)*(det v0' v1' v2' v3')
 318
 319 -- | The barycentric coordinates of a point with respect to v0.
 320 {-# INLINE b0 #-}
 321 b0 :: Tetrahedron -> (RealFunction Point)
 322 b0 t point = (volume inner_tetrahedron) / (precomputed_volume t)
 323              where
 324                inner_tetrahedron = t { v0 = point }
 325
 326
 327 -- | The barycentric coordinates of a point with respect to v1.
 328 {-# INLINE b1 #-}
 329 b1 :: Tetrahedron -> (RealFunction Point)
 330 b1 t point = (volume inner_tetrahedron) / (precomputed_volume t)
 331              where
 332                inner_tetrahedron = t { v1 = point }
 333
 334
 335 -- | The barycentric coordinates of a point with respect to v2.
 336 {-# INLINE b2 #-}
 337 b2 :: Tetrahedron -> (RealFunction Point)
 338 b2 t point = (volume inner_tetrahedron) / (precomputed_volume t)
 339              where
 340                inner_tetrahedron = t { v2 = point }
 341
 342
 343 -- | The barycentric coordinates of a point with respect to v3.
 344 {-# INLINE b3 #-}
 345 b3 :: Tetrahedron -> (RealFunction Point)
 346 b3 t point = (volume inner_tetrahedron) / (precomputed_volume t)
 347              where
 348                inner_tetrahedron = t { v3 = point }
 349
 350
 351
 352
 353 -- Tests
 354
 355
 356 -- | Check the volume of a particular tetrahedron (computed by hand)
 357 --   and whether or not it contains a specific point chosen to be
 358 --   outside of it. Its vertices are in clockwise order, so the volume
 359 --   should be negative.
 360 tetrahedron1_geometry_tests :: Test.Framework.Test
 361 tetrahedron1_geometry_tests =
 362   testGroup "tetrahedron1 geometry"
 363               [ testCase "volume1" volume1,
 364                 testCase "doesn't contain point1" doesnt_contain_point1]
 365   where
 366     p0 = Point 0 (-0.5) 0
 367     p1 = Point 0 0.5 0
 368     p2 = Point 2 0 0
 369     p3 = Point 1 0 1
 370     t = Tetrahedron { v0 = p0,
 371                       v1 = p1,
 372                       v2 = p2,
 373                       v3 = p3,
 374                       function_values = empty_values,
 375                       precomputed_volume = 0 }
 376
 377     volume1 :: Assertion
 378     volume1 =
 379       assertEqual "volume is correct" True (vol ~= (-1/3))
 380       where
 381         vol = volume t
 382
 383     doesnt_contain_point1 :: Assertion
 384     doesnt_contain_point1 =
 385       assertEqual "doesn't contain an exterior point" False contained
 386       where
 387         exterior_point = Point 5 2 (-9.0212)
 388         contained = contains_point t exterior_point
 389
 390
 391 -- | Check the volume of a particular tetrahedron (computed by hand)
 392 --   and whether or not it contains a specific point chosen to be
 393 --   inside of it. Its vertices are in counter-clockwise order, so the
 394 --   volume should be positive.
 395 tetrahedron2_geometry_tests :: Test.Framework.Test
 396 tetrahedron2_geometry_tests =
 397   testGroup "tetrahedron2 geometry"
 398               [ testCase "volume1" volume1,
 399                 testCase "contains point1" contains_point1]
 400   where
 401     p0 = Point 0 (-0.5) 0
 402     p1 = Point 2 0 0
 403     p2 = Point 0 0.5 0
 404     p3 = Point 1 0 1
 405     t = Tetrahedron { v0 = p0,
 406                       v1 = p1,
 407                       v2 = p2,
 408                       v3 = p3,
 409                       function_values = empty_values,
 410                       precomputed_volume = 0 }
 411
 412     volume1 :: Assertion
 413     volume1 = assertEqual "volume1 is correct" True (vol ~= (1/3))
 414       where
 415         vol = volume t
 416
 417     contains_point1 :: Assertion
 418     contains_point1 = assertEqual "contains an inner point" True contained
 419         where
 420           inner_point = Point 1 0 0.5
 421           contained = contains_point t inner_point
 422
 423
 424 -- | Ensure that tetrahedra do not contain a particular point chosen to
 425 --   be outside of them.
 426 containment_tests :: Test.Framework.Test
 427 containment_tests =
 428   testGroup "containment tests"
 429               [ testCase "doesn't contain point2" doesnt_contain_point2,
 430                 testCase "doesn't contain point3" doesnt_contain_point3,
 431                 testCase "doesn't contain point4" doesnt_contain_point4,
 432                 testCase "doesn't contain point5" doesnt_contain_point5]
 433   where
 434     p2 = Point 0.5 0.5 1
 435     p3 = Point 0.5 0.5 0.5
 436     exterior_point = Point 0 0 0
 437
 438     doesnt_contain_point2 :: Assertion
 439     doesnt_contain_point2 =
 440       assertEqual "doesn't contain an exterior point" False contained
 441       where
 442         p0 = Point 0 1 1
 443         p1 = Point 1 1 1
 444         t = Tetrahedron { v0 = p0,
 445                           v1 = p1,
 446                           v2 = p2,
 447                           v3 = p3,
 448                           function_values = empty_values,
 449                           precomputed_volume = 0 }
 450         contained = contains_point t exterior_point
 451
 452
 453     doesnt_contain_point3 :: Assertion
 454     doesnt_contain_point3 =
 455       assertEqual "doesn't contain an exterior point" False contained
 456       where
 457         p0 = Point 1 1 1
 458         p1 = Point 1 0 1
 459         t = Tetrahedron { v0 = p0,
 460                           v1 = p1,
 461                           v2 = p2,
 462                           v3 = p3,
 463                           function_values = empty_values,
 464                           precomputed_volume = 0 }
 465         contained = contains_point t exterior_point
 466
 467
 468     doesnt_contain_point4 :: Assertion
 469     doesnt_contain_point4 =
 470       assertEqual "doesn't contain an exterior point" False contained
 471       where
 472         p0 = Point 1 0 1
 473         p1 = Point 0 0 1
 474         t = Tetrahedron { v0 = p0,
 475                           v1 = p1,
 476                           v2 = p2,
 477                           v3 = p3,
 478                           function_values = empty_values,
 479                           precomputed_volume = 0 }
 480         contained = contains_point t exterior_point
 481
 482
 483     doesnt_contain_point5 :: Assertion
 484     doesnt_contain_point5 =
 485       assertEqual "doesn't contain an exterior point" False contained
 486       where
 487         p0 = Point 0 0 1
 488         p1 = Point 0 1 1
 489         t = Tetrahedron { v0 = p0,
 490                           v1 = p1,
 491                           v2 = p2,
 492                           v3 = p3,
 493                           function_values = empty_values,
 494                           precomputed_volume = 0 }
 495         contained = contains_point t exterior_point
 496
 497
 498 -- | The barycentric coordinate of v0 with respect to itself should
 499 --   be one.
 500 prop_b0_v0_always_unity :: Tetrahedron -> Property
 501 prop_b0_v0_always_unity t =
 502     (volume t) > 0 ==> (b0 t) (v0 t) ~= 1.0
 503
 504 -- | The barycentric coordinate of v1 with respect to v0 should
 505 --   be zero.
 506 prop_b0_v1_always_zero :: Tetrahedron -> Property
 507 prop_b0_v1_always_zero t =
 508     (volume t) > 0 ==> (b0 t) (v1 t) ~= 0
 509
 510 -- | The barycentric coordinate of v2 with respect to v0 should
 511 --   be zero.
 512 prop_b0_v2_always_zero :: Tetrahedron -> Property
 513 prop_b0_v2_always_zero t =
 514     (volume t) > 0 ==> (b0 t) (v2 t) ~= 0
 515
 516 -- | The barycentric coordinate of v3 with respect to v0 should
 517 --   be zero.
 518 prop_b0_v3_always_zero :: Tetrahedron -> Property
 519 prop_b0_v3_always_zero t =
 520     (volume t) > 0 ==> (b0 t) (v3 t) ~= 0
 521
 522 -- | The barycentric coordinate of v1 with respect to itself should
 523 --   be one.
 524 prop_b1_v1_always_unity :: Tetrahedron -> Property
 525 prop_b1_v1_always_unity t =
 526     (volume t) > 0 ==> (b1 t) (v1 t) ~= 1.0
 527
 528 -- | The barycentric coordinate of v0 with respect to v1 should
 529 --   be zero.
 530 prop_b1_v0_always_zero :: Tetrahedron -> Property
 531 prop_b1_v0_always_zero t =
 532     (volume t) > 0 ==> (b1 t) (v0 t) ~= 0
 533
 534 -- | The barycentric coordinate of v2 with respect to v1 should
 535 --   be zero.
 536 prop_b1_v2_always_zero :: Tetrahedron -> Property
 537 prop_b1_v2_always_zero t =
 538     (volume t) > 0 ==> (b1 t) (v2 t) ~= 0
 539
 540 -- | The barycentric coordinate of v3 with respect to v1 should
 541 --   be zero.
 542 prop_b1_v3_always_zero :: Tetrahedron -> Property
 543 prop_b1_v3_always_zero t =
 544     (volume t) > 0 ==> (b1 t) (v3 t) ~= 0
 545
 546 -- | The barycentric coordinate of v2 with respect to itself should
 547 --   be one.
 548 prop_b2_v2_always_unity :: Tetrahedron -> Property
 549 prop_b2_v2_always_unity t =
 550     (volume t) > 0 ==> (b2 t) (v2 t) ~= 1.0
 551
 552 -- | The barycentric coordinate of v0 with respect to v2 should
 553 --   be zero.
 554 prop_b2_v0_always_zero :: Tetrahedron -> Property
 555 prop_b2_v0_always_zero t =
 556     (volume t) > 0 ==> (b2 t) (v0 t) ~= 0
 557
 558 -- | The barycentric coordinate of v1 with respect to v2 should
 559 --   be zero.
 560 prop_b2_v1_always_zero :: Tetrahedron -> Property
 561 prop_b2_v1_always_zero t =
 562     (volume t) > 0 ==> (b2 t) (v1 t) ~= 0
 563
 564 -- | The barycentric coordinate of v3 with respect to v2 should
 565 --   be zero.
 566 prop_b2_v3_always_zero :: Tetrahedron -> Property
 567 prop_b2_v3_always_zero t =
 568     (volume t) > 0 ==> (b2 t) (v3 t) ~= 0
 569
 570 -- | The barycentric coordinate of v3 with respect to itself should
 571 --   be one.
 572 prop_b3_v3_always_unity :: Tetrahedron -> Property
 573 prop_b3_v3_always_unity t =
 574     (volume t) > 0 ==> (b3 t) (v3 t) ~= 1.0
 575
 576 -- | The barycentric coordinate of v0 with respect to v3 should
 577 --   be zero.
 578 prop_b3_v0_always_zero :: Tetrahedron -> Property
 579 prop_b3_v0_always_zero t =
 580     (volume t) > 0 ==> (b3 t) (v0 t) ~= 0
 581
 582 -- | The barycentric coordinate of v1 with respect to v3 should
 583 --   be zero.
 584 prop_b3_v1_always_zero :: Tetrahedron -> Property
 585 prop_b3_v1_always_zero t =
 586     (volume t) > 0 ==> (b3 t) (v1 t) ~= 0
 587
 588 -- | The barycentric coordinate of v2 with respect to v3 should
 589 --   be zero.
 590 prop_b3_v2_always_zero :: Tetrahedron -> Property
 591 prop_b3_v2_always_zero t =
 592     (volume t) > 0 ==> (b3 t) (v2 t) ~= 0
 593
 594
 595 prop_swapping_vertices_doesnt_affect_coefficients1 :: Tetrahedron -> Bool
 596 prop_swapping_vertices_doesnt_affect_coefficients1 t =
 597       c t 0 0 1 2 == c t' 0 0 1 2
 598         where
 599           t' = t { v0 = (v1 t), v1 = (v0 t) }
 600
 601 prop_swapping_vertices_doesnt_affect_coefficients2 :: Tetrahedron -> Bool
 602 prop_swapping_vertices_doesnt_affect_coefficients2 t =
 603       c t 0 1 1 1 == c t' 0 1 1 1
 604         where
 605           t' = t { v2 = (v3 t), v3 = (v2 t) }
 606
 607 prop_swapping_vertices_doesnt_affect_coefficients3 :: Tetrahedron -> Bool
 608 prop_swapping_vertices_doesnt_affect_coefficients3 t =
 609       c t 2 1 0 0 == c t' 2 1 0 0
 610         where
 611           t' = t { v2 = (v3 t), v3 = (v2 t) }
 612
 613 prop_swapping_vertices_doesnt_affect_coefficients4 :: Tetrahedron -> Bool
 614 prop_swapping_vertices_doesnt_affect_coefficients4 t =
 615       c t 2 0 0 1 == c t' 2 0 0 1
 616         where
 617           t' = t { v0 = (v3 t), v3 = (v0 t) }
 618
 619
 620
 621
 622 tetrahedron_tests :: Test.Framework.Test
 623 tetrahedron_tests =
 624     testGroup "Tetrahedron Tests" [
 625       tetrahedron1_geometry_tests,
 626       tetrahedron2_geometry_tests,
 627       containment_tests ]
 628
 629
 630
 631 p78_24_properties :: Test.Framework.Test
 632 p78_24_properties =
 633   testGroup "p. 78, Section (2.4) Properties" [
 634     testProperty "c3000 identity" prop_c3000_identity,
 635     testProperty "c2100 identity" prop_c2100_identity,
 636     testProperty "c1110 identity" prop_c1110_identity]
 637   where
 638     -- | Returns the domain point of t with indices i,j,k,l.
 639     domain_point :: Tetrahedron -> Int -> Int -> Int -> Int -> Point
 640     domain_point t i j k l =
 641       weighted_sum `scale` (1/3)
 642       where
 643         v0' = (v0 t) `scale` (fromIntegral i)
 644         v1' = (v1 t) `scale` (fromIntegral j)
 645         v2' = (v2 t) `scale` (fromIntegral k)
 646         v3' = (v3 t) `scale` (fromIntegral l)
 647         weighted_sum = v0' + v1' + v2' + v3'
 648
 649
 650     -- | Used for convenience in the next few tests.
 651     p :: Tetrahedron -> Int -> Int -> Int -> Int -> Double
 652     p t i j k l = (polynomial t) (domain_point t i j k l)
 653
 654
 655     -- | Given in Sorokina and Zeilfelder, p. 78.
 656     prop_c3000_identity :: Tetrahedron -> Property
 657     prop_c3000_identity t =
 658       (volume t) > 0 ==>
 659         c t 3 0 0 0 ~= p t 3 0 0 0
 660
 661     -- | Given in Sorokina and Zeilfelder, p. 78.
 662     prop_c2100_identity :: Tetrahedron -> Property
 663     prop_c2100_identity t =
 664       (volume t) > 0 ==>
 665         c t 2 1 0 0 ~= (term1 - term2 + term3 - term4)
 666         where
 667           term1 = (1/3)*(p t 0 3 0 0)
 668           term2 = (5/6)*(p t 3 0 0 0)
 669           term3 = 3*(p t 2 1 0 0)
 670           term4 = (3/2)*(p t 1 2 0 0)
 671
 672     -- | Given in Sorokina and Zeilfelder, p. 78.
 673     prop_c1110_identity :: Tetrahedron -> Property
 674     prop_c1110_identity t =
 675       (volume t) > 0 ==>
 676         c t 1 1 1 0 ~= (term1 + term2 - term3 - term4)
 677         where
 678           term1 = (1/3)*((p t 3 0 0 0) + (p t 0 3 0 0) + (p t 0 0 3 0))
 679           term2 = (9/2)*(p t 1 1 1 0)
 680           term3 = (3/4)*((p t 2 1 0 0) + (p t 1 2 0 0) + (p t 2 0 1 0))
 681           term4 = (3/4)*((p t 1 0 2 0) + (p t 0 2 1 0) + (p t 0 1 2 0))
 682
 683
 684
 685 tetrahedron_properties :: Test.Framework.Test
 686 tetrahedron_properties =
 687   testGroup "Tetrahedron Properties" [
 688     p78_24_properties,
 689     testProperty "b0_v0_always_unity" prop_b0_v0_always_unity,
 690     testProperty "b0_v1_always_zero" prop_b0_v1_always_zero,
 691     testProperty "b0_v2_always_zero" prop_b0_v2_always_zero,
 692     testProperty "b0_v3_always_zero" prop_b0_v3_always_zero,
 693     testProperty "b1_v1_always_unity" prop_b1_v1_always_unity,
 694     testProperty "b1_v0_always_zero" prop_b1_v0_always_zero,
 695     testProperty "b1_v2_always_zero" prop_b1_v2_always_zero,
 696     testProperty "b1_v3_always_zero" prop_b1_v3_always_zero,
 697     testProperty "b2_v2_always_unity" prop_b2_v2_always_unity,
 698     testProperty "b2_v0_always_zero" prop_b2_v0_always_zero,
 699     testProperty "b2_v1_always_zero" prop_b2_v1_always_zero,
 700     testProperty "b2_v3_always_zero" prop_b2_v3_always_zero,
 701     testProperty "b3_v3_always_unity" prop_b3_v3_always_unity,
 702     testProperty "b3_v0_always_zero" prop_b3_v0_always_zero,
 703     testProperty "b3_v1_always_zero" prop_b3_v1_always_zero,
 704     testProperty "b3_v2_always_zero" prop_b3_v2_always_zero,
 705     testProperty "swapping_vertices_doesnt_affect_coefficients1" $
 706       prop_swapping_vertices_doesnt_affect_coefficients1,
 707     testProperty "swapping_vertices_doesnt_affect_coefficients2" $
 708       prop_swapping_vertices_doesnt_affect_coefficients2,
 709     testProperty "swapping_vertices_doesnt_affect_coefficients3" $
 710       prop_swapping_vertices_doesnt_affect_coefficients3,
 711     testProperty "swapping_vertices_doesnt_affect_coefficients4" $
 712       prop_swapping_vertices_doesnt_affect_coefficients4 ]