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