src/Cube.hs

   1 module Cube
   2 where
   3
   4 import Data.List ( (\\) )
   5 import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
   6
   7 import Cardinal
   8 import qualified Face (Face(Face, v0, v1, v2, v3))
   9 import FunctionValues
  10 import Point
  11 import Tetrahedron hiding (c)
  12 import ThreeDimensional
  13
  14 data Cube = Cube { h :: Double,
  15                    i :: Int,
  16                    j :: Int,
  17                    k :: Int,
  18                    fv :: FunctionValues }
  19             deriving (Eq)
  20
  21
  22 instance Arbitrary Cube where
  23     arbitrary = do
  24       (Positive h') <- arbitrary :: Gen (Positive Double)
  25       i' <- choose (coordmin, coordmax)
  26       j' <- choose (coordmin, coordmax)
  27       k' <- choose (coordmin, coordmax)
  28       fv' <- arbitrary :: Gen FunctionValues
  29       return (Cube h' i' j' k' fv')
  30         where
  31           coordmin = -268435456 -- -(2^29 / 2)
  32           coordmax = 268435456  -- +(2^29 / 2)
  33
  34
  35 instance Show Cube where
  36     show c =
  37         "Cube_" ++ subscript ++ "\n" ++
  38         " h: " ++ (show (h c)) ++ "\n" ++
  39         " Center: " ++ (show (center c)) ++ "\n" ++
  40         " xmin: " ++ (show (xmin c)) ++ "\n" ++
  41         " xmax: " ++ (show (xmax c)) ++ "\n" ++
  42         " ymin: " ++ (show (ymin c)) ++ "\n" ++
  43         " ymax: " ++ (show (ymax c)) ++ "\n" ++
  44         " zmin: " ++ (show (zmin c)) ++ "\n" ++
  45         " zmax: " ++ (show (zmax c)) ++ "\n" ++
  46         " fv: " ++ (show (Cube.fv c)) ++ "\n"
  47         where
  48           subscript =
  49               (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k c))
  50
  51
  52 -- | Returns an empty 'Cube'.
  53 empty_cube :: Cube
  54 empty_cube = Cube 0 0 0 0 empty_values
  55
  56
  57 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
  58 --   p. 76.
  59 xmin :: Cube -> Double
  60 xmin c = (2*i' - 1)*delta / 2
  61     where
  62       i' = fromIntegral (i c) :: Double
  63       delta = h c
  64
  65 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
  66 --   p. 76.
  67 xmax :: Cube -> Double
  68 xmax c = (2*i' + 1)*delta / 2
  69     where
  70       i' = fromIntegral (i c) :: Double
  71       delta = h c
  72
  73 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
  74 --   p. 76.
  75 ymin :: Cube -> Double
  76 ymin c = (2*j' - 1)*delta / 2
  77     where
  78       j' = fromIntegral (j c) :: Double
  79       delta = h c
  80
  81 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
  82 --   p. 76.
  83 ymax :: Cube -> Double
  84 ymax c = (2*j' + 1)*delta / 2
  85     where
  86       j' = fromIntegral (j c) :: Double
  87       delta = h c
  88
  89 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
  90 --   p. 76.
  91 zmin :: Cube -> Double
  92 zmin c = (2*k' - 1)*delta / 2
  93     where
  94       k' = fromIntegral (k c) :: Double
  95       delta = h c
  96
  97 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
  98 --   p. 76.
  99 zmax :: Cube -> Double
 100 zmax c = (2*k' + 1)*delta / 2
 101     where
 102       k' = fromIntegral (k c) :: Double
 103       delta = h c
 104
 105 instance ThreeDimensional Cube where
 106     -- | The center of Cube_ijk coincides with v_ijk at
 107     --   (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
 108     center c = (x, y, z)
 109            where
 110              delta = h c
 111              i' = fromIntegral (i c) :: Double
 112              j' = fromIntegral (j c) :: Double
 113              k' = fromIntegral (k c) :: Double
 114              x = delta * i'
 115              y = delta * j'
 116              z = delta * k'
 117
 118     -- | It's easy to tell if a point is within a cube; just make sure
 119     --   that it falls on the proper side of each of the cube's faces.
 120     contains_point c (x, y, z)
 121         | x < (xmin c) = False
 122         | x > (xmax c) = False
 123         | y < (ymin c) = False
 124         | y > (ymax c) = False
 125         | z < (zmin c) = False
 126         | z > (zmax c) = False
 127         | otherwise = True
 128
 129
 130
 131 -- Face stuff.
 132
 133 -- | The top (in the direction of z) face of the cube.
 134 top_face :: Cube -> Face.Face
 135 top_face c = Face.Face v0' v1' v2' v3'
 136     where
 137       delta = (1/2)*(h c)
 138       v0' = (center c) + (delta, -delta, delta)
 139       v1' = (center c) + (delta, delta, delta)
 140       v2' = (center c) + (-delta, delta, delta)
 141       v3' = (center c) + (-delta, -delta, delta)
 142
 143
 144
 145 -- | The back (in the direction of x) face of the cube.
 146 back_face :: Cube -> Face.Face
 147 back_face c = Face.Face v0' v1' v2' v3'
 148     where
 149       delta = (1/2)*(h c)
 150       v0' = (center c) + (delta, -delta, -delta)
 151       v1' = (center c) + (delta, delta, -delta)
 152       v2' = (center c) + (delta, delta, delta)
 153       v3' = (center c) + (delta, -delta, delta)
 154
 155
 156 -- The bottom face (in the direction of -z) of the cube.
 157 down_face :: Cube -> Face.Face
 158 down_face c = Face.Face v0' v1' v2' v3'
 159     where
 160       delta = (1/2)*(h c)
 161       v0' = (center c) + (-delta, -delta, -delta)
 162       v1' = (center c) + (-delta, delta, -delta)
 163       v2' = (center c) + (delta, delta, -delta)
 164       v3' = (center c) + (delta, -delta, -delta)
 165
 166
 167
 168 -- | The front (in the direction of -x) face of the cube.
 169 front_face :: Cube -> Face.Face
 170 front_face c = Face.Face v0' v1' v2' v3'
 171     where
 172       delta = (1/2)*(h c)
 173       v0' = (center c) + (-delta, -delta, delta)
 174       v1' = (center c) + (-delta, delta, delta)
 175       v2' = (center c) + (-delta, delta, -delta)
 176       v3' = (center c) + (-delta, -delta, -delta)
 177
 178 -- | The left (in the direction of -y) face of the cube.
 179 left_face :: Cube -> Face.Face
 180 left_face c = Face.Face v0' v1' v2' v3'
 181     where
 182       delta = (1/2)*(h c)
 183       v0' = (center c) + (delta, -delta, delta)
 184       v1' = (center c) + (-delta, -delta, delta)
 185       v2' = (center c) + (-delta, -delta, -delta)
 186       v3' = (center c) + (delta, -delta, -delta)
 187
 188
 189 -- | The right (in the direction of y) face of the cube.
 190 right_face :: Cube -> Face.Face
 191 right_face c = Face.Face v0' v1' v2' v3'
 192     where
 193       delta = (1/2)*(h c)
 194       v0' = (center c) + (-delta, delta, delta)
 195       v1' = (center c) + (delta, delta, delta)
 196       v2' = (center c) + (delta, delta, -delta)
 197       v3' = (center c) + (-delta, delta, -delta)
 198
 199
 200 tetrahedron0 :: Cube -> Tetrahedron
 201 tetrahedron0 c =
 202     Tetrahedron (Cube.fv c) v0' v1' v2' v3' vol
 203     where
 204       v0' = center c
 205       v1' = center (front_face c)
 206       v2' = Face.v0 (front_face c)
 207       v3' = Face.v1 (front_face c)
 208       vol = 0
 209
 210 tetrahedron1 :: Cube -> Tetrahedron
 211 tetrahedron1 c =
 212     Tetrahedron fv' v0' v1' v2' v3' vol
 213     where
 214       v0' = center c
 215       v1' = center (front_face c)
 216       v2' = Face.v1 (front_face c)
 217       v3' = Face.v2 (front_face c)
 218       fv' = rotate ccwx (Cube.fv c)
 219       vol = 0
 220
 221 tetrahedron2 :: Cube -> Tetrahedron
 222 tetrahedron2 c =
 223     Tetrahedron fv' v0' v1' v2' v3' vol
 224     where
 225       v0' = center c
 226       v1' = center (front_face c)
 227       v2' = Face.v2 (front_face c)
 228       v3' = Face.v3 (front_face c)
 229       fv' = rotate ccwx $ rotate ccwx $ Cube.fv c
 230       vol = 0
 231
 232 tetrahedron3 :: Cube -> Tetrahedron
 233 tetrahedron3 c =
 234     Tetrahedron fv' v0' v1' v2' v3' vol
 235     where
 236       v0' = center c
 237       v1' = center (front_face c)
 238       v2' = Face.v3 (front_face c)
 239       v3' = Face.v0 (front_face c)
 240       fv' = rotate cwx (Cube.fv c)
 241       vol = 0
 242
 243 tetrahedron4 :: Cube -> Tetrahedron
 244 tetrahedron4 c =
 245     Tetrahedron fv' v0' v1' v2' v3' vol
 246     where
 247       v0' = center c
 248       v1' = center (top_face c)
 249       v2' = Face.v0 (top_face c)
 250       v3' = Face.v1 (top_face c)
 251       fv' = rotate cwy (Cube.fv c)
 252       vol = 0
 253
 254 tetrahedron5 :: Cube -> Tetrahedron
 255 tetrahedron5 c =
 256     Tetrahedron fv' v0' v1' v2' v3' vol
 257     where
 258       v0' = center c
 259       v1' = center (top_face c)
 260       v2' = Face.v1 (top_face c)
 261       v3' = Face.v2 (top_face c)
 262       fv' = rotate cwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
 263       vol = 0
 264
 265 tetrahedron6 :: Cube -> Tetrahedron
 266 tetrahedron6 c =
 267     Tetrahedron fv' v0' v1' v2' v3' vol
 268     where
 269       v0' = center c
 270       v1' = center (top_face c)
 271       v2' = Face.v2 (top_face c)
 272       v3' = Face.v3 (top_face c)
 273       fv' = rotate cwy $ rotate cwz
 274                        $ rotate cwz
 275                        $ Tetrahedron.fv (tetrahedron0 c)
 276       vol = 0
 277
 278 tetrahedron7 :: Cube -> Tetrahedron
 279 tetrahedron7 c =
 280     Tetrahedron fv' v0' v1' v2' v3' vol
 281     where
 282       v0' = center c
 283       v1' = center (top_face c)
 284       v2' = Face.v3 (top_face c)
 285       v3' = Face.v0 (top_face c)
 286       fv' = rotate cwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
 287       vol = 0
 288
 289 tetrahedron8 :: Cube -> Tetrahedron
 290 tetrahedron8 c =
 291     Tetrahedron fv' v0' v1' v2' v3' vol
 292     where
 293       v0' = center c
 294       v1' = center (back_face c)
 295       v2' = Face.v0 (back_face c)
 296       v3' = Face.v1 (back_face c)
 297       fv' = rotate cwy $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)
 298       vol = 0
 299
 300 tetrahedron9 :: Cube -> Tetrahedron
 301 tetrahedron9 c =
 302     Tetrahedron fv' v0' v1' v2' v3' vol
 303     where
 304       v0' = center c
 305       v1' = center (back_face c)
 306       v2' = Face.v1 (back_face c)
 307       v3' = Face.v2 (back_face c)
 308       fv' = rotate cwy $ rotate cwy
 309                        $ rotate cwx
 310                        $ Tetrahedron.fv (tetrahedron0 c)
 311       vol = 0
 312
 313 tetrahedron10 :: Cube -> Tetrahedron
 314 tetrahedron10 c =
 315     Tetrahedron fv' v0' v1' v2' v3' vol
 316     where
 317       v0' = center c
 318       v1' = center (back_face c)
 319       v2' = Face.v2 (back_face c)
 320       v3' = Face.v3 (back_face c)
 321       fv' = rotate cwy $ rotate cwy
 322                        $ rotate cwx
 323                        $ rotate cwx
 324                        $ Tetrahedron.fv (tetrahedron0 c)
 325
 326       vol = 0
 327
 328 tetrahedron11 :: Cube -> Tetrahedron
 329 tetrahedron11 c =
 330     Tetrahedron fv' v0' v1' v2' v3' vol
 331     where
 332       v0' = center c
 333       v1' = center (back_face c)
 334       v2' = Face.v3 (back_face c)
 335       v3' = Face.v0 (back_face c)
 336       fv' = rotate cwy $ rotate cwy
 337                        $ rotate ccwx
 338                        $ Tetrahedron.fv (tetrahedron0 c)
 339       vol = 0
 340
 341
 342 tetrahedron12 :: Cube -> Tetrahedron
 343 tetrahedron12 c =
 344     Tetrahedron fv' v0' v1' v2' v3' vol
 345     where
 346       v0' = center c
 347       v1' = center (down_face c)
 348       v2' = Face.v0 (down_face c)
 349       v3' = Face.v1 (down_face c)
 350       fv' = rotate ccwy (Tetrahedron.fv (tetrahedron0 c))
 351       vol = 0
 352
 353
 354 tetrahedron13 :: Cube -> Tetrahedron
 355 tetrahedron13 c =
 356     Tetrahedron fv' v0' v1' v2' v3' vol
 357     where
 358       v0' = center c
 359       v1' = center (down_face c)
 360       v2' = Face.v1 (down_face c)
 361       v3' = Face.v2 (down_face c)
 362       fv' = rotate ccwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
 363       vol = 0
 364
 365
 366 tetrahedron14 :: Cube -> Tetrahedron
 367 tetrahedron14 c =
 368     Tetrahedron fv' v0' v1' v2' v3' vol
 369     where
 370       v0' = center c
 371       v1' = center (down_face c)
 372       v2' = Face.v2 (down_face c)
 373       v3' = Face.v3 (down_face c)
 374       fv' = rotate ccwy $ rotate ccwz
 375                         $ rotate ccwz
 376                         $ Tetrahedron.fv (tetrahedron0 c)
 377       vol = 0
 378
 379
 380 tetrahedron15 :: Cube -> Tetrahedron
 381 tetrahedron15 c =
 382     Tetrahedron fv' v0' v1' v2' v3' vol
 383     where
 384       v0' = center c
 385       v1' = center (down_face c)
 386       v2' = Face.v3 (down_face c)
 387       v3' = Face.v0 (down_face c)
 388       fv' = rotate ccwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
 389       vol = 0
 390
 391
 392 tetrahedron16 :: Cube -> Tetrahedron
 393 tetrahedron16 c =
 394     Tetrahedron fv' v0' v1' v2' v3' vol
 395     where
 396       v0' = center c
 397       v1' = center (right_face c)
 398       v2' = Face.v0 (right_face c)
 399       v3' = Face.v1 (right_face c)
 400       fv' = rotate ccwz (Tetrahedron.fv (tetrahedron0 c))
 401       vol = 0
 402
 403
 404 tetrahedron17 :: Cube -> Tetrahedron
 405 tetrahedron17 c =
 406     Tetrahedron fv' v0' v1' v2' v3' vol
 407     where
 408       v0' = center c
 409       v1' = center (right_face c)
 410       v2' = Face.v1 (right_face c)
 411       v3' = Face.v2 (right_face c)
 412       fv' = rotate ccwz $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)
 413       vol = 0
 414
 415
 416 tetrahedron18 :: Cube -> Tetrahedron
 417 tetrahedron18 c =
 418     Tetrahedron fv' v0' v1' v2' v3' vol
 419     where
 420       v0' = center c
 421       v1' = center (right_face c)
 422       v2' = Face.v2 (right_face c)
 423       v3' = Face.v3 (right_face c)
 424       fv' = rotate ccwz $ rotate cwy
 425                         $ rotate cwy
 426                         $ Tetrahedron.fv (tetrahedron0 c)
 427       vol = 0
 428
 429
 430 tetrahedron19 :: Cube -> Tetrahedron
 431 tetrahedron19 c =
 432     Tetrahedron fv' v0' v1' v2' v3' vol
 433     where
 434       v0' = center c
 435       v1' = center (right_face c)
 436       v2' = Face.v3 (right_face c)
 437       v3' = Face.v0 (right_face c)
 438       fv' = rotate ccwz $ rotate ccwy
 439                         $ Tetrahedron.fv (tetrahedron0 c)
 440       vol = 0
 441
 442
 443 tetrahedron20 :: Cube -> Tetrahedron
 444 tetrahedron20 c =
 445     Tetrahedron fv' v0' v1' v2' v3' vol
 446     where
 447       v0' = center c
 448       v1' = center (left_face c)
 449       v2' = Face.v0 (left_face c)
 450       v3' = Face.v1 (left_face c)
 451       fv' = rotate cwz (Tetrahedron.fv (tetrahedron0 c))
 452       vol = 0
 453
 454
 455 tetrahedron21 :: Cube -> Tetrahedron
 456 tetrahedron21 c =
 457     Tetrahedron fv' v0' v1' v2' v3' vol
 458     where
 459       v0' = center c
 460       v1' = center (left_face c)
 461       v2' = Face.v1 (left_face c)
 462       v3' = Face.v2 (left_face c)
 463       fv' = rotate cwz $ rotate ccwy $ Tetrahedron.fv (tetrahedron0 c)
 464       vol = 0
 465
 466
 467 tetrahedron22 :: Cube -> Tetrahedron
 468 tetrahedron22 c =
 469     Tetrahedron fv' v0' v1' v2' v3' vol
 470     where
 471       v0' = center c
 472       v1' = center (left_face c)
 473       v2' = Face.v2 (left_face c)
 474       v3' = Face.v3 (left_face c)
 475       fv' = rotate cwz $ rotate ccwy
 476                        $ rotate ccwy
 477                        $ Tetrahedron.fv (tetrahedron0 c)
 478       vol = 0
 479
 480
 481 tetrahedron23 :: Cube -> Tetrahedron
 482 tetrahedron23 c =
 483     Tetrahedron fv' v0' v1' v2' v3' vol
 484     where
 485       v0' = center c
 486       v1' = center (left_face c)
 487       v2' = Face.v3 (left_face c)
 488       v3' = Face.v0 (left_face c)
 489       fv' = rotate cwz $ rotate cwy
 490                        $ Tetrahedron.fv (tetrahedron0 c)
 491       vol = 0
 492
 493
 494 tetrahedra :: Cube -> [Tetrahedron]
 495 tetrahedra c =
 496     [tetrahedron0 c,
 497      tetrahedron1 c,
 498      tetrahedron2 c,
 499      tetrahedron3 c,
 500      tetrahedron4 c,
 501      tetrahedron5 c,
 502      tetrahedron6 c,
 503      tetrahedron7 c,
 504      tetrahedron8 c,
 505      tetrahedron9 c,
 506      tetrahedron10 c,
 507      tetrahedron11 c,
 508      tetrahedron12 c,
 509      tetrahedron13 c,
 510      tetrahedron14 c,
 511      tetrahedron15 c,
 512      tetrahedron16 c,
 513      tetrahedron17 c,
 514      tetrahedron18 c,
 515      tetrahedron19 c,
 516      tetrahedron20 c,
 517      tetrahedron21 c,
 518      tetrahedron22 c,
 519      tetrahedron23 c]
 520
 521 -- | All completely contained in the front half of the cube.
 522 front_half_tetrahedra :: Cube -> [Tetrahedron]
 523 front_half_tetrahedra c =
 524   [tetrahedron0 c,
 525    tetrahedron1 c,
 526    tetrahedron2 c,
 527    tetrahedron3 c,
 528    tetrahedron6 c,
 529    tetrahedron12 c,
 530    tetrahedron19 c,
 531    tetrahedron21 c]
 532
 533
 534 -- | All tetrahedra completely contained in the top half of the cube.
 535 top_half_tetrahedra :: Cube -> [Tetrahedron]
 536 top_half_tetrahedra c =
 537   [tetrahedron4 c,
 538    tetrahedron5 c,
 539    tetrahedron6 c,
 540    tetrahedron7 c,
 541    tetrahedron0 c,
 542    tetrahedron10 c,
 543    tetrahedron16 c,
 544    tetrahedron20 c]
 545
 546
 547 -- | All tetrahedra completely contained in the back half of the cube.
 548 back_half_tetrahedra :: Cube -> [Tetrahedron]
 549 back_half_tetrahedra c =
 550   [tetrahedron8 c,
 551    tetrahedron9 c,
 552    tetrahedron10 c,
 553    tetrahedron11 c,
 554    tetrahedron4 c,
 555    tetrahedron14 c,
 556    tetrahedron17 c,
 557    tetrahedron23 c]
 558
 559
 560 -- | All tetrahedra completely contained in the down half of the cube.
 561 down_half_tetrahedra :: Cube -> [Tetrahedron]
 562 down_half_tetrahedra c =
 563   [tetrahedron12 c,
 564    tetrahedron13 c,
 565    tetrahedron14 c,
 566    tetrahedron15 c,
 567    tetrahedron2 c,
 568    tetrahedron8 c,
 569    tetrahedron18 c,
 570    tetrahedron22 c]
 571
 572
 573 -- | All tetrahedra completely contained in the right half of the cube.
 574 right_half_tetrahedra :: Cube -> [Tetrahedron]
 575 right_half_tetrahedra c =
 576   [tetrahedron16 c,
 577    tetrahedron17 c,
 578    tetrahedron18 c,
 579    tetrahedron19 c,
 580    tetrahedron1 c,
 581    tetrahedron5 c,
 582    tetrahedron9 c,
 583    tetrahedron13 c]
 584
 585
 586 -- | All tetrahedra completely contained in the left half of the cube.
 587 left_half_tetrahedra :: Cube -> [Tetrahedron]
 588 left_half_tetrahedra c =
 589   [tetrahedron20 c,
 590    tetrahedron21 c,
 591    tetrahedron22 c,
 592    tetrahedron23 c,
 593    tetrahedron3 c,
 594    tetrahedron7 c,
 595    tetrahedron11 c,
 596    tetrahedron15 c]
 597
 598
 599 in_top_half :: Cube -> Point -> Bool
 600 in_top_half c (_,_,z) =
 601   distance_from_top <= distance_from_bottom
 602   where
 603     distance_from_top = abs $ (zmax c) - z
 604     distance_from_bottom = abs $ (zmin c) - z
 605
 606 in_front_half :: Cube -> Point -> Bool
 607 in_front_half c (x,_,_) =
 608     distance_from_front <= distance_from_back
 609   where
 610     distance_from_front = abs $ (xmin c) - x
 611     distance_from_back = abs $ (xmax c) - x
 612
 613
 614 in_left_half :: Cube -> Point -> Bool
 615 in_left_half c (_,y,_) =
 616     distance_from_left <= distance_from_right
 617   where
 618     distance_from_left = abs $ (ymin c) - y
 619     distance_from_right = abs $ (ymax c) - y
 620
 621
 622 -- | Takes a 'Cube', and returns the Tetrahedra belonging to it that
 623 --   contain the given 'Point'. This should be faster than checking
 624 --   every tetrahedron individually, since we determine which half
 625 --   (hemisphere?) of the cube the point lies in three times: once in
 626 --   each dimension. This allows us to eliminate non-candidates
 627 --   quickly.
 628 --
 629 --   This can throw an exception, but the use of 'head' might
 630 --   save us some unnecessary computations.
 631 --
 632 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
 633 find_containing_tetrahedron c p =
 634   head containing_tetrahedra
 635   where
 636     candidates = tetrahedra c
 637     non_candidates_x =
 638         if (in_front_half c p) then
 639           back_half_tetrahedra c
 640         else
 641           front_half_tetrahedra c
 642
 643     candidates_x = candidates \\ non_candidates_x
 644
 645     non_candidates_y =
 646       if (in_left_half c p) then
 647         right_half_tetrahedra c
 648       else
 649         left_half_tetrahedra c
 650
 651     candidates_xy = candidates_x \\ non_candidates_y
 652
 653     non_candidates_z =
 654       if (in_top_half c p) then
 655         down_half_tetrahedra c
 656       else
 657         top_half_tetrahedra c
 658
 659     candidates_xyz = candidates_xy \\ non_candidates_z
 660
 661     contains_our_point = flip contains_point p
 662     containing_tetrahedra = filter contains_our_point candidates_xyz