src/Tetrahedron.hs

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