src/Cube.hs

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