src/Cube.hs

   1 -- The "tetrahedron" function pattern matches on the integers zero
   2 -- through twenty-three, but doesn't handle the "otherwise" case, for
   3 -- performance reasons.
   4 {-# OPTIONS_GHC -Wno-incomplete-patterns #-}
   5
   6 module Cube (
   7   Cube(..),
   8   cube_properties,
   9   find_containing_tetrahedron,
  10   tetrahedra,
  11   tetrahedron )
  12 where
  13
  14 import Data.Maybe ( fromJust )
  15 import qualified Data.Vector as V (
  16   Vector,
  17   findIndex,
  18   map,
  19   minimum,
  20   singleton,
  21   snoc,
  22   unsafeIndex)
  23 import Prelude(
  24   Bool,
  25   Double,
  26   Int,
  27   Eq( (==) ),
  28   Fractional( (/) ),
  29   Maybe,
  30   Num( (+), (-), (*) ),
  31   Ord( (>=), (<=) ),
  32   Show( show ),
  33   ($),
  34   (.),
  35   (&&),
  36   (++),
  37   abs,
  38   all,
  39   and,
  40   fromIntegral,
  41   head,
  42   map,
  43   otherwise,
  44   return,
  45   tail )
  46 import Test.Tasty ( TestTree, testGroup )
  47 import Test.Tasty.QuickCheck (
  48   Arbitrary( arbitrary ),
  49   Gen,
  50   Positive( Positive ),
  51   choose,
  52   testProperty )
  53 import Cardinal (
  54   Cardinal(F, B, L, R, D, T, FL, FR, FD, FT,
  55            BL, BR, BD, BT, LD, LT, RD, RT, I),
  56   ccwx,
  57   ccwy,
  58   ccwz,
  59   cwx,
  60   cwy,
  61   cwz )
  62 import Comparisons ( (~=), (~~=) )
  63 import qualified Face ( Face(..), center )
  64 import FunctionValues ( FunctionValues, eval, rotate )
  65 import Misc ( all_equal, disjoint )
  66 import Point ( Point( Point ), dot )
  67 import Tetrahedron (
  68   Tetrahedron(Tetrahedron, function_values, v0, v1, v2, v3),
  69   barycenter,
  70   c,
  71   volume )
  72
  73 data Cube = Cube { i  :: !Int,
  74                    j  :: !Int,
  75                    k  :: !Int,
  76                    fv :: !FunctionValues,
  77                    tetrahedra_volume :: !Double }
  78             deriving (Eq)
  79
  80
  81 instance Arbitrary Cube where
  82     arbitrary = do
  83       i' <- choose (coordmin, coordmax)
  84       j' <- choose (coordmin, coordmax)
  85       k' <- choose (coordmin, coordmax)
  86       fv' <- arbitrary :: Gen FunctionValues
  87       (Positive tet_vol) <- arbitrary :: Gen (Positive Double)
  88       return (Cube i' j' k' fv' tet_vol)
  89       where
  90         -- The idea here is that, when cubed in the volume formula,
  91         -- these numbers don't overflow 64 bits. This number is not
  92         -- magic in any other sense than that it does not cause test
  93         -- failures, while 2^23 does.
  94         coordmax = 4194304 :: Int -- 2^22
  95         coordmin = -coordmax
  96
  97
  98 instance Show Cube where
  99     show cube =
 100         "Cube_" ++ subscript ++ "\n" ++
 101         " Center: " ++ (show (center cube)) ++ "\n" ++
 102         " xmin: " ++ (show (xmin cube)) ++ "\n" ++
 103         " xmax: " ++ (show (xmax cube)) ++ "\n" ++
 104         " ymin: " ++ (show (ymin cube)) ++ "\n" ++
 105         " ymax: " ++ (show (ymax cube)) ++ "\n" ++
 106         " zmin: " ++ (show (zmin cube)) ++ "\n" ++
 107         " zmax: " ++ (show (zmax cube)) ++ "\n"
 108         where
 109           subscript =
 110               (show (i cube)) ++ "," ++ (show (j cube)) ++ "," ++ (show (k cube))
 111
 112
 113 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
 114 --   p. 76.
 115 xmin :: Cube -> Double
 116 xmin cube = (i' - 1/2)
 117     where
 118       i' = fromIntegral (i cube) :: Double
 119
 120 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
 121 --   p. 76.
 122 xmax :: Cube -> Double
 123 xmax cube = (i' + 1/2)
 124     where
 125       i' = fromIntegral (i cube) :: Double
 126
 127 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
 128 --   p. 76.
 129 ymin :: Cube -> Double
 130 ymin cube = (j' - 1/2)
 131     where
 132       j' = fromIntegral (j cube) :: Double
 133
 134 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
 135 --   p. 76.
 136 ymax :: Cube -> Double
 137 ymax cube = (j' + 1/2)
 138     where
 139       j' = fromIntegral (j cube) :: Double
 140
 141 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
 142 --   p. 76.
 143 zmin :: Cube -> Double
 144 zmin cube = (k' - 1/2)
 145     where
 146       k' = fromIntegral (k cube) :: Double
 147
 148 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
 149 --   p. 76.
 150 zmax :: Cube -> Double
 151 zmax cube = (k' + 1/2)
 152     where
 153       k' = fromIntegral (k cube) :: Double
 154
 155
 156 -- | The center of Cube_ijk coincides with v_ijk at
 157 --   (i, j, k). See Sorokina and Zeilfelder, p. 76.
 158 center :: Cube -> Point
 159 center cube =
 160   Point x y z
 161   where
 162     x = fromIntegral (i cube) :: Double
 163     y = fromIntegral (j cube) :: Double
 164     z = fromIntegral (k cube) :: Double
 165
 166
 167 -- Face stuff.
 168
 169 -- | The top (in the direction of z) face of the cube.
 170 top_face :: Cube -> Face.Face
 171 top_face cube = Face.Face v0' v1' v2' v3'
 172     where
 173       delta = (1/2) :: Double
 174       cc  = center cube
 175       v0' = cc + ( Point delta (-delta) delta )
 176       v1' = cc + ( Point delta delta delta )
 177       v2' = cc + ( Point (-delta) delta delta )
 178       v3' = cc + ( Point (-delta) (-delta) delta )
 179
 180
 181
 182 -- | The back (in the direction of x) face of the cube.
 183 back_face :: Cube -> Face.Face
 184 back_face cube = Face.Face v0' v1' v2' v3'
 185     where
 186       delta = (1/2) :: Double
 187       cc  = center cube
 188       v0' = cc + ( Point delta (-delta) (-delta) )
 189       v1' = cc + ( Point delta delta (-delta) )
 190       v2' = cc + ( Point delta delta delta )
 191       v3' = cc + ( Point delta (-delta) delta )
 192
 193
 194 -- The bottom face (in the direction of -z) of the cube.
 195 down_face :: Cube -> Face.Face
 196 down_face cube = Face.Face v0' v1' v2' v3'
 197     where
 198       delta = (1/2) :: Double
 199       cc  = center cube
 200       v0' = cc + ( Point (-delta) (-delta) (-delta) )
 201       v1' = cc + ( Point (-delta) delta (-delta) )
 202       v2' = cc + ( Point delta delta (-delta) )
 203       v3' = cc + ( Point delta (-delta) (-delta) )
 204
 205
 206
 207 -- | The front (in the direction of -x) face of the cube.
 208 front_face :: Cube -> Face.Face
 209 front_face cube = Face.Face v0' v1' v2' v3'
 210     where
 211       delta = (1/2) :: Double
 212       cc  = center cube
 213       v0' = cc + ( Point (-delta) (-delta) delta )
 214       v1' = cc + ( Point (-delta) delta delta )
 215       v2' = cc + ( Point (-delta) delta (-delta) )
 216       v3' = cc + ( Point (-delta) (-delta) (-delta) )
 217
 218 -- | The left (in the direction of -y) face of the cube.
 219 left_face :: Cube -> Face.Face
 220 left_face cube = Face.Face v0' v1' v2' v3'
 221     where
 222       delta = (1/2) :: Double
 223       cc  = center cube
 224       v0' = cc + ( Point delta (-delta) delta )
 225       v1' = cc + ( Point (-delta) (-delta) delta )
 226       v2' = cc + ( Point (-delta) (-delta) (-delta) )
 227       v3' = cc + ( Point delta (-delta) (-delta) )
 228
 229
 230 -- | The right (in the direction of y) face of the cube.
 231 right_face :: Cube -> Face.Face
 232 right_face cube = Face.Face v0' v1' v2' v3'
 233     where
 234       delta = (1/2) :: Double
 235       cc  = center cube
 236       v0' = cc + ( Point (-delta) delta delta)
 237       v1' = cc + ( Point delta  delta delta )
 238       v2' = cc + ( Point delta delta (-delta) )
 239       v3' = cc + ( Point (-delta) delta (-delta) )
 240
 241
 242 tetrahedron :: Cube -> Int -> Tetrahedron
 243
 244 tetrahedron cube 0 =
 245     Tetrahedron (fv cube) v0' v1' v2' v3' vol
 246     where
 247       v0' = center cube
 248       ff  = front_face cube
 249       v1' = Face.center ff
 250       v2' = Face.v0 ff
 251       v3' = Face.v1 ff
 252       vol = tetrahedra_volume cube
 253
 254 tetrahedron cube 1 =
 255     Tetrahedron fv' v0' v1' v2' v3' vol
 256     where
 257       v0' = center cube
 258       ff  = front_face cube
 259       v1' = Face.center ff
 260       v2' = Face.v1 ff
 261       v3' = Face.v2 ff
 262       fv' = rotate ccwx (fv cube)
 263       vol = tetrahedra_volume cube
 264
 265 tetrahedron cube 2 =
 266     Tetrahedron fv' v0' v1' v2' v3' vol
 267     where
 268       v0' = center cube
 269       ff  = front_face cube
 270       v1' = Face.center ff
 271       v2' = Face.v2 ff
 272       v3' = Face.v3 ff
 273       fv' = rotate ccwx $ rotate ccwx $ fv cube
 274       vol = tetrahedra_volume cube
 275
 276 tetrahedron cube 3 =
 277     Tetrahedron fv' v0' v1' v2' v3' vol
 278     where
 279       v0' = center cube
 280       ff  = front_face cube
 281       v1' = Face.center ff
 282       v2' = Face.v3 ff
 283       v3' = Face.v0 ff
 284       fv' = rotate cwx (fv cube)
 285       vol = tetrahedra_volume cube
 286
 287 tetrahedron cube 4 =
 288     Tetrahedron fv' v0' v1' v2' v3' vol
 289     where
 290       v0' = center cube
 291       tf  = top_face cube
 292       v1' = Face.center tf
 293       v2' = Face.v0 tf
 294       v3' = Face.v1 tf
 295       fv' = rotate cwy (fv cube)
 296       vol = tetrahedra_volume cube
 297
 298 tetrahedron cube 5 =
 299     Tetrahedron fv' v0' v1' v2' v3' vol
 300     where
 301       v0' = center cube
 302       tf  = top_face cube
 303       v1' = Face.center tf
 304       v2' = Face.v1 tf
 305       v3' = Face.v2 tf
 306       fv' = rotate cwy $ rotate cwz $ fv cube
 307       vol = tetrahedra_volume cube
 308
 309 tetrahedron cube 6 =
 310     Tetrahedron fv' v0' v1' v2' v3' vol
 311     where
 312       v0' = center cube
 313       tf  = top_face cube
 314       v1' = Face.center tf
 315       v2' = Face.v2 tf
 316       v3' = Face.v3 tf
 317       fv' = rotate cwy $ rotate cwz
 318                        $ rotate cwz
 319                        $ fv cube
 320       vol = tetrahedra_volume cube
 321
 322 tetrahedron cube 7 =
 323     Tetrahedron fv' v0' v1' v2' v3' vol
 324     where
 325       v0' = center cube
 326       tf  = top_face cube
 327       v1' = Face.center tf
 328       v2' = Face.v3 tf
 329       v3' = Face.v0 tf
 330       fv' = rotate cwy $ rotate ccwz $ fv cube
 331       vol = tetrahedra_volume cube
 332
 333 tetrahedron cube 8 =
 334     Tetrahedron fv' v0' v1' v2' v3' vol
 335     where
 336       v0' = center cube
 337       bf  = back_face cube
 338       v1' = Face.center bf
 339       v2' = Face.v0 bf
 340       v3' = Face.v1 bf
 341       fv' = rotate cwy $ rotate cwy $ fv cube
 342       vol = tetrahedra_volume cube
 343
 344 tetrahedron cube 9 =
 345     Tetrahedron fv' v0' v1' v2' v3' vol
 346     where
 347       v0' = center cube
 348       bf  = back_face cube
 349       v1' = Face.center bf
 350       v2' = Face.v1 bf
 351       v3' = Face.v2 bf
 352       fv' = rotate cwy $ rotate cwy
 353                        $ rotate cwx
 354                        $ fv cube
 355       vol = tetrahedra_volume cube
 356
 357 tetrahedron cube 10 =
 358     Tetrahedron fv' v0' v1' v2' v3' vol
 359     where
 360       v0' = center cube
 361       bf  = back_face cube
 362       v1' = Face.center bf
 363       v2' = Face.v2 bf
 364       v3' = Face.v3 bf
 365       fv' = rotate cwy $ rotate cwy
 366                        $ rotate cwx
 367                        $ rotate cwx
 368                        $ fv cube
 369
 370       vol = tetrahedra_volume cube
 371
 372 tetrahedron cube 11 =
 373     Tetrahedron fv' v0' v1' v2' v3' vol
 374     where
 375       v0' = center cube
 376       bf  = back_face cube
 377       v1' = Face.center bf
 378       v2' = Face.v3 bf
 379       v3' = Face.v0 bf
 380       fv' = rotate cwy $ rotate cwy
 381                        $ rotate ccwx
 382                        $ fv cube
 383       vol = tetrahedra_volume cube
 384
 385 tetrahedron cube 12 =
 386     Tetrahedron fv' v0' v1' v2' v3' vol
 387     where
 388       v0' = center cube
 389       df  = down_face cube
 390       v1' = Face.center df
 391       v2' = Face.v0 df
 392       v3' = Face.v1 df
 393       fv' = rotate ccwy $ fv cube
 394       vol = tetrahedra_volume cube
 395
 396 tetrahedron cube 13 =
 397     Tetrahedron fv' v0' v1' v2' v3' vol
 398     where
 399       v0' = center cube
 400       df  = down_face cube
 401       v1' = Face.center df
 402       v2' = Face.v1 df
 403       v3' = Face.v2 df
 404       fv' = rotate ccwy $ rotate ccwz $ fv cube
 405       vol = tetrahedra_volume cube
 406
 407 tetrahedron cube 14 =
 408     Tetrahedron fv' v0' v1' v2' v3' vol
 409     where
 410       v0' = center cube
 411       df  = down_face cube
 412       v1' = Face.center df
 413       v2' = Face.v2 df
 414       v3' = Face.v3 df
 415       fv' = rotate ccwy $ rotate ccwz
 416                         $ rotate ccwz
 417                         $ fv cube
 418       vol = tetrahedra_volume cube
 419
 420 tetrahedron cube 15 =
 421     Tetrahedron fv' v0' v1' v2' v3' vol
 422     where
 423       v0' = center cube
 424       df  = down_face cube
 425       v1' = Face.center df
 426       v2' = Face.v3 df
 427       v3' = Face.v0 df
 428       fv' = rotate ccwy $ rotate cwz $ fv cube
 429       vol = tetrahedra_volume cube
 430
 431 tetrahedron cube 16 =
 432     Tetrahedron fv' v0' v1' v2' v3' vol
 433     where
 434       v0' = center cube
 435       rf  = right_face cube
 436       v1' = Face.center rf
 437       v2' = Face.v0 rf
 438       v3' = Face.v1 rf
 439       fv' = rotate ccwz $ fv cube
 440       vol = tetrahedra_volume cube
 441
 442 tetrahedron cube 17 =
 443     Tetrahedron fv' v0' v1' v2' v3' vol
 444     where
 445       v0' = center cube
 446       rf  = right_face cube
 447       v1' = Face.center rf
 448       v2' = Face.v1 rf
 449       v3' = Face.v2 rf
 450       fv' = rotate ccwz $ rotate cwy $ fv cube
 451       vol = tetrahedra_volume cube
 452
 453 tetrahedron cube 18 =
 454     Tetrahedron fv' v0' v1' v2' v3' vol
 455     where
 456       v0' = center cube
 457       rf  = right_face cube
 458       v1' = Face.center rf
 459       v2' = Face.v2 rf
 460       v3' = Face.v3 rf
 461       fv' = rotate ccwz $ rotate cwy
 462                         $ rotate cwy
 463                         $ fv cube
 464       vol = tetrahedra_volume cube
 465
 466 tetrahedron cube 19 =
 467     Tetrahedron fv' v0' v1' v2' v3' vol
 468     where
 469       v0' = center cube
 470       rf  = right_face cube
 471       v1' = Face.center rf
 472       v2' = Face.v3 rf
 473       v3' = Face.v0 rf
 474       fv' = rotate ccwz $ rotate ccwy
 475                         $ fv cube
 476       vol = tetrahedra_volume cube
 477
 478 tetrahedron cube 20 =
 479     Tetrahedron fv' v0' v1' v2' v3' vol
 480     where
 481       v0' = center cube
 482       lf  = left_face cube
 483       v1' = Face.center lf
 484       v2' = Face.v0 lf
 485       v3' = Face.v1 lf
 486       fv' = rotate cwz $ fv cube
 487       vol = tetrahedra_volume cube
 488
 489 tetrahedron cube 21 =
 490     Tetrahedron fv' v0' v1' v2' v3' vol
 491     where
 492       v0' = center cube
 493       lf  = left_face cube
 494       v1' = Face.center lf
 495       v2' = Face.v1 lf
 496       v3' = Face.v2 lf
 497       fv' = rotate cwz $ rotate ccwy $ fv cube
 498       vol = tetrahedra_volume cube
 499
 500 tetrahedron cube 22 =
 501     Tetrahedron fv' v0' v1' v2' v3' vol
 502     where
 503       v0' = center cube
 504       lf  = left_face cube
 505       v1' = Face.center lf
 506       v2' = Face.v2 lf
 507       v3' = Face.v3 lf
 508       fv' = rotate cwz $ rotate ccwy
 509                        $ rotate ccwy
 510                        $ fv cube
 511       vol = tetrahedra_volume cube
 512
 513 tetrahedron cube 23 =
 514     Tetrahedron fv' v0' v1' v2' v3' vol
 515     where
 516       v0' = center cube
 517       lf  = left_face cube
 518       v1' = Face.center lf
 519       v2' = Face.v3 lf
 520       v3' = Face.v0 lf
 521       fv' = rotate cwz $ rotate cwy
 522                        $ fv cube
 523       vol = tetrahedra_volume cube
 524
 525
 526 -- Only used in tests, so we don't need the added speed
 527 -- of Data.Vector.
 528 tetrahedra :: Cube -> [Tetrahedron]
 529 tetrahedra cube = [ tetrahedron cube n | n <- [0..23] ]
 530
 531 front_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
 532 front_left_top_tetrahedra  cube =
 533   V.singleton (tetrahedron cube 0) `V.snoc`
 534     (tetrahedron cube 3) `V.snoc`
 535     (tetrahedron cube 6) `V.snoc`
 536     (tetrahedron cube 7) `V.snoc`
 537     (tetrahedron cube 20) `V.snoc`
 538     (tetrahedron cube 21)
 539
 540 front_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
 541 front_left_down_tetrahedra  cube =
 542   V.singleton (tetrahedron cube 0) `V.snoc`
 543     (tetrahedron cube 2) `V.snoc`
 544     (tetrahedron cube 3) `V.snoc`
 545     (tetrahedron cube 12) `V.snoc`
 546     (tetrahedron cube 15) `V.snoc`
 547     (tetrahedron cube 21)
 548
 549 front_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
 550 front_right_top_tetrahedra  cube =
 551   V.singleton (tetrahedron cube 0) `V.snoc`
 552     (tetrahedron cube 1) `V.snoc`
 553     (tetrahedron cube 5) `V.snoc`
 554     (tetrahedron cube 6) `V.snoc`
 555     (tetrahedron cube 16) `V.snoc`
 556     (tetrahedron cube 19)
 557
 558 front_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
 559 front_right_down_tetrahedra  cube =
 560   V.singleton (tetrahedron cube 1) `V.snoc`
 561     (tetrahedron cube 2) `V.snoc`
 562     (tetrahedron cube 12) `V.snoc`
 563     (tetrahedron cube 13) `V.snoc`
 564     (tetrahedron cube 18) `V.snoc`
 565     (tetrahedron cube 19)
 566
 567 back_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
 568 back_left_top_tetrahedra  cube =
 569   V.singleton (tetrahedron cube 0) `V.snoc`
 570     (tetrahedron cube 3) `V.snoc`
 571     (tetrahedron cube 6) `V.snoc`
 572     (tetrahedron cube 7) `V.snoc`
 573     (tetrahedron cube 20) `V.snoc`
 574     (tetrahedron cube 21)
 575
 576 back_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
 577 back_left_down_tetrahedra  cube =
 578   V.singleton (tetrahedron cube 8) `V.snoc`
 579     (tetrahedron cube 11) `V.snoc`
 580     (tetrahedron cube 14) `V.snoc`
 581     (tetrahedron cube 15) `V.snoc`
 582     (tetrahedron cube 22) `V.snoc`
 583     (tetrahedron cube 23)
 584
 585 back_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
 586 back_right_top_tetrahedra  cube =
 587   V.singleton (tetrahedron cube 4) `V.snoc`
 588     (tetrahedron cube 5) `V.snoc`
 589     (tetrahedron cube 9) `V.snoc`
 590     (tetrahedron cube 10) `V.snoc`
 591     (tetrahedron cube 16) `V.snoc`
 592     (tetrahedron cube 17)
 593
 594 back_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
 595 back_right_down_tetrahedra  cube =
 596   V.singleton (tetrahedron cube 8) `V.snoc`
 597     (tetrahedron cube 9) `V.snoc`
 598     (tetrahedron cube 13) `V.snoc`
 599     (tetrahedron cube 14) `V.snoc`
 600     (tetrahedron cube 17) `V.snoc`
 601     (tetrahedron cube 18)
 602
 603 in_top_half :: Cube -> Point -> Bool
 604 in_top_half cube (Point _ _ z) =
 605   distance_from_top <= distance_from_bottom
 606   where
 607     distance_from_top = abs $ (zmax cube) - z
 608     distance_from_bottom = abs $ (zmin cube) - z
 609
 610 in_front_half :: Cube -> Point -> Bool
 611 in_front_half cube (Point x _ _) =
 612     distance_from_front <= distance_from_back
 613   where
 614     distance_from_front = abs $ (xmin cube) - x
 615     distance_from_back = abs $ (xmax cube) - x
 616
 617
 618 in_left_half :: Cube -> Point -> Bool
 619 in_left_half cube (Point _ y _) =
 620     distance_from_left <= distance_from_right
 621   where
 622     distance_from_left = abs $ (ymin cube) - y
 623     distance_from_right = abs $ (ymax cube) - y
 624
 625
 626 -- | Takes a 'Cube', and returns the Tetrahedra belonging to it that
 627 --   contain the given 'Point'. This should be faster than checking
 628 --   every tetrahedron individually, since we determine which half
 629 --   (hemisphere?) of the cube the point lies in three times: once in
 630 --   each dimension. This allows us to eliminate non-candidates
 631 --   quickly.
 632 --
 633 --   This can throw an exception, but the use of 'head' might
 634 --   save us some unnecessary computations.
 635 --
 636 {-# INLINE find_containing_tetrahedron #-}
 637 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
 638 find_containing_tetrahedron cube p =
 639   candidates `V.unsafeIndex` (fromJust lucky_idx)
 640   where
 641     front_half = in_front_half cube p
 642     top_half = in_top_half cube p
 643     left_half = in_left_half cube p
 644
 645     candidates :: V.Vector Tetrahedron
 646     candidates
 647       | front_half =
 648         if left_half then
 649           if top_half then
 650             front_left_top_tetrahedra cube
 651           else
 652             front_left_down_tetrahedra cube
 653         else
 654           if top_half then
 655             front_right_top_tetrahedra cube
 656           else
 657             front_right_down_tetrahedra cube
 658
 659       | otherwise = -- back half
 660         if left_half then
 661           if top_half then
 662             back_left_top_tetrahedra cube
 663           else
 664             back_left_down_tetrahedra cube
 665         else
 666           if top_half then
 667             back_right_top_tetrahedra cube
 668           else
 669             back_right_down_tetrahedra cube
 670
 671     -- Use the dot product instead of Euclidean distance here to save
 672     -- a sqrt(). So, "distances" below really means "distances
 673     -- squared."
 674     distances :: V.Vector Double
 675     distances = V.map ((dot p) . barycenter) candidates
 676
 677     shortest_distance :: Double
 678     shortest_distance = V.minimum distances
 679
 680     -- Compute the index of the tetrahedron with the center closest to
 681     -- p. This is a bad algorithm, but don't change it! If you make it
 682     -- smarter by finding the index of shortest_distance in distances
 683     -- (this should give the same answer and avoids recomputing the
 684     -- dot product), the program gets slower. Seriously!
 685     lucky_idx :: Maybe Int
 686     lucky_idx = V.findIndex
 687                   (\t -> (barycenter t) `dot` p == shortest_distance)
 688                   candidates
 689
 690
 691
 692
 693
 694
 695 -- * Tests
 696
 697 prop_opposite_octant_tetrahedra_disjoint1 :: Cube -> Bool
 698 prop_opposite_octant_tetrahedra_disjoint1 cube =
 699   disjoint (front_left_top_tetrahedra cube) (front_right_down_tetrahedra cube)
 700
 701 prop_opposite_octant_tetrahedra_disjoint2 :: Cube -> Bool
 702 prop_opposite_octant_tetrahedra_disjoint2 cube =
 703   disjoint (back_left_top_tetrahedra cube) (back_right_down_tetrahedra cube)
 704
 705 prop_opposite_octant_tetrahedra_disjoint3 :: Cube -> Bool
 706 prop_opposite_octant_tetrahedra_disjoint3 cube =
 707   disjoint (front_left_top_tetrahedra cube) (back_right_top_tetrahedra cube)
 708
 709 prop_opposite_octant_tetrahedra_disjoint4 :: Cube -> Bool
 710 prop_opposite_octant_tetrahedra_disjoint4 cube =
 711   disjoint (front_left_down_tetrahedra cube) (back_right_down_tetrahedra cube)
 712
 713 prop_opposite_octant_tetrahedra_disjoint5 :: Cube -> Bool
 714 prop_opposite_octant_tetrahedra_disjoint5 cube =
 715   disjoint (front_left_top_tetrahedra cube) (back_left_down_tetrahedra cube)
 716
 717 prop_opposite_octant_tetrahedra_disjoint6 :: Cube -> Bool
 718 prop_opposite_octant_tetrahedra_disjoint6 cube =
 719   disjoint (front_right_top_tetrahedra cube) (back_right_down_tetrahedra cube)
 720
 721
 722 -- | Since the grid size is necessarily positive, all tetrahedra
 723 --   (which comprise cubes of positive volume) must have positive
 724 --   volume as well.
 725 prop_all_volumes_positive :: Cube -> Bool
 726 prop_all_volumes_positive cube =
 727     all (>= 0) volumes
 728     where
 729       ts = tetrahedra cube
 730       volumes = map volume ts
 731
 732
 733 -- | In fact, since all of the tetrahedra are identical, we should
 734 --   already know their volumes. There's 24 tetrahedra to a cube, so
 735 --   we'd expect the volume of each one to be 1/24.
 736 prop_all_volumes_exact :: Cube -> Bool
 737 prop_all_volumes_exact cube =
 738     and [volume t ~~= 1/24 | t <- tetrahedra cube]
 739
 740 -- | All tetrahedron should have their v0 located at the center of the
 741 --   cube.
 742 prop_v0_all_equal :: Cube -> Bool
 743 prop_v0_all_equal cube = (v0 t0) == (v0 t1)
 744     where
 745       t0 = head (tetrahedra cube) -- Doesn't matter which two we choose.
 746       t1 = head $ tail (tetrahedra cube)
 747
 748
 749 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Note that the
 750 --   third and fourth indices of c-t3 have been switched. This is
 751 --   because we store the triangles oriented such that their volume is
 752 --   positive. If T and T-tilde share \<v0,v1,v2\> and v3,v3-tilde point
 753 --   in opposite directions, one of them has to have negative volume!
 754 prop_c0120_identity1 :: Cube -> Bool
 755 prop_c0120_identity1 cube =
 756    c t0 0 1 2 0 ~= (c t0 0 0 2 1 + c t3 0 0 1 2) / 2
 757      where
 758        t0 = tetrahedron cube 0
 759        t3 = tetrahedron cube 3
 760
 761
 762 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
 763 --   'prop_c0120_identity1' with tetrahedrons 1 and 2.
 764 prop_c0120_identity2 :: Cube -> Bool
 765 prop_c0120_identity2 cube =
 766    c t1 0 1 2 0 ~= (c t1 0 0 2 1 + c t0 0 0 1 2) / 2
 767      where
 768        t0 = tetrahedron cube 0
 769        t1 = tetrahedron cube 1
 770
 771 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
 772 --   'prop_c0120_identity1' with tetrahedrons 1 and 2.
 773 prop_c0120_identity3 :: Cube -> Bool
 774 prop_c0120_identity3 cube =
 775    c t2 0 1 2 0 ~= (c t2 0 0 2 1 + c t1 0 0 1 2) / 2
 776      where
 777        t1 = tetrahedron cube 1
 778        t2 = tetrahedron cube 2
 779
 780 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
 781 --   'prop_c0120_identity1' with tetrahedrons 2 and 3.
 782 prop_c0120_identity4 :: Cube -> Bool
 783 prop_c0120_identity4 cube =
 784    c t3 0 1 2 0 ~= (c t3 0 0 2 1 + c t2 0 0 1 2) / 2
 785      where
 786        t2 = tetrahedron cube 2
 787        t3 = tetrahedron cube 3
 788
 789
 790 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
 791 --   'prop_c0120_identity1' with tetrahedrons 4 and 5.
 792 prop_c0120_identity5 :: Cube -> Bool
 793 prop_c0120_identity5 cube =
 794     c t5 0 1 2 0 ~= (c t5 0 0 2 1 + c t4 0 0 1 2) / 2
 795     where
 796       t4 = tetrahedron cube 4
 797       t5 = tetrahedron cube 5
 798
 799 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
 800 --   'prop_c0120_identity1' with tetrahedrons 5 and 6.
 801 prop_c0120_identity6 :: Cube -> Bool
 802 prop_c0120_identity6 cube =
 803    c t6 0 1 2 0 ~= (c t6 0 0 2 1 + c t5 0 0 1 2) / 2
 804      where
 805        t5 = tetrahedron cube 5
 806        t6 = tetrahedron cube 6
 807
 808
 809 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
 810 --   'prop_c0120_identity1' with tetrahedrons 6 and 7.
 811 prop_c0120_identity7 :: Cube -> Bool
 812 prop_c0120_identity7 cube =
 813    c t7 0 1 2 0 ~= (c t7 0 0 2 1 + c t6 0 0 1 2) / 2
 814      where
 815        t6 = tetrahedron cube 6
 816        t7 = tetrahedron cube 7
 817
 818
 819 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
 820 --   'prop_c0120_identity1'.
 821 prop_c0210_identity1 :: Cube -> Bool
 822 prop_c0210_identity1 cube =
 823     c t0 0 2 1 0 ~= (c t0 0 1 1 1 + c t3 0 1 1 1) / 2
 824       where
 825         t0 = tetrahedron cube 0
 826         t3 = tetrahedron cube 3
 827
 828
 829 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
 830 --   'prop_c0120_identity1'.
 831 prop_c0300_identity1 :: Cube -> Bool
 832 prop_c0300_identity1 cube =
 833     c t0 0 3 0 0 ~= (c t0 0 2 0 1 + c t3 0 2 1 0) / 2
 834       where
 835         t0 = tetrahedron cube 0
 836         t3 = tetrahedron cube 3
 837
 838
 839 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
 840 --   'prop_c0120_identity1'.
 841 prop_c1110_identity :: Cube -> Bool
 842 prop_c1110_identity cube =
 843     c t0 1 1 1 0 ~= (c t0 1 0 1 1 + c t3 1 0 1 1) / 2
 844       where
 845         t0 = tetrahedron cube 0
 846         t3 = tetrahedron cube 3
 847
 848
 849 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
 850 --   'prop_c0120_identity1'.
 851 prop_c1200_identity1 :: Cube -> Bool
 852 prop_c1200_identity1 cube =
 853     c t0 1 2 0 0 ~= (c t0 1 1 0 1 + c t3 1 1 1 0) / 2
 854       where
 855         t0 = tetrahedron cube 0
 856         t3 = tetrahedron cube 3
 857
 858
 859 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
 860 --   'prop_c0120_identity1'.
 861 prop_c2100_identity1 :: Cube -> Bool
 862 prop_c2100_identity1 cube =
 863     c t0 2 1 0 0 ~= (c t0 2 0 0 1 + c t3 2 0 1 0) / 2
 864       where
 865         t0 = tetrahedron cube 0
 866         t3 = tetrahedron cube 3
 867
 868
 869
 870 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). Note that the
 871 --   third and fourth indices of c-t3 have been switched. This is
 872 --   because we store the triangles oriented such that their volume is
 873 --   positive. If T and T-tilde share \<v0,v1,v2\> and v3,v3-tilde
 874 --   point in opposite directions, one of them has to have negative
 875 --   volume!
 876 prop_c0102_identity1 :: Cube -> Bool
 877 prop_c0102_identity1 cube =
 878     c t0 0 1 0 2 ~= (c t0 0 0 1 2 + c t1 0 0 2 1) / 2
 879       where
 880         t0 = tetrahedron cube 0
 881         t1 = tetrahedron cube 1
 882
 883
 884 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
 885 --   'prop_c0102_identity1'.
 886 prop_c0201_identity1 :: Cube -> Bool
 887 prop_c0201_identity1 cube =
 888     c t0 0 2 0 1 ~= (c t0 0 1 1 1 + c t1 0 1 1 1) / 2
 889       where
 890         t0 = tetrahedron cube 0
 891         t1 = tetrahedron cube 1
 892
 893
 894 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
 895 --   'prop_c0102_identity1'.
 896 prop_c0300_identity2 :: Cube -> Bool
 897 prop_c0300_identity2 cube =
 898     c t0 0 3 0 0 ~= (c t0 0 2 1 0 + c t1 0 2 0 1) / 2
 899       where
 900         t0 = tetrahedron cube 0
 901         t1 = tetrahedron cube 1
 902
 903
 904 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
 905 --   'prop_c0102_identity1'.
 906 prop_c1101_identity :: Cube -> Bool
 907 prop_c1101_identity cube =
 908     c t0 1 1 0 1 ~= (c t0 1 0 1 1 + c t1 1 0 1 1) / 2
 909       where
 910         t0 = tetrahedron cube 0
 911         t1 = tetrahedron cube 1
 912
 913
 914 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
 915 --   'prop_c0102_identity1'.
 916 prop_c1200_identity2 :: Cube -> Bool
 917 prop_c1200_identity2 cube =
 918     c t0 1 2 0 0 ~= (c t0 1 1 1 0 + c t1 1 1 0 1) / 2
 919       where
 920         t0 = tetrahedron cube 0
 921         t1 = tetrahedron cube 1
 922
 923
 924 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
 925 --   'prop_c0102_identity1'.
 926 prop_c2100_identity2 :: Cube -> Bool
 927 prop_c2100_identity2 cube =
 928     c t0 2 1 0 0 ~= (c t0 2 0 1 0 + c t1 2 0 0 1) / 2
 929       where
 930         t0 = tetrahedron cube 0
 931         t1 = tetrahedron cube 1
 932
 933
 934 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). The third and
 935 --   fourth indices of c-t6 have been switched. This is because we
 936 --   store the triangles oriented such that their volume is
 937 --   positive. If T and T-tilde share \<v0,v1,v2\> and v3,v3-tilde
 938 --   point in opposite directions, one of them has to have negative
 939 --   volume!
 940 prop_c3000_identity :: Cube -> Bool
 941 prop_c3000_identity cube =
 942     c t0 3 0 0 0 ~= c t0 2 1 0 0 + c t6 2 1 0 0
 943                     - ((c t0 2 0 1 0 + c t0 2 0 0 1)/ 2)
 944       where
 945         t0 = tetrahedron cube 0
 946         t6 = tetrahedron cube 6
 947
 948
 949 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
 950 --   'prop_c3000_identity'.
 951 prop_c2010_identity :: Cube -> Bool
 952 prop_c2010_identity cube =
 953     c t0 2 0 1 0 ~= c t0 1 1 1 0 + c t6 1 1 0 1
 954                     - ((c t0 1 0 2 0 + c t0 1 0 1 1)/ 2)
 955       where
 956         t0 = tetrahedron cube 0
 957         t6 = tetrahedron cube 6
 958
 959
 960 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
 961 --   'prop_c3000_identity'.
 962 prop_c2001_identity :: Cube -> Bool
 963 prop_c2001_identity cube =
 964     c t0 2 0 0 1 ~= c t0 1 1 0 1 + c t6 1 1 1 0
 965                     - ((c t0 1 0 0 2 + c t0 1 0 1 1)/ 2)
 966       where
 967         t0 = tetrahedron cube 0
 968         t6 = tetrahedron cube 6
 969
 970
 971 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
 972 --   'prop_c3000_identity'.
 973 prop_c1020_identity :: Cube -> Bool
 974 prop_c1020_identity cube =
 975     c t0 1 0 2 0 ~= c t0 0 1 2 0 + c t6 0 1 0 2
 976                     - ((c t0 0 0 3 0 + c t0 0 0 2 1)/ 2)
 977       where
 978         t0 = tetrahedron cube 0
 979         t6 = tetrahedron cube 6
 980
 981
 982 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
 983 --   'prop_c3000_identity'.
 984 prop_c1002_identity :: Cube -> Bool
 985 prop_c1002_identity cube =
 986     c t0 1 0 0 2 ~= c t0 0 1 0 2 + c t6 0 1 2 0
 987                     - ((c t0 0 0 0 3 + c t0 0 0 1 2)/ 2)
 988       where
 989         t0 = tetrahedron cube 0
 990         t6 = tetrahedron cube 6
 991
 992
 993 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
 994 --   'prop_c3000_identity'.
 995 prop_c1011_identity :: Cube -> Bool
 996 prop_c1011_identity cube =
 997     c t0 1 0 1 1 ~= c t0 0 1 1 1 + c t6 0 1 1 1 -
 998                     ((c t0 0 0 1 2 + c t0 0 0 2 1)/ 2)
 999       where
1000         t0 = tetrahedron cube 0
1001         t6 = tetrahedron cube 6
1002
1003
1004 -- | The function values at the interior should be the same for all
1005 --   tetrahedra.
1006 prop_interior_values_all_identical :: Cube -> Bool
1007 prop_interior_values_all_identical cube =
1008     all_equal [ eval (function_values tet) I | tet <- tetrahedra cube ]
1009
1010
1011 -- | We know what (c t6 2 1 0 0) should be from Sorokina and Zeilfelder, p. 87.
1012 --   This test checks the rotation works as expected.
1013 prop_c_tilde_2100_rotation_correct :: Cube -> Bool
1014 prop_c_tilde_2100_rotation_correct cube =
1015     expr1 ~= expr2
1016     where
1017       t0 = tetrahedron cube 0
1018       t6 = tetrahedron cube 6
1019
1020       -- What gets computed for c2100 of t6.
1021       expr1 = eval (function_values t6) $
1022           (3/8)*I +
1023           (1/12)*(T + R + L + D) +
1024           (1/64)*(FT + FR + FL + FD) +
1025           (7/48)*F +
1026           (1/48)*B +
1027           (1/96)*(RT + LD + LT + RD) +
1028           (1/192)*(BT + BR + BL + BD)
1029
1030       -- What should be computed for c2100 of t6.
1031       expr2 = eval (function_values t0) $
1032               (3/8)*I +
1033               (1/12)*(F + R + L + B) +
1034               (1/64)*(FT + RT + LT + BT) +
1035               (7/48)*T +
1036               (1/48)*D +
1037               (1/96)*(FR + FL + BR + BL) +
1038               (1/192)*(FD + RD + LD + BD)
1039
1040
1041 -- | We know what (c t6 2 1 0 0) should be from Sorokina and
1042 --   Zeilfelder, p. 87.  This test checks the actual value based on
1043 --   the FunctionValues of the cube.
1044 --
1045 --   If 'prop_c_tilde_2100_rotation_correct' passes, then this test is
1046 --   even meaningful!
1047 prop_c_tilde_2100_correct :: Cube -> Bool
1048 prop_c_tilde_2100_correct cube =
1049     c t6 2 1 0 0 ~= expected
1050     where
1051       t0 = tetrahedron cube 0
1052       t6 = tetrahedron cube 6
1053       fvs = function_values t0
1054       expected = eval fvs $
1055                   (3/8)*I +
1056                   (1/12)*(F + R + L + B) +
1057                   (1/64)*(FT + RT + LT + BT) +
1058                   (7/48)*T +
1059                   (1/48)*D +
1060                   (1/96)*(FR + FL + BR + BL) +
1061                   (1/192)*(FD + RD + LD + BD)
1062
1063
1064 -- Tests to check that the correct edges are incidental.
1065 prop_t0_shares_edge_with_t1 :: Cube -> Bool
1066 prop_t0_shares_edge_with_t1 cube =
1067     (v1 t0) == (v1 t1) && (v3 t0) == (v2 t1)
1068       where
1069         t0 = tetrahedron cube 0
1070         t1 = tetrahedron cube 1
1071
1072 prop_t0_shares_edge_with_t3 :: Cube -> Bool
1073 prop_t0_shares_edge_with_t3 cube =
1074     (v1 t0) == (v1 t3) && (v2 t0) == (v3 t3)
1075       where
1076         t0 = tetrahedron cube 0
1077         t3 = tetrahedron cube 3
1078
1079 prop_t0_shares_edge_with_t6 :: Cube -> Bool
1080 prop_t0_shares_edge_with_t6 cube =
1081     (v2 t0) == (v3 t6) && (v3 t0) == (v2 t6)
1082       where
1083         t0 = tetrahedron cube 0
1084         t6 = tetrahedron cube 6
1085
1086 prop_t1_shares_edge_with_t2 :: Cube -> Bool
1087 prop_t1_shares_edge_with_t2 cube =
1088     (v1 t1) == (v1 t2) && (v3 t1) == (v2 t2)
1089       where
1090         t1 = tetrahedron cube 1
1091         t2 = tetrahedron cube 2
1092
1093 prop_t1_shares_edge_with_t19 :: Cube -> Bool
1094 prop_t1_shares_edge_with_t19 cube =
1095     (v2 t1) == (v3 t19) && (v3 t1) == (v2 t19)
1096       where
1097         t1 = tetrahedron cube 1
1098         t19 = tetrahedron cube 19
1099
1100 prop_t2_shares_edge_with_t3 :: Cube -> Bool
1101 prop_t2_shares_edge_with_t3 cube =
1102     (v1 t1) == (v1 t2) && (v3 t1) == (v2 t2)
1103       where
1104         t1 = tetrahedron cube 1
1105         t2 = tetrahedron cube 2
1106
1107 prop_t2_shares_edge_with_t12 :: Cube -> Bool
1108 prop_t2_shares_edge_with_t12 cube =
1109     (v2 t2) == (v3 t12) && (v3 t2) == (v2 t12)
1110       where
1111         t2 = tetrahedron cube 2
1112         t12 = tetrahedron cube 12
1113
1114 prop_t3_shares_edge_with_t21 :: Cube -> Bool
1115 prop_t3_shares_edge_with_t21 cube =
1116     (v2 t3) == (v3 t21) && (v3 t3) == (v2 t21)
1117       where
1118         t3 = tetrahedron cube 3
1119         t21 = tetrahedron cube 21
1120
1121 prop_t4_shares_edge_with_t5 :: Cube -> Bool
1122 prop_t4_shares_edge_with_t5 cube =
1123     (v1 t4) == (v1 t5) && (v3 t4) == (v2 t5)
1124       where
1125         t4 = tetrahedron cube 4
1126         t5 = tetrahedron cube 5
1127
1128 prop_t4_shares_edge_with_t7 :: Cube -> Bool
1129 prop_t4_shares_edge_with_t7 cube =
1130     (v1 t4) == (v1 t7) && (v2 t4) == (v3 t7)
1131       where
1132         t4 = tetrahedron cube 4
1133         t7 = tetrahedron cube 7
1134
1135 prop_t4_shares_edge_with_t10 :: Cube -> Bool
1136 prop_t4_shares_edge_with_t10 cube =
1137     (v2 t4) == (v3 t10) && (v3 t4) == (v2 t10)
1138       where
1139         t4 = tetrahedron cube 4
1140         t10 = tetrahedron cube 10
1141
1142 prop_t5_shares_edge_with_t6 :: Cube -> Bool
1143 prop_t5_shares_edge_with_t6 cube =
1144     (v1 t5) == (v1 t6) && (v3 t5) == (v2 t6)
1145       where
1146         t5 = tetrahedron cube 5
1147         t6 = tetrahedron cube 6
1148
1149 prop_t5_shares_edge_with_t16 :: Cube -> Bool
1150 prop_t5_shares_edge_with_t16 cube =
1151     (v2 t5) == (v3 t16) && (v3 t5) == (v2 t16)
1152       where
1153         t5 = tetrahedron cube 5
1154         t16 = tetrahedron cube 16
1155
1156 prop_t6_shares_edge_with_t7 :: Cube -> Bool
1157 prop_t6_shares_edge_with_t7 cube =
1158     (v1 t6) == (v1 t7) && (v3 t6) == (v2 t7)
1159       where
1160         t6 = tetrahedron cube 6
1161         t7 = tetrahedron cube 7
1162
1163 prop_t7_shares_edge_with_t20 :: Cube -> Bool
1164 prop_t7_shares_edge_with_t20 cube =
1165     (v2 t7) == (v3 t20) && (v2 t7) == (v3 t20)
1166       where
1167         t7 = tetrahedron cube 7
1168         t20 = tetrahedron cube 20
1169
1170
1171 p79_26_properties :: TestTree
1172 p79_26_properties =
1173     testGroup "p. 79, Section (2.6) properties" [
1174       testProperty "c0120 identity1" prop_c0120_identity1,
1175       testProperty "c0120 identity2" prop_c0120_identity2,
1176       testProperty "c0120 identity3" prop_c0120_identity3,
1177       testProperty "c0120 identity4" prop_c0120_identity4,
1178       testProperty "c0120 identity5" prop_c0120_identity5,
1179       testProperty "c0120 identity6" prop_c0120_identity6,
1180       testProperty "c0120 identity7" prop_c0120_identity7,
1181       testProperty "c0210 identity1" prop_c0210_identity1,
1182       testProperty "c0300 identity1" prop_c0300_identity1,
1183       testProperty "c1110 identity" prop_c1110_identity,
1184       testProperty "c1200 identity1" prop_c1200_identity1,
1185       testProperty "c2100 identity1" prop_c2100_identity1]
1186
1187 p79_27_properties :: TestTree
1188 p79_27_properties =
1189     testGroup "p. 79, Section (2.7) properties" [
1190       testProperty "c0102 identity1" prop_c0102_identity1,
1191       testProperty "c0201 identity1" prop_c0201_identity1,
1192       testProperty "c0300 identity2" prop_c0300_identity2,
1193       testProperty "c1101 identity" prop_c1101_identity,
1194       testProperty "c1200 identity2" prop_c1200_identity2,
1195       testProperty "c2100 identity2" prop_c2100_identity2 ]
1196
1197
1198 p79_28_properties :: TestTree
1199 p79_28_properties =
1200   testGroup "p. 79, Section (2.8) properties" [
1201     testProperty "c3000 identity" prop_c3000_identity,
1202     testProperty "c2010 identity" prop_c2010_identity,
1203     testProperty "c2001 identity" prop_c2001_identity,
1204     testProperty "c1020 identity" prop_c1020_identity,
1205     testProperty "c1002 identity" prop_c1002_identity,
1206     testProperty "c1011 identity" prop_c1011_identity ]
1207
1208
1209 edge_incidence_tests :: TestTree
1210 edge_incidence_tests =
1211     testGroup "Edge incidence tests" [
1212       testProperty "t0 shares edge with t6" prop_t0_shares_edge_with_t6,
1213       testProperty "t0 shares edge with t1" prop_t0_shares_edge_with_t1,
1214       testProperty "t0 shares edge with t3" prop_t0_shares_edge_with_t3,
1215       testProperty "t1 shares edge with t2" prop_t1_shares_edge_with_t2,
1216       testProperty "t1 shares edge with t19" prop_t1_shares_edge_with_t19,
1217       testProperty "t2 shares edge with t3" prop_t2_shares_edge_with_t3,
1218       testProperty "t2 shares edge with t12" prop_t2_shares_edge_with_t12,
1219       testProperty "t3 shares edge with t21" prop_t3_shares_edge_with_t21,
1220       testProperty "t4 shares edge with t5" prop_t4_shares_edge_with_t5,
1221       testProperty "t4 shares edge with t7" prop_t4_shares_edge_with_t7,
1222       testProperty "t4 shares edge with t10" prop_t4_shares_edge_with_t10,
1223       testProperty "t5 shares edge with t6" prop_t5_shares_edge_with_t6,
1224       testProperty "t5 shares edge with t16" prop_t5_shares_edge_with_t16,
1225       testProperty "t6 shares edge with t7" prop_t6_shares_edge_with_t7,
1226       testProperty "t7 shares edge with t20" prop_t7_shares_edge_with_t20 ]
1227
1228 cube_properties :: TestTree
1229 cube_properties =
1230   testGroup "Cube properties" [
1231     p79_26_properties,
1232     p79_27_properties,
1233     p79_28_properties,
1234     edge_incidence_tests,
1235     testProperty "opposite octant tetrahedra are disjoint (1)"
1236       prop_opposite_octant_tetrahedra_disjoint1,
1237     testProperty "opposite octant tetrahedra are disjoint (2)"
1238       prop_opposite_octant_tetrahedra_disjoint2,
1239     testProperty "opposite octant tetrahedra are disjoint (3)"
1240      prop_opposite_octant_tetrahedra_disjoint3,
1241     testProperty "opposite octant tetrahedra are disjoint (4)"
1242       prop_opposite_octant_tetrahedra_disjoint4,
1243     testProperty "opposite octant tetrahedra are disjoint (5)"
1244       prop_opposite_octant_tetrahedra_disjoint5,
1245     testProperty "opposite octant tetrahedra are disjoint (6)"
1246       prop_opposite_octant_tetrahedra_disjoint6,
1247     testProperty "all volumes positive" prop_all_volumes_positive,
1248     testProperty "all volumes exact" prop_all_volumes_exact,
1249     testProperty "v0 all equal" prop_v0_all_equal,
1250     testProperty "interior values all identical"
1251       prop_interior_values_all_identical,
1252     testProperty "c-tilde_2100 rotation correct"
1253       prop_c_tilde_2100_rotation_correct,
1254     testProperty "c-tilde_2100 correct"
1255       prop_c_tilde_2100_correct ]