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'
 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
 209 tetrahedron1 :: Cube -> Tetrahedron
 210 tetrahedron1 c =
 211     Tetrahedron fv' v0' v1' v2' v3'
 212     where
 213       v0' = center c
 214       v1' = center (front_face c)
 215       v2' = Face.v1 (front_face c)
 216       v3' = Face.v2 (front_face c)
 217       fv' = rotate ccwx (Cube.fv c)
 218
 219 tetrahedron2 :: Cube -> Tetrahedron
 220 tetrahedron2 c =
 221     Tetrahedron fv' v0' v1' v2' v3'
 222     where
 223       v0' = center c
 224       v1' = center (front_face c)
 225       v2' = Face.v2 (front_face c)
 226       v3' = Face.v3 (front_face c)
 227       fv' = rotate ccwx $ rotate ccwx $ Cube.fv c
 228
 229 tetrahedron3 :: Cube -> Tetrahedron
 230 tetrahedron3 c =
 231     Tetrahedron fv' v0' v1' v2' v3'
 232     where
 233       v0' = center c
 234       v1' = center (front_face c)
 235       v2' = Face.v3 (front_face c)
 236       v3' = Face.v0 (front_face c)
 237       fv' = rotate cwx (Cube.fv c)
 238
 239 tetrahedron4 :: Cube -> Tetrahedron
 240 tetrahedron4 c =
 241     Tetrahedron fv' v0' v1' v2' v3'
 242     where
 243       v0' = center c
 244       v1' = center (top_face c)
 245       v2' = Face.v0 (top_face c)
 246       v3' = Face.v1 (top_face c)
 247       fv' = rotate cwy (Cube.fv c)
 248
 249 tetrahedron5 :: Cube -> Tetrahedron
 250 tetrahedron5 c =
 251     Tetrahedron fv' v0' v1' v2' v3'
 252     where
 253       v0' = center c
 254       v1' = center (top_face c)
 255       v2' = Face.v1 (top_face c)
 256       v3' = Face.v2 (top_face c)
 257       fv' = rotate cwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
 258
 259 tetrahedron6 :: Cube -> Tetrahedron
 260 tetrahedron6 c =
 261     Tetrahedron fv' v0' v1' v2' v3'
 262     where
 263       v0' = center c
 264       v1' = center (top_face c)
 265       v2' = Face.v2 (top_face c)
 266       v3' = Face.v3 (top_face c)
 267       fv' = rotate cwy $ rotate cwz
 268                        $ rotate cwz
 269                        $ Tetrahedron.fv (tetrahedron0 c)
 270
 271 tetrahedron7 :: Cube -> Tetrahedron
 272 tetrahedron7 c =
 273     Tetrahedron fv' v0' v1' v2' v3'
 274     where
 275       v0' = center c
 276       v1' = center (top_face c)
 277       v2' = Face.v3 (top_face c)
 278       v3' = Face.v0 (top_face c)
 279       fv' = rotate cwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
 280
 281 tetrahedron8 :: Cube -> Tetrahedron
 282 tetrahedron8 c =
 283     Tetrahedron fv' v0' v1' v2' v3'
 284     where
 285       v0' = center c
 286       v1' = center (back_face c)
 287       v2' = Face.v0 (back_face c)
 288       v3' = Face.v1 (back_face c)
 289       fv' = rotate cwy $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)
 290
 291 tetrahedron9 :: Cube -> Tetrahedron
 292 tetrahedron9 c =
 293     Tetrahedron fv' v0' v1' v2' v3'
 294     where
 295       v0' = center c
 296       v1' = center (back_face c)
 297       v2' = Face.v1 (back_face c)
 298       v3' = Face.v2 (back_face c)
 299       fv' = rotate cwy $ rotate cwy
 300                        $ rotate cwx
 301                        $ Tetrahedron.fv (tetrahedron0 c)
 302
 303 tetrahedron10 :: Cube -> Tetrahedron
 304 tetrahedron10 c =
 305     Tetrahedron fv' v0' v1' v2' v3'
 306     where
 307       v0' = center c
 308       v1' = center (back_face c)
 309       v2' = Face.v2 (back_face c)
 310       v3' = Face.v3 (back_face c)
 311       fv' = rotate cwy $ rotate cwy
 312                        $ rotate cwx
 313                        $ rotate cwx
 314                        $ Tetrahedron.fv (tetrahedron0 c)
 315
 316
 317 tetrahedron11 :: Cube -> Tetrahedron
 318 tetrahedron11 c =
 319     Tetrahedron fv' v0' v1' v2' v3'
 320     where
 321       v0' = center c
 322       v1' = center (back_face c)
 323       v2' = Face.v3 (back_face c)
 324       v3' = Face.v0 (back_face c)
 325       fv' = rotate cwy $ rotate cwy
 326                        $ rotate ccwx
 327                        $ Tetrahedron.fv (tetrahedron0 c)
 328
 329
 330 tetrahedron12 :: Cube -> Tetrahedron
 331 tetrahedron12 c =
 332     Tetrahedron fv' v0' v1' v2' v3'
 333     where
 334       v0' = center c
 335       v1' = center (down_face c)
 336       v2' = Face.v0 (down_face c)
 337       v3' = Face.v1 (down_face c)
 338       fv' = rotate ccwy (Tetrahedron.fv (tetrahedron0 c))
 339
 340
 341 tetrahedron13 :: Cube -> Tetrahedron
 342 tetrahedron13 c =
 343     Tetrahedron fv' v0' v1' v2' v3'
 344     where
 345       v0' = center c
 346       v1' = center (down_face c)
 347       v2' = Face.v1 (down_face c)
 348       v3' = Face.v2 (down_face c)
 349       fv' = rotate ccwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
 350
 351
 352 tetrahedron14 :: Cube -> Tetrahedron
 353 tetrahedron14 c =
 354     Tetrahedron fv' v0' v1' v2' v3'
 355     where
 356       v0' = center c
 357       v1' = center (down_face c)
 358       v2' = Face.v2 (down_face c)
 359       v3' = Face.v3 (down_face c)
 360       fv' = rotate ccwy $ rotate ccwz
 361                         $ rotate ccwz
 362                         $ Tetrahedron.fv (tetrahedron0 c)
 363
 364
 365 tetrahedron15 :: Cube -> Tetrahedron
 366 tetrahedron15 c =
 367     Tetrahedron fv' v0' v1' v2' v3'
 368     where
 369       v0' = center c
 370       v1' = center (down_face c)
 371       v2' = Face.v3 (down_face c)
 372       v3' = Face.v0 (down_face c)
 373       fv' = rotate ccwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
 374
 375
 376 tetrahedron16 :: Cube -> Tetrahedron
 377 tetrahedron16 c =
 378     Tetrahedron fv' v0' v1' v2' v3'
 379     where
 380       v0' = center c
 381       v1' = center (right_face c)
 382       v2' = Face.v0 (right_face c)
 383       v3' = Face.v1 (right_face c)
 384       fv' = rotate ccwz (Tetrahedron.fv (tetrahedron0 c))
 385
 386
 387 tetrahedron17 :: Cube -> Tetrahedron
 388 tetrahedron17 c =
 389     Tetrahedron fv' v0' v1' v2' v3'
 390     where
 391       v0' = center c
 392       v1' = center (right_face c)
 393       v2' = Face.v1 (right_face c)
 394       v3' = Face.v2 (right_face c)
 395       fv' = rotate ccwz $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)
 396
 397
 398 tetrahedron18 :: Cube -> Tetrahedron
 399 tetrahedron18 c =
 400     Tetrahedron fv' v0' v1' v2' v3'
 401     where
 402       v0' = center c
 403       v1' = center (right_face c)
 404       v2' = Face.v2 (right_face c)
 405       v3' = Face.v3 (right_face c)
 406       fv' = rotate ccwz $ rotate cwy
 407                         $ rotate cwy
 408                         $ Tetrahedron.fv (tetrahedron0 c)
 409
 410
 411 tetrahedron19 :: Cube -> Tetrahedron
 412 tetrahedron19 c =
 413     Tetrahedron fv' v0' v1' v2' v3'
 414     where
 415       v0' = center c
 416       v1' = center (right_face c)
 417       v2' = Face.v3 (right_face c)
 418       v3' = Face.v0 (right_face c)
 419       fv' = rotate ccwz $ rotate ccwy
 420                         $ Tetrahedron.fv (tetrahedron0 c)
 421
 422
 423 tetrahedron20 :: Cube -> Tetrahedron
 424 tetrahedron20 c =
 425     Tetrahedron fv' v0' v1' v2' v3'
 426     where
 427       v0' = center c
 428       v1' = center (left_face c)
 429       v2' = Face.v0 (left_face c)
 430       v3' = Face.v1 (left_face c)
 431       fv' = rotate cwz (Tetrahedron.fv (tetrahedron0 c))
 432
 433
 434 tetrahedron21 :: Cube -> Tetrahedron
 435 tetrahedron21 c =
 436     Tetrahedron fv' v0' v1' v2' v3'
 437     where
 438       v0' = center c
 439       v1' = center (left_face c)
 440       v2' = Face.v1 (left_face c)
 441       v3' = Face.v2 (left_face c)
 442       fv' = rotate cwz $ rotate ccwy $ Tetrahedron.fv (tetrahedron0 c)
 443
 444
 445 tetrahedron22 :: Cube -> Tetrahedron
 446 tetrahedron22 c =
 447     Tetrahedron fv' v0' v1' v2' v3'
 448     where
 449       v0' = center c
 450       v1' = center (left_face c)
 451       v2' = Face.v2 (left_face c)
 452       v3' = Face.v3 (left_face c)
 453       fv' = rotate cwz $ rotate ccwy
 454                        $ rotate ccwy
 455                        $ Tetrahedron.fv (tetrahedron0 c)
 456
 457
 458 tetrahedron23 :: Cube -> Tetrahedron
 459 tetrahedron23 c =
 460     Tetrahedron fv' v0' v1' v2' v3'
 461     where
 462       v0' = center c
 463       v1' = center (left_face c)
 464       v2' = Face.v3 (left_face c)
 465       v3' = Face.v0 (left_face c)
 466       fv' = rotate cwz $ rotate cwy
 467                        $ Tetrahedron.fv (tetrahedron0 c)
 468
 469
 470 tetrahedra :: Cube -> [Tetrahedron]
 471 tetrahedra c =
 472     [tetrahedron0 c,
 473      tetrahedron1 c,
 474      tetrahedron2 c,
 475      tetrahedron3 c,
 476      tetrahedron4 c,
 477      tetrahedron5 c,
 478      tetrahedron6 c,
 479      tetrahedron7 c,
 480      tetrahedron8 c,
 481      tetrahedron9 c,
 482      tetrahedron10 c,
 483      tetrahedron11 c,
 484      tetrahedron12 c,
 485      tetrahedron13 c,
 486      tetrahedron14 c,
 487      tetrahedron15 c,
 488      tetrahedron16 c,
 489      tetrahedron17 c,
 490      tetrahedron18 c,
 491      tetrahedron19 c,
 492      tetrahedron20 c,
 493      tetrahedron21 c,
 494      tetrahedron22 c,
 495      tetrahedron23 c]
 496
 497 -- | All completely contained in the front half of the cube.
 498 front_half_tetrahedra :: Cube -> [Tetrahedron]
 499 front_half_tetrahedra c =
 500   [tetrahedron0 c,
 501    tetrahedron1 c,
 502    tetrahedron2 c,
 503    tetrahedron3 c,
 504    tetrahedron6 c,
 505    tetrahedron12 c,
 506    tetrahedron19 c,
 507    tetrahedron21 c]
 508
 509
 510 -- | All tetrahedra completely contained in the top half of the cube.
 511 top_half_tetrahedra :: Cube -> [Tetrahedron]
 512 top_half_tetrahedra c =
 513   [tetrahedron4 c,
 514    tetrahedron5 c,
 515    tetrahedron6 c,
 516    tetrahedron7 c,
 517    tetrahedron0 c,
 518    tetrahedron10 c,
 519    tetrahedron16 c,
 520    tetrahedron20 c]
 521
 522
 523 -- | All tetrahedra completely contained in the back half of the cube.
 524 back_half_tetrahedra :: Cube -> [Tetrahedron]
 525 back_half_tetrahedra c =
 526   [tetrahedron8 c,
 527    tetrahedron9 c,
 528    tetrahedron10 c,
 529    tetrahedron11 c,
 530    tetrahedron4 c,
 531    tetrahedron14 c,
 532    tetrahedron17 c,
 533    tetrahedron23 c]
 534
 535
 536 -- | All tetrahedra completely contained in the down half of the cube.
 537 down_half_tetrahedra :: Cube -> [Tetrahedron]
 538 down_half_tetrahedra c =
 539   [tetrahedron12 c,
 540    tetrahedron13 c,
 541    tetrahedron14 c,
 542    tetrahedron15 c,
 543    tetrahedron2 c,
 544    tetrahedron8 c,
 545    tetrahedron18 c,
 546    tetrahedron22 c]
 547
 548
 549 -- | All tetrahedra completely contained in the right half of the cube.
 550 right_half_tetrahedra :: Cube -> [Tetrahedron]
 551 right_half_tetrahedra c =
 552   [tetrahedron16 c,
 553    tetrahedron17 c,
 554    tetrahedron18 c,
 555    tetrahedron19 c,
 556    tetrahedron1 c,
 557    tetrahedron5 c,
 558    tetrahedron9 c,
 559    tetrahedron13 c]
 560
 561
 562 -- | All tetrahedra completely contained in the left half of the cube.
 563 left_half_tetrahedra :: Cube -> [Tetrahedron]
 564 left_half_tetrahedra c =
 565   [tetrahedron20 c,
 566    tetrahedron21 c,
 567    tetrahedron22 c,
 568    tetrahedron23 c,
 569    tetrahedron3 c,
 570    tetrahedron7 c,
 571    tetrahedron11 c,
 572    tetrahedron15 c]
 573
 574
 575 in_top_half :: Cube -> Point -> Bool
 576 in_top_half c (_,_,z) =
 577   distance_from_top <= distance_from_bottom
 578   where
 579     distance_from_top = abs $ (zmax c) - z
 580     distance_from_bottom = abs $ (zmin c) - z
 581
 582 in_front_half :: Cube -> Point -> Bool
 583 in_front_half c (x,_,_) =
 584     distance_from_front <= distance_from_back
 585   where
 586     distance_from_front = abs $ (xmin c) - x
 587     distance_from_back = abs $ (xmax c) - x
 588
 589
 590 in_left_half :: Cube -> Point -> Bool
 591 in_left_half c (_,y,_) =
 592     distance_from_left <= distance_from_right
 593   where
 594     distance_from_left = abs $ (ymin c) - y
 595     distance_from_right = abs $ (ymax c) - y
 596
 597
 598 -- | Takes a 'Cube', and returns the Tetrahedra belonging to it that
 599 --   contain the given 'Point'. This should be faster than checking
 600 --   every tetrahedron individually, since we determine which half
 601 --   (hemisphere?) of the cube the point lies in three times: once in
 602 --   each dimension. This allows us to eliminate non-candidates
 603 --   quickly.
 604 --
 605 --   This can throw an exception, but the use of 'head' might
 606 --   save us some unnecessary computations.
 607 --
 608 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
 609 find_containing_tetrahedron c p =
 610   head containing_tetrahedra
 611   where
 612     candidates = tetrahedra c
 613     non_candidates_x =
 614         if (in_front_half c p) then
 615           back_half_tetrahedra c
 616         else
 617           front_half_tetrahedra c
 618
 619     candidates_x = candidates \\ non_candidates_x
 620
 621     non_candidates_y =
 622       if (in_left_half c p) then
 623         right_half_tetrahedra c
 624       else
 625         left_half_tetrahedra c
 626
 627     candidates_xy = candidates_x \\ non_candidates_y
 628
 629     non_candidates_z =
 630       if (in_top_half c p) then
 631         down_half_tetrahedra c
 632       else
 633         top_half_tetrahedra c
 634
 635     candidates_xyz = candidates_xy \\ non_candidates_z
 636
 637     contains_our_point = flip contains_point p
 638     containing_tetrahedra = filter contains_our_point candidates_xyz