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 t
 317        | v0' == v1' = 0
 318        | v0' == v2' = 0
 319        | v0' == v3' = 0
 320        | v1' == v2' = 0
 321        | v1' == v3' = 0
 322        | v2' == v3' = 0
 323        | otherwise = (1/6)*(det v0' v1' v2' v3')
 324   where
 325     v0' = v0 t
 326     v1' = v1 t
 327     v2' = v2 t
 328     v3' = v3 t
 329
 330
 331 -- | The barycentric coordinates of a point with respect to v0.
 332 {-# INLINE b0 #-}
 333 b0 :: Tetrahedron -> (RealFunction Point)
 334 b0 t point = (volume inner_tetrahedron) / (precomputed_volume t)
 335              where
 336                inner_tetrahedron = t { v0 = point }
 337
 338
 339 -- | The barycentric coordinates of a point with respect to v1.
 340 {-# INLINE b1 #-}
 341 b1 :: Tetrahedron -> (RealFunction Point)
 342 b1 t point = (volume inner_tetrahedron) / (precomputed_volume t)
 343              where
 344                inner_tetrahedron = t { v1 = point }
 345
 346
 347 -- | The barycentric coordinates of a point with respect to v2.
 348 {-# INLINE b2 #-}
 349 b2 :: Tetrahedron -> (RealFunction Point)
 350 b2 t point = (volume inner_tetrahedron) / (precomputed_volume t)
 351              where
 352                inner_tetrahedron = t { v2 = point }
 353
 354
 355 -- | The barycentric coordinates of a point with respect to v3.
 356 {-# INLINE b3 #-}
 357 b3 :: Tetrahedron -> (RealFunction Point)
 358 b3 t point = (volume inner_tetrahedron) / (precomputed_volume t)
 359              where
 360                inner_tetrahedron = t { v3 = point }
 361
 362
 363
 364
 365 -- Tests
 366
 367
 368 -- | Check the volume of a particular tetrahedron (computed by hand)
 369 --   and whether or not it contains a specific point chosen to be
 370 --   outside of it. Its vertices are in clockwise order, so the volume
 371 --   should be negative.
 372 tetrahedron1_geometry_tests :: Test.Framework.Test
 373 tetrahedron1_geometry_tests =
 374   testGroup "tetrahedron1 geometry"
 375               [ testCase "volume1" volume1,
 376                 testCase "doesn't contain point1" doesnt_contain_point1]
 377   where
 378     p0 = Point 0 (-0.5) 0
 379     p1 = Point 0 0.5 0
 380     p2 = Point 2 0 0
 381     p3 = Point 1 0 1
 382     t = Tetrahedron { v0 = p0,
 383                       v1 = p1,
 384                       v2 = p2,
 385                       v3 = p3,
 386                       function_values = empty_values,
 387                       precomputed_volume = 0 }
 388
 389     volume1 :: Assertion
 390     volume1 =
 391       assertEqual "volume is correct" True (vol ~= (-1/3))
 392       where
 393         vol = volume t
 394
 395     doesnt_contain_point1 :: Assertion
 396     doesnt_contain_point1 =
 397       assertEqual "doesn't contain an exterior point" False contained
 398       where
 399         exterior_point = Point 5 2 (-9.0212)
 400         contained = contains_point t exterior_point
 401
 402
 403 -- | Check the volume of a particular tetrahedron (computed by hand)
 404 --   and whether or not it contains a specific point chosen to be
 405 --   inside of it. Its vertices are in counter-clockwise order, so the
 406 --   volume should be positive.
 407 tetrahedron2_geometry_tests :: Test.Framework.Test
 408 tetrahedron2_geometry_tests =
 409   testGroup "tetrahedron2 geometry"
 410               [ testCase "volume1" volume1,
 411                 testCase "contains point1" contains_point1]
 412   where
 413     p0 = Point 0 (-0.5) 0
 414     p1 = Point 2 0 0
 415     p2 = Point 0 0.5 0
 416     p3 = Point 1 0 1
 417     t = Tetrahedron { v0 = p0,
 418                       v1 = p1,
 419                       v2 = p2,
 420                       v3 = p3,
 421                       function_values = empty_values,
 422                       precomputed_volume = 0 }
 423
 424     volume1 :: Assertion
 425     volume1 = assertEqual "volume1 is correct" True (vol ~= (1/3))
 426       where
 427         vol = volume t
 428
 429     contains_point1 :: Assertion
 430     contains_point1 = assertEqual "contains an inner point" True contained
 431         where
 432           inner_point = Point 1 0 0.5
 433           contained = contains_point t inner_point
 434
 435
 436 -- | Ensure that tetrahedra do not contain a particular point chosen to
 437 --   be outside of them.
 438 containment_tests :: Test.Framework.Test
 439 containment_tests =
 440   testGroup "containment tests"
 441               [ testCase "doesn't contain point2" doesnt_contain_point2,
 442                 testCase "doesn't contain point3" doesnt_contain_point3,
 443                 testCase "doesn't contain point4" doesnt_contain_point4,
 444                 testCase "doesn't contain point5" doesnt_contain_point5]
 445   where
 446     p2 = Point 0.5 0.5 1
 447     p3 = Point 0.5 0.5 0.5
 448     exterior_point = Point 0 0 0
 449
 450     doesnt_contain_point2 :: Assertion
 451     doesnt_contain_point2 =
 452       assertEqual "doesn't contain an exterior point" False contained
 453       where
 454         p0 = Point 0 1 1
 455         p1 = Point 1 1 1
 456         t = Tetrahedron { v0 = p0,
 457                           v1 = p1,
 458                           v2 = p2,
 459                           v3 = p3,
 460                           function_values = empty_values,
 461                           precomputed_volume = 0 }
 462         contained = contains_point t exterior_point
 463
 464
 465     doesnt_contain_point3 :: Assertion
 466     doesnt_contain_point3 =
 467       assertEqual "doesn't contain an exterior point" False contained
 468       where
 469         p0 = Point 1 1 1
 470         p1 = Point 1 0 1
 471         t = Tetrahedron { v0 = p0,
 472                           v1 = p1,
 473                           v2 = p2,
 474                           v3 = p3,
 475                           function_values = empty_values,
 476                           precomputed_volume = 0 }
 477         contained = contains_point t exterior_point
 478
 479
 480     doesnt_contain_point4 :: Assertion
 481     doesnt_contain_point4 =
 482       assertEqual "doesn't contain an exterior point" False contained
 483       where
 484         p0 = Point 1 0 1
 485         p1 = Point 0 0 1
 486         t = Tetrahedron { v0 = p0,
 487                           v1 = p1,
 488                           v2 = p2,
 489                           v3 = p3,
 490                           function_values = empty_values,
 491                           precomputed_volume = 0 }
 492         contained = contains_point t exterior_point
 493
 494
 495     doesnt_contain_point5 :: Assertion
 496     doesnt_contain_point5 =
 497       assertEqual "doesn't contain an exterior point" False contained
 498       where
 499         p0 = Point 0 0 1
 500         p1 = Point 0 1 1
 501         t = Tetrahedron { v0 = p0,
 502                           v1 = p1,
 503                           v2 = p2,
 504                           v3 = p3,
 505                           function_values = empty_values,
 506                           precomputed_volume = 0 }
 507         contained = contains_point t exterior_point
 508
 509
 510 -- | The barycentric coordinate of v0 with respect to itself should
 511 --   be one.
 512 prop_b0_v0_always_unity :: Tetrahedron -> Property
 513 prop_b0_v0_always_unity t =
 514     (volume t) > 0 ==> (b0 t) (v0 t) ~= 1.0
 515
 516 -- | The barycentric coordinate of v1 with respect to v0 should
 517 --   be zero.
 518 prop_b0_v1_always_zero :: Tetrahedron -> Property
 519 prop_b0_v1_always_zero t =
 520     (volume t) > 0 ==> (b0 t) (v1 t) ~= 0
 521
 522 -- | The barycentric coordinate of v2 with respect to v0 should
 523 --   be zero.
 524 prop_b0_v2_always_zero :: Tetrahedron -> Property
 525 prop_b0_v2_always_zero t =
 526     (volume t) > 0 ==> (b0 t) (v2 t) ~= 0
 527
 528 -- | The barycentric coordinate of v3 with respect to v0 should
 529 --   be zero.
 530 prop_b0_v3_always_zero :: Tetrahedron -> Property
 531 prop_b0_v3_always_zero t =
 532     (volume t) > 0 ==> (b0 t) (v3 t) ~= 0
 533
 534 -- | The barycentric coordinate of v1 with respect to itself should
 535 --   be one.
 536 prop_b1_v1_always_unity :: Tetrahedron -> Property
 537 prop_b1_v1_always_unity t =
 538     (volume t) > 0 ==> (b1 t) (v1 t) ~= 1.0
 539
 540 -- | The barycentric coordinate of v0 with respect to v1 should
 541 --   be zero.
 542 prop_b1_v0_always_zero :: Tetrahedron -> Property
 543 prop_b1_v0_always_zero t =
 544     (volume t) > 0 ==> (b1 t) (v0 t) ~= 0
 545
 546 -- | The barycentric coordinate of v2 with respect to v1 should
 547 --   be zero.
 548 prop_b1_v2_always_zero :: Tetrahedron -> Property
 549 prop_b1_v2_always_zero t =
 550     (volume t) > 0 ==> (b1 t) (v2 t) ~= 0
 551
 552 -- | The barycentric coordinate of v3 with respect to v1 should
 553 --   be zero.
 554 prop_b1_v3_always_zero :: Tetrahedron -> Property
 555 prop_b1_v3_always_zero t =
 556     (volume t) > 0 ==> (b1 t) (v3 t) ~= 0
 557
 558 -- | The barycentric coordinate of v2 with respect to itself should
 559 --   be one.
 560 prop_b2_v2_always_unity :: Tetrahedron -> Property
 561 prop_b2_v2_always_unity t =
 562     (volume t) > 0 ==> (b2 t) (v2 t) ~= 1.0
 563
 564 -- | The barycentric coordinate of v0 with respect to v2 should
 565 --   be zero.
 566 prop_b2_v0_always_zero :: Tetrahedron -> Property
 567 prop_b2_v0_always_zero t =
 568     (volume t) > 0 ==> (b2 t) (v0 t) ~= 0
 569
 570 -- | The barycentric coordinate of v1 with respect to v2 should
 571 --   be zero.
 572 prop_b2_v1_always_zero :: Tetrahedron -> Property
 573 prop_b2_v1_always_zero t =
 574     (volume t) > 0 ==> (b2 t) (v1 t) ~= 0
 575
 576 -- | The barycentric coordinate of v3 with respect to v2 should
 577 --   be zero.
 578 prop_b2_v3_always_zero :: Tetrahedron -> Property
 579 prop_b2_v3_always_zero t =
 580     (volume t) > 0 ==> (b2 t) (v3 t) ~= 0
 581
 582 -- | The barycentric coordinate of v3 with respect to itself should
 583 --   be one.
 584 prop_b3_v3_always_unity :: Tetrahedron -> Property
 585 prop_b3_v3_always_unity t =
 586     (volume t) > 0 ==> (b3 t) (v3 t) ~= 1.0
 587
 588 -- | The barycentric coordinate of v0 with respect to v3 should
 589 --   be zero.
 590 prop_b3_v0_always_zero :: Tetrahedron -> Property
 591 prop_b3_v0_always_zero t =
 592     (volume t) > 0 ==> (b3 t) (v0 t) ~= 0
 593
 594 -- | The barycentric coordinate of v1 with respect to v3 should
 595 --   be zero.
 596 prop_b3_v1_always_zero :: Tetrahedron -> Property
 597 prop_b3_v1_always_zero t =
 598     (volume t) > 0 ==> (b3 t) (v1 t) ~= 0
 599
 600 -- | The barycentric coordinate of v2 with respect to v3 should
 601 --   be zero.
 602 prop_b3_v2_always_zero :: Tetrahedron -> Property
 603 prop_b3_v2_always_zero t =
 604     (volume t) > 0 ==> (b3 t) (v2 t) ~= 0
 605
 606
 607 prop_swapping_vertices_doesnt_affect_coefficients1 :: Tetrahedron -> Bool
 608 prop_swapping_vertices_doesnt_affect_coefficients1 t =
 609       c t 0 0 1 2 == c t' 0 0 1 2
 610         where
 611           t' = t { v0 = (v1 t), v1 = (v0 t) }
 612
 613 prop_swapping_vertices_doesnt_affect_coefficients2 :: Tetrahedron -> Bool
 614 prop_swapping_vertices_doesnt_affect_coefficients2 t =
 615       c t 0 1 1 1 == c t' 0 1 1 1
 616         where
 617           t' = t { v2 = (v3 t), v3 = (v2 t) }
 618
 619 prop_swapping_vertices_doesnt_affect_coefficients3 :: Tetrahedron -> Bool
 620 prop_swapping_vertices_doesnt_affect_coefficients3 t =
 621       c t 2 1 0 0 == c t' 2 1 0 0
 622         where
 623           t' = t { v2 = (v3 t), v3 = (v2 t) }
 624
 625 prop_swapping_vertices_doesnt_affect_coefficients4 :: Tetrahedron -> Bool
 626 prop_swapping_vertices_doesnt_affect_coefficients4 t =
 627       c t 2 0 0 1 == c t' 2 0 0 1
 628         where
 629           t' = t { v0 = (v3 t), v3 = (v0 t) }
 630
 631
 632
 633
 634 tetrahedron_tests :: Test.Framework.Test
 635 tetrahedron_tests =
 636     testGroup "Tetrahedron Tests" [
 637       tetrahedron1_geometry_tests,
 638       tetrahedron2_geometry_tests,
 639       containment_tests ]
 640
 641
 642
 643 p78_24_properties :: Test.Framework.Test
 644 p78_24_properties =
 645   testGroup "p. 78, Section (2.4) Properties" [
 646     testProperty "c3000 identity" prop_c3000_identity,
 647     testProperty "c2100 identity" prop_c2100_identity,
 648     testProperty "c1110 identity" prop_c1110_identity]
 649   where
 650     -- | Returns the domain point of t with indices i,j,k,l.
 651     domain_point :: Tetrahedron -> Int -> Int -> Int -> Int -> Point
 652     domain_point t i j k l =
 653       weighted_sum `scale` (1/3)
 654       where
 655         v0' = (v0 t) `scale` (fromIntegral i)
 656         v1' = (v1 t) `scale` (fromIntegral j)
 657         v2' = (v2 t) `scale` (fromIntegral k)
 658         v3' = (v3 t) `scale` (fromIntegral l)
 659         weighted_sum = v0' + v1' + v2' + v3'
 660
 661
 662     -- | Used for convenience in the next few tests.
 663     p :: Tetrahedron -> Int -> Int -> Int -> Int -> Double
 664     p t i j k l = (polynomial t) (domain_point t i j k l)
 665
 666
 667     -- | Given in Sorokina and Zeilfelder, p. 78.
 668     prop_c3000_identity :: Tetrahedron -> Property
 669     prop_c3000_identity t =
 670       (volume t) > 0 ==>
 671         c t 3 0 0 0 ~= p t 3 0 0 0
 672
 673     -- | Given in Sorokina and Zeilfelder, p. 78.
 674     prop_c2100_identity :: Tetrahedron -> Property
 675     prop_c2100_identity t =
 676       (volume t) > 0 ==>
 677         c t 2 1 0 0 ~= (term1 - term2 + term3 - term4)
 678         where
 679           term1 = (1/3)*(p t 0 3 0 0)
 680           term2 = (5/6)*(p t 3 0 0 0)
 681           term3 = 3*(p t 2 1 0 0)
 682           term4 = (3/2)*(p t 1 2 0 0)
 683
 684     -- | Given in Sorokina and Zeilfelder, p. 78.
 685     prop_c1110_identity :: Tetrahedron -> Property
 686     prop_c1110_identity t =
 687       (volume t) > 0 ==>
 688         c t 1 1 1 0 ~= (term1 + term2 - term3 - term4)
 689         where
 690           term1 = (1/3)*((p t 3 0 0 0) + (p t 0 3 0 0) + (p t 0 0 3 0))
 691           term2 = (9/2)*(p t 1 1 1 0)
 692           term3 = (3/4)*((p t 2 1 0 0) + (p t 1 2 0 0) + (p t 2 0 1 0))
 693           term4 = (3/4)*((p t 1 0 2 0) + (p t 0 2 1 0) + (p t 0 1 2 0))
 694
 695
 696
 697 tetrahedron_properties :: Test.Framework.Test
 698 tetrahedron_properties =
 699   testGroup "Tetrahedron Properties" [
 700     p78_24_properties,
 701     testProperty "b0_v0_always_unity" prop_b0_v0_always_unity,
 702     testProperty "b0_v1_always_zero" prop_b0_v1_always_zero,
 703     testProperty "b0_v2_always_zero" prop_b0_v2_always_zero,
 704     testProperty "b0_v3_always_zero" prop_b0_v3_always_zero,
 705     testProperty "b1_v1_always_unity" prop_b1_v1_always_unity,
 706     testProperty "b1_v0_always_zero" prop_b1_v0_always_zero,
 707     testProperty "b1_v2_always_zero" prop_b1_v2_always_zero,
 708     testProperty "b1_v3_always_zero" prop_b1_v3_always_zero,
 709     testProperty "b2_v2_always_unity" prop_b2_v2_always_unity,
 710     testProperty "b2_v0_always_zero" prop_b2_v0_always_zero,
 711     testProperty "b2_v1_always_zero" prop_b2_v1_always_zero,
 712     testProperty "b2_v3_always_zero" prop_b2_v3_always_zero,
 713     testProperty "b3_v3_always_unity" prop_b3_v3_always_unity,
 714     testProperty "b3_v0_always_zero" prop_b3_v0_always_zero,
 715     testProperty "b3_v1_always_zero" prop_b3_v1_always_zero,
 716     testProperty "b3_v2_always_zero" prop_b3_v2_always_zero,
 717     testProperty "swapping_vertices_doesnt_affect_coefficients1" $
 718       prop_swapping_vertices_doesnt_affect_coefficients1,
 719     testProperty "swapping_vertices_doesnt_affect_coefficients2" $
 720       prop_swapping_vertices_doesnt_affect_coefficients2,
 721     testProperty "swapping_vertices_doesnt_affect_coefficients3" $
 722       prop_swapping_vertices_doesnt_affect_coefficients3,
 723     testProperty "swapping_vertices_doesnt_affect_coefficients4" $
 724       prop_swapping_vertices_doesnt_affect_coefficients4 ]