src/Cube.hs

   1 module Cube
   2 where
   3
   4 import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
   5
   6 import Cardinal
   7 import qualified Face (Face(Face, v0, v1, v2, v3))
   8 import FunctionValues
   9 import Point
  10 import Tetrahedron hiding (c)
  11 import ThreeDimensional
  12
  13 data Cube = Cube { h :: Double,
  14                    i :: Int,
  15                    j :: Int,
  16                    k :: Int,
  17                    fv :: FunctionValues }
  18             deriving (Eq)
  19
  20
  21 instance Arbitrary Cube where
  22     arbitrary = do
  23       (Positive h') <- arbitrary :: Gen (Positive Double)
  24       i' <- choose (coordmin, coordmax)
  25       j' <- choose (coordmin, coordmax)
  26       k' <- choose (coordmin, coordmax)
  27       fv' <- arbitrary :: Gen FunctionValues
  28       return (Cube h' i' j' k' fv')
  29         where
  30           coordmin = -268435456 -- -(2^29 / 2)
  31           coordmax = 268435456  -- +(2^29 / 2)
  32
  33
  34 instance Show Cube where
  35     show c =
  36         "Cube_" ++ subscript ++ "\n" ++
  37         " h: " ++ (show (h c)) ++ "\n" ++
  38         " Center: " ++ (show (center c)) ++ "\n" ++
  39         " xmin: " ++ (show (xmin c)) ++ "\n" ++
  40         " xmax: " ++ (show (xmax c)) ++ "\n" ++
  41         " ymin: " ++ (show (ymin c)) ++ "\n" ++
  42         " ymax: " ++ (show (ymax c)) ++ "\n" ++
  43         " zmin: " ++ (show (zmin c)) ++ "\n" ++
  44         " zmax: " ++ (show (zmax c)) ++ "\n" ++
  45         " fv: " ++ (show (Cube.fv c)) ++ "\n"
  46         where
  47           subscript =
  48               (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k c))
  49
  50
  51 -- | Returns an empty 'Cube'.
  52 empty_cube :: Cube
  53 empty_cube = Cube 0 0 0 0 empty_values
  54
  55
  56 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
  57 --   p. 76.
  58 xmin :: Cube -> Double
  59 xmin c = (2*i' - 1)*delta / 2
  60     where
  61       i' = fromIntegral (i c) :: Double
  62       delta = h c
  63
  64 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
  65 --   p. 76.
  66 xmax :: Cube -> Double
  67 xmax c = (2*i' + 1)*delta / 2
  68     where
  69       i' = fromIntegral (i c) :: Double
  70       delta = h c
  71
  72 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
  73 --   p. 76.
  74 ymin :: Cube -> Double
  75 ymin c = (2*j' - 1)*delta / 2
  76     where
  77       j' = fromIntegral (j c) :: Double
  78       delta = h c
  79
  80 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
  81 --   p. 76.
  82 ymax :: Cube -> Double
  83 ymax c = (2*j' + 1)*delta / 2
  84     where
  85       j' = fromIntegral (j c) :: Double
  86       delta = h c
  87
  88 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
  89 --   p. 76.
  90 zmin :: Cube -> Double
  91 zmin c = (2*k' - 1)*delta / 2
  92     where
  93       k' = fromIntegral (k c) :: Double
  94       delta = h c
  95
  96 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
  97 --   p. 76.
  98 zmax :: Cube -> Double
  99 zmax c = (2*k' + 1)*delta / 2
 100     where
 101       k' = fromIntegral (k c) :: Double
 102       delta = h c
 103
 104 instance ThreeDimensional Cube where
 105     -- | The center of Cube_ijk coincides with v_ijk at
 106     --   (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
 107     center c = (x, y, z)
 108            where
 109              delta = h c
 110              i' = fromIntegral (i c) :: Double
 111              j' = fromIntegral (j c) :: Double
 112              k' = fromIntegral (k c) :: Double
 113              x = delta * i'
 114              y = delta * j'
 115              z = delta * k'
 116
 117     -- | It's easy to tell if a point is within a cube; just make sure
 118     --   that it falls on the proper side of each of the cube's faces.
 119     contains_point c p
 120         | (x_coord p) < (xmin c) = False
 121         | (x_coord p) > (xmax c) = False
 122         | (y_coord p) < (ymin c) = False
 123         | (y_coord p) > (ymax c) = False
 124         | (z_coord p) < (zmin c) = False
 125         | (z_coord p) > (zmax c) = False
 126         | otherwise = True
 127
 128
 129
 130 -- Face stuff.
 131
 132 -- | The top (in the direction of z) face of the cube.
 133 top_face :: Cube -> Face.Face
 134 top_face c = Face.Face v0' v1' v2' v3'
 135     where
 136       delta = (1/2)*(h c)
 137       v0' = (center c) + (delta, -delta, delta)
 138       v1' = (center c) + (delta, delta, delta)
 139       v2' = (center c) + (-delta, delta, delta)
 140       v3' = (center c) + (-delta, -delta, delta)
 141
 142
 143
 144 -- | The back (in the direction of x) face of the cube.
 145 back_face :: Cube -> Face.Face
 146 back_face c = Face.Face v0' v1' v2' v3'
 147     where
 148       delta = (1/2)*(h c)
 149       v0' = (center c) + (delta, -delta, -delta)
 150       v1' = (center c) + (delta, delta, -delta)
 151       v2' = (center c) + (delta, delta, delta)
 152       v3' = (center c) + (delta, -delta, delta)
 153
 154
 155 -- The bottom face (in the direction of -z) of the cube.
 156 down_face :: Cube -> Face.Face
 157 down_face c = Face.Face v0' v1' v2' v3'
 158     where
 159       delta = (1/2)*(h c)
 160       v0' = (center c) + (-delta, -delta, -delta)
 161       v1' = (center c) + (-delta, delta, -delta)
 162       v2' = (center c) + (delta, delta, -delta)
 163       v3' = (center c) + (delta, -delta, -delta)
 164
 165
 166
 167 -- | The front (in the direction of -x) face of the cube.
 168 front_face :: Cube -> Face.Face
 169 front_face c = Face.Face v0' v1' v2' v3'
 170     where
 171       delta = (1/2)*(h c)
 172       v0' = (center c) + (-delta, -delta, delta)
 173       v1' = (center c) + (-delta, delta, delta)
 174       v2' = (center c) + (-delta, delta, -delta)
 175       v3' = (center c) + (-delta, -delta, -delta)
 176
 177 -- | The left (in the direction of -y) face of the cube.
 178 left_face :: Cube -> Face.Face
 179 left_face c = Face.Face v0' v1' v2' v3'
 180     where
 181       delta = (1/2)*(h c)
 182       v0' = (center c) + (delta, -delta, delta)
 183       v1' = (center c) + (-delta, -delta, delta)
 184       v2' = (center c) + (-delta, -delta, -delta)
 185       v3' = (center c) + (delta, -delta, -delta)
 186
 187
 188 -- | The right (in the direction of y) face of the cube.
 189 right_face :: Cube -> Face.Face
 190 right_face c = Face.Face v0' v1' v2' v3'
 191     where
 192       delta = (1/2)*(h c)
 193       v0' = (center c) + (-delta, delta, delta)
 194       v1' = (center c) + (delta, delta, delta)
 195       v2' = (center c) + (delta, delta, -delta)
 196       v3' = (center c) + (-delta, delta, -delta)
 197
 198
 199 tetrahedron0 :: Cube -> Tetrahedron
 200 tetrahedron0 c =
 201     Tetrahedron (Cube.fv c) v0' v1' v2' v3'
 202     where
 203       v0' = center c
 204       v1' = center (front_face c)
 205       v2' = Face.v0 (front_face c)
 206       v3' = Face.v1 (front_face c)
 207
 208 tetrahedron1 :: Cube -> Tetrahedron
 209 tetrahedron1 c =
 210     Tetrahedron fv' v0' v1' v2' v3'
 211     where
 212       v0' = center c
 213       v1' = center (front_face c)
 214       v2' = Face.v1 (front_face c)
 215       v3' = Face.v2 (front_face c)
 216       fv' = rotate ccwx (Cube.fv c)
 217
 218 tetrahedron2 :: Cube -> Tetrahedron
 219 tetrahedron2 c =
 220     Tetrahedron fv' v0' v1' v2' v3'
 221     where
 222       v0' = center c
 223       v1' = center (front_face c)
 224       v2' = Face.v2 (front_face c)
 225       v3' = Face.v3 (front_face c)
 226       fv' = rotate ccwx $ rotate ccwx $ Cube.fv c
 227
 228 tetrahedron3 :: Cube -> Tetrahedron
 229 tetrahedron3 c =
 230     Tetrahedron fv' v0' v1' v2' v3'
 231     where
 232       v0' = center c
 233       v1' = center (front_face c)
 234       v2' = Face.v3 (front_face c)
 235       v3' = Face.v0 (front_face c)
 236       fv' = rotate cwx (Cube.fv c)
 237
 238 tetrahedron4 :: Cube -> Tetrahedron
 239 tetrahedron4 c =
 240     Tetrahedron fv' v0' v1' v2' v3'
 241     where
 242       v0' = center c
 243       v1' = center (top_face c)
 244       v2' = Face.v0 (top_face c)
 245       v3' = Face.v1 (top_face c)
 246       fv' = rotate cwy (Cube.fv c)
 247
 248 tetrahedron5 :: Cube -> Tetrahedron
 249 tetrahedron5 c =
 250     Tetrahedron fv' v0' v1' v2' v3'
 251     where
 252       v0' = center c
 253       v1' = center (top_face c)
 254       v2' = Face.v1 (top_face c)
 255       v3' = Face.v2 (top_face c)
 256       fv' = rotate cwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
 257
 258 tetrahedron6 :: Cube -> Tetrahedron
 259 tetrahedron6 c =
 260     Tetrahedron fv' v0' v1' v2' v3'
 261     where
 262       v0' = center c
 263       v1' = center (top_face c)
 264       v2' = Face.v2 (top_face c)
 265       v3' = Face.v3 (top_face c)
 266       fv' = rotate cwy $ rotate cwz $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
 267
 268 tetrahedron7 :: Cube -> Tetrahedron
 269 tetrahedron7 c =
 270     Tetrahedron fv' v0' v1' v2' v3'
 271     where
 272       v0' = center c
 273       v1' = center (top_face c)
 274       v2' = Face.v3 (top_face c)
 275       v3' = Face.v0 (top_face c)
 276       fv' = rotate cwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
 277
 278 tetrahedron8 :: Cube -> Tetrahedron
 279 tetrahedron8 c =
 280     Tetrahedron fv' v0' v1' v2' v3'
 281     where
 282       v0' = center c
 283       v1' = center (back_face c)
 284       v2' = Face.v0 (back_face c)
 285       v3' = Face.v1 (back_face c)
 286       fv' = rotate cwy $ rotate cwy $ (Tetrahedron.fv (tetrahedron0 c))
 287
 288 tetrahedron9 :: Cube -> Tetrahedron
 289 tetrahedron9 c =
 290     Tetrahedron fv' v0' v1' v2' v3'
 291     where
 292       v0' = center c
 293       v1' = center (back_face c)
 294       v2' = Face.v1 (back_face c)
 295       v3' = Face.v2 (back_face c)
 296       fv' = rotate cwy $ rotate cwy $ rotate cwx $ Tetrahedron.fv (tetrahedron0 c)
 297
 298 tetrahedron10 :: Cube -> Tetrahedron
 299 tetrahedron10 c =
 300     Tetrahedron fv' v0' v1' v2' v3'
 301     where
 302       v0' = center c
 303       v1' = center (back_face c)
 304       v2' = Face.v2 (back_face c)
 305       v3' = Face.v3 (back_face c)
 306       fv' = rotate cwy $ rotate cwy
 307                        $ rotate cwx
 308                        $ rotate cwx
 309                        $ Tetrahedron.fv (tetrahedron0 c)
 310
 311
 312 tetrahedron11 :: Cube -> Tetrahedron
 313 tetrahedron11 c =
 314     Tetrahedron fv' v0' v1' v2' v3'
 315     where
 316       v0' = center c
 317       v1' = center (back_face c)
 318       v2' = Face.v3 (back_face c)
 319       v3' = Face.v0 (back_face c)
 320       fv' = rotate cwy $ rotate cwy
 321                        $ rotate ccwx
 322                        $ Tetrahedron.fv (tetrahedron0 c)
 323
 324
 325 tetrahedron12 :: Cube -> Tetrahedron
 326 tetrahedron12 c =
 327     Tetrahedron fv' v0' v1' v2' v3'
 328     where
 329       v0' = center c
 330       v1' = center (down_face c)
 331       v2' = Face.v0 (down_face c)
 332       v3' = Face.v1 (down_face c)
 333       fv' = rotate ccwy (Tetrahedron.fv (tetrahedron0 c))
 334
 335
 336 tetrahedron13 :: Cube -> Tetrahedron
 337 tetrahedron13 c =
 338     Tetrahedron fv' v0' v1' v2' v3'
 339     where
 340       v0' = center c
 341       v1' = center (down_face c)
 342       v2' = Face.v1 (down_face c)
 343       v3' = Face.v2 (down_face c)
 344       fv' = rotate ccwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
 345
 346
 347 tetrahedron14 :: Cube -> Tetrahedron
 348 tetrahedron14 c =
 349     Tetrahedron fv' v0' v1' v2' v3'
 350     where
 351       v0' = center c
 352       v1' = center (down_face c)
 353       v2' = Face.v2 (down_face c)
 354       v3' = Face.v3 (down_face c)
 355       fv' = rotate ccwy $ rotate ccwz
 356                         $ rotate ccwz
 357                         $ Tetrahedron.fv (tetrahedron0 c)
 358
 359
 360 tetrahedron15 :: Cube -> Tetrahedron
 361 tetrahedron15 c =
 362     Tetrahedron fv' v0' v1' v2' v3'
 363     where
 364       v0' = center c
 365       v1' = center (down_face c)
 366       v2' = Face.v3 (down_face c)
 367       v3' = Face.v0 (down_face c)
 368       fv' = rotate ccwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
 369
 370
 371 tetrahedron16 :: Cube -> Tetrahedron
 372 tetrahedron16 c =
 373     Tetrahedron fv' v0' v1' v2' v3'
 374     where
 375       v0' = center c
 376       v1' = center (right_face c)
 377       v2' = Face.v0 (right_face c)
 378       v3' = Face.v1 (right_face c)
 379       fv' = rotate ccwz (Tetrahedron.fv (tetrahedron0 c))
 380
 381
 382 tetrahedron17 :: Cube -> Tetrahedron
 383 tetrahedron17 c =
 384     Tetrahedron fv' v0' v1' v2' v3'
 385     where
 386       v0' = center c
 387       v1' = center (right_face c)
 388       v2' = Face.v1 (right_face c)
 389       v3' = Face.v2 (right_face c)
 390       fv' = rotate ccwz $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)
 391
 392
 393 tetrahedron18 :: Cube -> Tetrahedron
 394 tetrahedron18 c =
 395     Tetrahedron fv' v0' v1' v2' v3'
 396     where
 397       v0' = center c
 398       v1' = center (right_face c)
 399       v2' = Face.v2 (right_face c)
 400       v3' = Face.v3 (right_face c)
 401       fv' = rotate ccwz $ rotate cwy
 402                         $ rotate cwy
 403                         $ Tetrahedron.fv (tetrahedron0 c)
 404
 405
 406 tetrahedron19 :: Cube -> Tetrahedron
 407 tetrahedron19 c =
 408     Tetrahedron fv' v0' v1' v2' v3'
 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 (tetrahedron0 c)
 416
 417
 418 tetrahedron20 :: Cube -> Tetrahedron
 419 tetrahedron20 c =
 420     Tetrahedron fv' v0' v1' v2' v3'
 421     where
 422       v0' = center c
 423       v1' = center (left_face c)
 424       v2' = Face.v0 (left_face c)
 425       v3' = Face.v1 (left_face c)
 426       fv' = rotate cwz (Tetrahedron.fv (tetrahedron0 c))
 427
 428
 429 tetrahedron21 :: Cube -> Tetrahedron
 430 tetrahedron21 c =
 431     Tetrahedron fv' v0' v1' v2' v3'
 432     where
 433       v0' = center c
 434       v1' = center (left_face c)
 435       v2' = Face.v1 (left_face c)
 436       v3' = Face.v2 (left_face c)
 437       fv' = rotate cwz $ rotate ccwy $ Tetrahedron.fv (tetrahedron0 c)
 438
 439
 440 tetrahedron22 :: Cube -> Tetrahedron
 441 tetrahedron22 c =
 442     Tetrahedron fv' v0' v1' v2' v3'
 443     where
 444       v0' = center c
 445       v1' = center (left_face c)
 446       v2' = Face.v2 (left_face c)
 447       v3' = Face.v3 (left_face c)
 448       fv' = rotate cwz $ rotate ccwy
 449                        $ rotate ccwy
 450                        $ Tetrahedron.fv (tetrahedron0 c)
 451
 452
 453 tetrahedron23 :: Cube -> Tetrahedron
 454 tetrahedron23 c =
 455     Tetrahedron fv' v0' v1' v2' v3'
 456     where
 457       v0' = center c
 458       v1' = center (left_face c)
 459       v2' = Face.v3 (left_face c)
 460       v3' = Face.v0 (left_face c)
 461       fv' = rotate cwz $ rotate cwy
 462                        $ Tetrahedron.fv (tetrahedron0 c)
 463
 464
 465 tetrahedrons :: Cube -> [Tetrahedron]
 466 tetrahedrons c =
 467     [tetrahedron0 c,
 468      tetrahedron1 c,
 469      tetrahedron2 c,
 470      tetrahedron3 c,
 471      tetrahedron4 c,
 472      tetrahedron5 c,
 473      tetrahedron6 c,
 474      tetrahedron7 c,
 475      tetrahedron8 c,
 476      tetrahedron9 c,
 477      tetrahedron10 c,
 478      tetrahedron11 c,
 479      tetrahedron12 c,
 480      tetrahedron13 c,
 481      tetrahedron14 c,
 482      tetrahedron15 c,
 483      tetrahedron16 c,
 484      tetrahedron17 c,
 485      tetrahedron18 c,
 486      tetrahedron19 c,
 487      tetrahedron20 c,
 488      tetrahedron21 c,
 489      tetrahedron22 c,
 490      tetrahedron23 c]
 491
 492
 493 -- | Takes a 'Cube', and returns all Tetrahedra belonging to it that
 494 --   contain the given 'Point'.
 495 find_containing_tetrahedra :: Cube -> Point -> [Tetrahedron]
 496 find_containing_tetrahedra c p =
 497     filter contains_our_point all_tetrahedra
 498     where
 499       contains_our_point = flip contains_point p
 500       all_tetrahedra = tetrahedrons c