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 tetrahedron :: Cube -> Int -> Tetrahedron
 203
 204 tetrahedron c 0 =
 205     Tetrahedron (Cube.fv c) v0' v1' v2' v3' vol 0
 206     where
 207       v0' = center c
 208       v1' = center (front_face c)
 209       v2' = Face.v0 (front_face c)
 210       v3' = Face.v1 (front_face c)
 211       vol = tetrahedra_volume c
 212
 213 tetrahedron c 1 =
 214     Tetrahedron fv' v0' v1' v2' v3' vol 1
 215     where
 216       v0' = center c
 217       v1' = center (front_face c)
 218       v2' = Face.v1 (front_face c)
 219       v3' = Face.v2 (front_face c)
 220       fv' = rotate ccwx (Cube.fv c)
 221       vol = tetrahedra_volume c
 222
 223 tetrahedron c 2 =
 224     Tetrahedron fv' v0' v1' v2' v3' vol 2
 225     where
 226       v0' = center c
 227       v1' = center (front_face c)
 228       v2' = Face.v2 (front_face c)
 229       v3' = Face.v3 (front_face c)
 230       fv' = rotate ccwx $ rotate ccwx $ Cube.fv c
 231       vol = tetrahedra_volume c
 232
 233 tetrahedron c 3 =
 234     Tetrahedron fv' v0' v1' v2' v3' vol 3
 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 = tetrahedra_volume c
 242
 243 tetrahedron c 4 =
 244     Tetrahedron fv' v0' v1' v2' v3' vol 4
 245     where
 246       v0' = center c
 247       v1' = center (top_face c)
 248       v2' = Face.v0 (top_face c)
 249       v3' = Face.v1 (top_face c)
 250       fv' = rotate cwy (Cube.fv c)
 251       vol = tetrahedra_volume c
 252
 253 tetrahedron c 5 =
 254     Tetrahedron fv' v0' v1' v2' v3' vol 5
 255     where
 256       v0' = center c
 257       v1' = center (top_face c)
 258       v2' = Face.v1 (top_face c)
 259       v3' = Face.v2 (top_face c)
 260       fv' = rotate cwy $ rotate cwz $ Tetrahedron.fv (tetrahedron c 0)
 261       vol = tetrahedra_volume c
 262
 263 tetrahedron c 6 =
 264     Tetrahedron fv' v0' v1' v2' v3' vol 6
 265     where
 266       v0' = center c
 267       v1' = center (top_face c)
 268       v2' = Face.v2 (top_face c)
 269       v3' = Face.v3 (top_face c)
 270       fv' = rotate cwy $ rotate cwz
 271                        $ rotate cwz
 272                        $ Tetrahedron.fv (tetrahedron c 0)
 273       vol = tetrahedra_volume c
 274
 275 tetrahedron c 7 =
 276     Tetrahedron fv' v0' v1' v2' v3' vol 7
 277     where
 278       v0' = center c
 279       v1' = center (top_face c)
 280       v2' = Face.v3 (top_face c)
 281       v3' = Face.v0 (top_face c)
 282       fv' = rotate cwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron c 0)
 283       vol = tetrahedra_volume c
 284
 285 tetrahedron c 8 =
 286     Tetrahedron fv' v0' v1' v2' v3' vol 8
 287     where
 288       v0' = center c
 289       v1' = center (back_face c)
 290       v2' = Face.v0 (back_face c)
 291       v3' = Face.v1 (back_face c)
 292       fv' = rotate cwy $ rotate cwy $ Tetrahedron.fv (tetrahedron c 0)
 293       vol = tetrahedra_volume c
 294
 295 tetrahedron c 9 =
 296     Tetrahedron fv' v0' v1' v2' v3' vol 9
 297     where
 298       v0' = center c
 299       v1' = center (back_face c)
 300       v2' = Face.v1 (back_face c)
 301       v3' = Face.v2 (back_face c)
 302       fv' = rotate cwy $ rotate cwy
 303                        $ rotate cwx
 304                        $ Tetrahedron.fv (tetrahedron c 0)
 305       vol = tetrahedra_volume c
 306
 307 tetrahedron c 10 =
 308     Tetrahedron fv' v0' v1' v2' v3' vol 10
 309     where
 310       v0' = center c
 311       v1' = center (back_face c)
 312       v2' = Face.v2 (back_face c)
 313       v3' = Face.v3 (back_face c)
 314       fv' = rotate cwy $ rotate cwy
 315                        $ rotate cwx
 316                        $ rotate cwx
 317                        $ Tetrahedron.fv (tetrahedron c 0)
 318
 319       vol = tetrahedra_volume c
 320
 321 tetrahedron c 11 =
 322     Tetrahedron fv' v0' v1' v2' v3' vol 11
 323     where
 324       v0' = center c
 325       v1' = center (back_face c)
 326       v2' = Face.v3 (back_face c)
 327       v3' = Face.v0 (back_face c)
 328       fv' = rotate cwy $ rotate cwy
 329                        $ rotate ccwx
 330                        $ Tetrahedron.fv (tetrahedron c 0)
 331       vol = tetrahedra_volume c
 332
 333 tetrahedron c 12 =
 334     Tetrahedron fv' v0' v1' v2' v3' vol 12
 335     where
 336       v0' = center c
 337       v1' = center (down_face c)
 338       v2' = Face.v0 (down_face c)
 339       v3' = Face.v1 (down_face c)
 340       fv' = rotate ccwy (Tetrahedron.fv (tetrahedron c 0))
 341       vol = tetrahedra_volume c
 342
 343 tetrahedron c 13 =
 344     Tetrahedron fv' v0' v1' v2' v3' vol 13
 345     where
 346       v0' = center c
 347       v1' = center (down_face c)
 348       v2' = Face.v1 (down_face c)
 349       v3' = Face.v2 (down_face c)
 350       fv' = rotate ccwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron c 0)
 351       vol = tetrahedra_volume c
 352
 353 tetrahedron c 14 =
 354     Tetrahedron fv' v0' v1' v2' v3' vol 14
 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 (tetrahedron c 0)
 363       vol = tetrahedra_volume c
 364
 365 tetrahedron c 15 =
 366     Tetrahedron fv' v0' v1' v2' v3' vol 15
 367     where
 368       v0' = center c
 369       v1' = center (down_face c)
 370       v2' = Face.v3 (down_face c)
 371       v3' = Face.v0 (down_face c)
 372       fv' = rotate ccwy $ rotate cwz $ Tetrahedron.fv (tetrahedron c 0)
 373       vol = tetrahedra_volume c
 374
 375 tetrahedron c 16 =
 376     Tetrahedron fv' v0' v1' v2' v3' vol 16
 377     where
 378       v0' = center c
 379       v1' = center (right_face c)
 380       v2' = Face.v0 (right_face c)
 381       v3' = Face.v1 (right_face c)
 382       fv' = rotate ccwz (Tetrahedron.fv (tetrahedron c 0))
 383       vol = tetrahedra_volume c
 384
 385 tetrahedron c 17 =
 386     Tetrahedron fv' v0' v1' v2' v3' vol 17
 387     where
 388       v0' = center c
 389       v1' = center (right_face c)
 390       v2' = Face.v1 (right_face c)
 391       v3' = Face.v2 (right_face c)
 392       fv' = rotate ccwz $ rotate cwy $ Tetrahedron.fv (tetrahedron c 0)
 393       vol = tetrahedra_volume c
 394
 395 tetrahedron c 18 =
 396     Tetrahedron fv' v0' v1' v2' v3' vol 18
 397     where
 398       v0' = center c
 399       v1' = center (right_face c)
 400       v2' = Face.v2 (right_face c)
 401       v3' = Face.v3 (right_face c)
 402       fv' = rotate ccwz $ rotate cwy
 403                         $ rotate cwy
 404                         $ Tetrahedron.fv (tetrahedron c 0)
 405       vol = tetrahedra_volume c
 406
 407 tetrahedron c 19 =
 408     Tetrahedron fv' v0' v1' v2' v3' vol 19
 409     where
 410       v0' = center c
 411       v1' = center (right_face c)
 412       v2' = Face.v3 (right_face c)
 413       v3' = Face.v0 (right_face c)
 414       fv' = rotate ccwz $ rotate ccwy
 415                         $ Tetrahedron.fv (tetrahedron c 0)
 416       vol = tetrahedra_volume c
 417
 418 tetrahedron c 20 =
 419     Tetrahedron fv' v0' v1' v2' v3' vol 20
 420     where
 421       v0' = center c
 422       v1' = center (left_face c)
 423       v2' = Face.v0 (left_face c)
 424       v3' = Face.v1 (left_face c)
 425       fv' = rotate cwz (Tetrahedron.fv (tetrahedron c 0))
 426       vol = tetrahedra_volume c
 427
 428 tetrahedron c 21 =
 429     Tetrahedron fv' v0' v1' v2' v3' vol 21
 430     where
 431       v0' = center c
 432       v1' = center (left_face c)
 433       v2' = Face.v1 (left_face c)
 434       v3' = Face.v2 (left_face c)
 435       fv' = rotate cwz $ rotate ccwy $ Tetrahedron.fv (tetrahedron c 0)
 436       vol = tetrahedra_volume c
 437
 438 tetrahedron c 22 =
 439     Tetrahedron fv' v0' v1' v2' v3' vol 22
 440     where
 441       v0' = center c
 442       v1' = center (left_face c)
 443       v2' = Face.v2 (left_face c)
 444       v3' = Face.v3 (left_face c)
 445       fv' = rotate cwz $ rotate ccwy
 446                        $ rotate ccwy
 447                        $ Tetrahedron.fv (tetrahedron c 0)
 448       vol = tetrahedra_volume c
 449
 450 tetrahedron c 23 =
 451     Tetrahedron fv' v0' v1' v2' v3' vol 23
 452     where
 453       v0' = center c
 454       v1' = center (left_face c)
 455       v2' = Face.v3 (left_face c)
 456       v3' = Face.v0 (left_face c)
 457       fv' = rotate cwz $ rotate cwy
 458                        $ Tetrahedron.fv (tetrahedron c 0)
 459       vol = tetrahedra_volume c
 460
 461 -- Feels dirty, but whatever.
 462 tetrahedron _ _ = error "asked for a nonexistent tetrahedron"
 463
 464
 465 tetrahedra :: Cube -> [Tetrahedron]
 466 tetrahedra c =
 467     [ tetrahedron c n | n <- [0..23] ]
 468
 469 -- | All completely contained in the front half of the cube.
 470 front_half_tetrahedra :: Cube -> [Tetrahedron]
 471 front_half_tetrahedra c =
 472   [ tetrahedron c n | n <- [0,1,2,3,6,12,19,21] ]
 473
 474 -- | All tetrahedra completely contained in the top half of the cube.
 475 top_half_tetrahedra :: Cube -> [Tetrahedron]
 476 top_half_tetrahedra c =
 477   [ tetrahedron c n | n <- [4,5,6,7,0,10,16,20] ]
 478
 479 -- | All tetrahedra completely contained in the back half of the cube.
 480 back_half_tetrahedra :: Cube -> [Tetrahedron]
 481 back_half_tetrahedra c =
 482   [ tetrahedron c n | n <- [8,9,10,11,4,14,17,23] ]
 483
 484 -- | All tetrahedra completely contained in the down half of the cube.
 485 down_half_tetrahedra :: Cube -> [Tetrahedron]
 486 down_half_tetrahedra c =
 487   [ tetrahedron c n | n <- [12,13,14,15,2,8,18,22] ]
 488
 489 -- | All tetrahedra completely contained in the right half of the cube.
 490 right_half_tetrahedra :: Cube -> [Tetrahedron]
 491 right_half_tetrahedra c =
 492   [ tetrahedron c n | n <- [16,17,18,19,1,5,9,13] ]
 493
 494 -- | All tetrahedra completely contained in the left half of the cube.
 495 left_half_tetrahedra :: Cube -> [Tetrahedron]
 496 left_half_tetrahedra c =
 497   [ tetrahedron c n | n <- [20,21,22,23,3,7,11,15] ]
 498
 499 in_top_half :: Cube -> Point -> Bool
 500 in_top_half c (_,_,z) =
 501   distance_from_top <= distance_from_bottom
 502   where
 503     distance_from_top = abs $ (zmax c) - z
 504     distance_from_bottom = abs $ (zmin c) - z
 505
 506 in_front_half :: Cube -> Point -> Bool
 507 in_front_half c (x,_,_) =
 508     distance_from_front <= distance_from_back
 509   where
 510     distance_from_front = abs $ (xmin c) - x
 511     distance_from_back = abs $ (xmax c) - x
 512
 513
 514 in_left_half :: Cube -> Point -> Bool
 515 in_left_half c (_,y,_) =
 516     distance_from_left <= distance_from_right
 517   where
 518     distance_from_left = abs $ (ymin c) - y
 519     distance_from_right = abs $ (ymax c) - y
 520
 521
 522 -- | Takes a 'Cube', and returns the Tetrahedra belonging to it that
 523 --   contain the given 'Point'. This should be faster than checking
 524 --   every tetrahedron individually, since we determine which half
 525 --   (hemisphere?) of the cube the point lies in three times: once in
 526 --   each dimension. This allows us to eliminate non-candidates
 527 --   quickly.
 528 --
 529 --   This can throw an exception, but the use of 'head' might
 530 --   save us some unnecessary computations.
 531 --
 532 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
 533 find_containing_tetrahedron c p =
 534   head containing_tetrahedra
 535   where
 536     candidates = tetrahedra c
 537     non_candidates_x =
 538         if (in_front_half c p) then
 539           back_half_tetrahedra c
 540         else
 541           front_half_tetrahedra c
 542
 543     candidates_x = candidates \\ non_candidates_x
 544
 545     non_candidates_y =
 546       if (in_left_half c p) then
 547         right_half_tetrahedra c
 548       else
 549         left_half_tetrahedra c
 550
 551     candidates_xy = candidates_x \\ non_candidates_y
 552
 553     non_candidates_z =
 554       if (in_top_half c p) then
 555         down_half_tetrahedra c
 556       else
 557         top_half_tetrahedra c
 558
 559     candidates_xyz = candidates_xy \\ non_candidates_z
 560
 561     contains_our_point = flip contains_point p
 562     containing_tetrahedra = filter contains_our_point candidates_xyz
 563