src/Cube.hs

   1 module Cube
   2 where
   3
   4 import Cardinal
   5 import Face (Face(Face, v0, v1, v2, v3))
   6 import FunctionValues
   7 import Point
   8 import Tetrahedron (Tetrahedron(Tetrahedron), fv)
   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_" ++ (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k c)) ++
  22         " (Center: " ++ (show (center c)) ++ ")" ++
  23         " (xmin: " ++ (show (xmin c)) ++ ")" ++
  24         " (xmax: " ++ (show (xmax c)) ++ ")" ++
  25         " (ymin: " ++ (show (ymin c)) ++ ")" ++
  26         " (ymax: " ++ (show (ymax c)) ++ ")" ++
  27         " (zmin: " ++ (show (zmin c)) ++ ")" ++
  28         " (zmax: " ++ (show (zmax c)) ++ ")"
  29
  30 empty_cube :: Cube
  31 empty_cube = Cube 0 0 0 0 empty_values
  32
  33
  34 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
  35 --   p. 76.
  36 xmin :: Cube -> Double
  37 xmin c = (2*i' - 1)*delta / 2
  38     where
  39       i' = fromIntegral (i c) :: Double
  40       delta = h c
  41
  42 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
  43 --   p. 76.
  44 xmax :: Cube -> Double
  45 xmax c = (2*i' + 1)*delta / 2
  46     where
  47       i' = fromIntegral (i c) :: Double
  48       delta = h c
  49
  50 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
  51 --   p. 76.
  52 ymin :: Cube -> Double
  53 ymin c = (2*j' - 1)*delta / 2
  54     where
  55       j' = fromIntegral (j c) :: Double
  56       delta = h c
  57
  58 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
  59 --   p. 76.
  60 ymax :: Cube -> Double
  61 ymax c = (2*j' + 1)*delta / 2
  62     where
  63       j' = fromIntegral (j c) :: Double
  64       delta = h c
  65
  66 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
  67 --   p. 76.
  68 zmin :: Cube -> Double
  69 zmin c = (2*k' - 1)*delta / 2
  70     where
  71       k' = fromIntegral (k c) :: Double
  72       delta = h c
  73
  74 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
  75 --   p. 76.
  76 zmax :: Cube -> Double
  77 zmax c = (2*k' + 1)*delta / 2
  78     where
  79       k' = fromIntegral (k c) :: Double
  80       delta = h c
  81
  82 instance ThreeDimensional Cube where
  83     -- | The center of Cube_ijk coincides with v_ijk at
  84     --   (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
  85     center c = (x, y, z)
  86            where
  87              delta = h c
  88              i' = fromIntegral (i c) :: Double
  89              j' = fromIntegral (j c) :: Double
  90              k' = fromIntegral (k c) :: Double
  91              x = delta * i'
  92              y = delta * j'
  93              z = delta * k'
  94
  95     contains_point c p
  96         | (x_coord p) < (xmin c) = False
  97         | (x_coord p) > (xmax c) = False
  98         | (y_coord p) < (ymin c) = False
  99         | (y_coord p) > (ymax c) = False
 100         | (z_coord p) < (zmin c) = False
 101         | (z_coord p) > (zmax c) = False
 102         | otherwise = True
 103
 104
 105 -- instance Num Cube where
 106 --     (Cube g1 i1 j1 k1 d1) + (Cube _ i2 j2 k2 d2) =
 107 --         Cube g1 (i1 + i2) (j1 + j2) (k1 + k2) (d1 + d2)
 108
 109 --     (Cube g1 i1 j1 k1 d1) - (Cube _ i2 j2 k2 d2) =
 110 --         Cube g1 (i1 - i2) (j1 - j2) (k1 - k2) (d1 - d2)
 111
 112 --     (Cube g1 i1 j1 k1 d1) * (Cube _ i2 j2 k2 d2) =
 113 --         Cube g1 (i1 * i2) (j1 * j2) (k1 * k2) (d1 * d2)
 114
 115 --     abs (Cube g1 i1 j1 k1 d1) =
 116 --         Cube g1 (abs i1) (abs j1) (abs k1) (abs d1)
 117
 118 --     signum (Cube g1 i1 j1 k1 d1) =
 119 --         Cube g1 (signum i1) (signum j1) (signum k1) (signum d1)
 120
 121 --     fromInteger x = empty_cube { datum = (fromIntegral x) }
 122
 123 -- instance Fractional Cube where
 124 --     (Cube g1 i1 j1 k1 d1) / (Cube _ _ _ _ d2) =
 125 --         Cube g1 i1 j1 k1 (d1 / d2)
 126
 127 --     recip (Cube g1 i1 j1 k1 d1) =
 128 --         Cube g1 i1 j1 k1 (recip d1)
 129
 130 --     fromRational q = empty_cube { datum = fromRational q }
 131
 132
 133
 134 -- | Return the cube corresponding to the grid point i,j,k. The list
 135 --   of cubes is stored as [z][y][x] but we'll be requesting it by
 136 --   [x][y][z] so we flip the indices in the last line.
 137 -- cube_at :: Grid -> Int -> Int -> Int -> Cube
 138 -- cube_at g i' j' k'
 139 --         | i' >= length (function_values g) = Cube g i' j' k' 0
 140 --         | i' < 0                          = Cube g i' j' k' 0
 141 --         | j' >= length ((function_values g) !! i') = Cube g i' j' k' 0
 142 --         | j' < 0                                  = Cube g i' j' k' 0
 143 --         | k' >= length (((function_values g) !! i') !! j') = Cube g i' j' k' 0
 144 --         | k' < 0                                          = Cube g i' j' k' 0
 145 --         | otherwise =
 146 --             (((cubes g) !! k') !! j') !! i'
 147
 148
 149
 150
 151
 152
 153 -- Face stuff.
 154
 155 -- | The top (in the direction of z) face of the cube.
 156 top_face :: Cube -> Face
 157 top_face c = 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 back (in the direction of x) face of the cube.
 168 back_face :: Cube -> Face
 169 back_face c = 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
 178 -- The bottom face (in the direction of -z) of the cube.
 179 down_face :: Cube -> Face
 180 down_face c = 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
 190 -- | The front (in the direction of -x) face of the cube.
 191 front_face :: Cube -> Face
 192 front_face c = Face v0' v1' v2' v3'
 193     where
 194       delta = (1/2)*(h c)
 195       v0' = (center c) + (-delta, -delta, delta)
 196       v1' = (center c) + (-delta, delta, delta)
 197       v2' = (center c) + (-delta, delta, -delta)
 198       v3' = (center c) + (-delta, -delta, -delta)
 199
 200
 201 -- | The left (in the direction of -y) face of the cube.
 202 left_face :: Cube -> Face
 203 left_face c = Face v0' v1' v2' v3'
 204     where
 205       delta = (1/2)*(h c)
 206       v0' = (center c) + (-delta, -delta, delta)
 207       v1' = (center c) + (delta, -delta, delta)
 208       v2' = (center c) + (delta, -delta, -delta)
 209       v3' = (center c) + (-delta, -delta, -delta)
 210
 211
 212 -- | The right (in the direction of y) face of the cube.
 213 right_face :: Cube -> Face
 214 right_face c = Face v0' v1' v2' v3'
 215     where
 216       delta = (1/2)*(h c)
 217       v0' = (center c) + (-delta, delta, -delta)
 218       v1' = (center c) + (delta, delta, -delta)
 219       v2' = (center c) + (delta, delta, delta)
 220       v3' = (center c) + (-delta, delta, delta)
 221
 222
 223
 224 tetrahedron0 :: Cube -> Tetrahedron
 225 tetrahedron0 c =
 226     Tetrahedron (Cube.fv c) v0' v1' v2' v3'
 227     where
 228       v0' = center c
 229       v1' = center (front_face c)
 230       v2' = v0 (front_face c)
 231       v3' = v1 (front_face c)
 232
 233 tetrahedron1 :: Cube -> Tetrahedron
 234 tetrahedron1 c =
 235     Tetrahedron fv' v0' v1' v2' v3'
 236     where
 237       v0' = center c
 238       v1' = center (front_face c)
 239       v2' = v1 (front_face c)
 240       v3' = v2 (front_face c)
 241       fv' = rotate (Cube.fv c) ccwx
 242
 243 tetrahedron2 :: Cube -> Tetrahedron
 244 tetrahedron2 c =
 245     Tetrahedron fv' v0' v1' v2' v3'
 246     where
 247       v0' = center c
 248       v1' = center (front_face c)
 249       v2' = v2 (front_face c)
 250       v3' = v3 (front_face c)
 251       fv' = rotate (Cube.fv c) (ccwx . ccwx)
 252
 253 tetrahedron3 :: Cube -> Tetrahedron
 254 tetrahedron3 c =
 255     Tetrahedron fv' v0' v1' v2' v3'
 256     where
 257       v0' = center c
 258       v1' = center (front_face c)
 259       v2' = v3 (front_face c)
 260       v3' = v1 (front_face c)
 261       fv' = rotate (Cube.fv c) cwx
 262
 263 tetrahedron4 :: Cube -> Tetrahedron
 264 tetrahedron4 c =
 265     Tetrahedron fv' v0' v1' v2' v3'
 266     where
 267       v0' = center c
 268       v1' = center (top_face c)
 269       v2' = v0 (front_face c)
 270       v3' = v1 (front_face c)
 271       fv' = rotate (Cube.fv c) cwy
 272
 273 tetrahedron5 :: Cube -> Tetrahedron
 274 tetrahedron5 c =
 275     Tetrahedron fv' v0' v1' v2' v3'
 276     where
 277       v0' = center c
 278       v1' = center (top_face c)
 279       v2' = v1 (top_face c)
 280       v3' = v2 (top_face c)
 281       fv' = rotate (Tetrahedron.fv (tetrahedron0 c)) ccwx
 282
 283 tetrahedron6 :: Cube -> Tetrahedron
 284 tetrahedron6 c =
 285     Tetrahedron fv' v0' v1' v2' v3'
 286     where
 287       v0' = center c
 288       v1' = center (top_face c)
 289       v2' = v2 (top_face c)
 290       v3' = v3 (top_face c)
 291       fv' = rotate (Tetrahedron.fv (tetrahedron0 c)) (ccwx . ccwx)
 292
 293 tetrahedron7 :: Cube -> Tetrahedron
 294 tetrahedron7 c =
 295     Tetrahedron fv' v0' v1' v2' v3'
 296     where
 297       v0' = center c
 298       v1' = center (top_face c)
 299       v2' = v3 (top_face c)
 300       v3' = v1 (top_face c)
 301       fv' = rotate (Tetrahedron.fv (tetrahedron0 c)) cwx
 302
 303 tetrahedrons :: Cube -> [Tetrahedron]
 304 tetrahedrons c =
 305     [tetrahedron0 c,
 306      tetrahedron1 c,
 307      tetrahedron2 c,
 308      tetrahedron3 c,
 309      tetrahedron4 c,
 310      tetrahedron5 c,
 311      tetrahedron6 c,
 312      tetrahedron7 c
 313                 -- ,
 314                 --  tetrahedron8 c,
 315                 --  tetrahedron9 c,
 316                 --  tetrahedron10 c,
 317                 --  tetrahedron11 c,
 318                 --  tetrahedron12 c,
 319                 --  tetrahedron13 c,
 320                 --  tetrahedron14 c,
 321                 --  tetrahedron15 c,
 322                 --  tetrahedron16 c,
 323                 --  tetrahedron17 c,
 324                 --  tetrahedron18 c,
 325                 --  tetrahedron19 c,
 326                 --  tetrahedron20 c,
 327                 --  tetrahedron21 c,
 328                 --  tetrahedron21 c,
 329                 --  tetrahedron22 c,
 330                 --  tetrahedron23 c,
 331                 --  tetrahedron24 c
 332     ]