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