src/Grid.hs

   1 -- | The Grid module just contains the Grid type and two constructors
   2 --   for it. We hide the main Grid constructor because we don't want
   3 --   to allow instantiation of a grid with h <= 0.
   4 module Grid (
   5   cube_at,
   6   grid_tests,
   7   make_grid,
   8   slow_tests,
   9   zoom
  10   )
  11 where
  12
  13 import Data.Array (Array, array, (!))
  14 import qualified Data.Array.Repa as R
  15 import Test.HUnit
  16 import Test.Framework (Test, testGroup)
  17 import Test.Framework.Providers.HUnit (testCase)
  18 import Test.Framework.Providers.QuickCheck2 (testProperty)
  19 import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
  20
  21 import Assertions
  22 import Comparisons
  23 import Cube (Cube(Cube),
  24              find_containing_tetrahedron,
  25              tetrahedra,
  26              tetrahedron)
  27 import Examples
  28 import FunctionValues
  29 import Point (Point)
  30 import ScaleFactor
  31 import Tetrahedron (Tetrahedron, c, polynomial, v0, v1, v2, v3)
  32 import ThreeDimensional
  33 import Values (Values3D, dims, empty3d, zoom_shape)
  34
  35
  36 type CubeGrid = Array (Int,Int,Int) Cube
  37
  38
  39 -- | Our problem is defined on a Grid. The grid size is given by the
  40 --   positive number h. The function values are the values of the
  41 --   function at the grid points, which are distance h from one
  42 --   another in each direction (x,y,z).
  43 data Grid = Grid { h :: Double, -- MUST BE GREATER THAN ZERO!
  44                    function_values :: Values3D,
  45                    cube_grid :: CubeGrid }
  46           deriving (Eq, Show)
  47
  48
  49 instance Arbitrary Grid where
  50     arbitrary = do
  51       (Positive h') <- arbitrary :: Gen (Positive Double)
  52       fvs <- arbitrary :: Gen Values3D
  53       return (make_grid h' fvs)
  54
  55
  56 -- | The constructor that we want people to use. If we're passed a
  57 --   non-positive grid size, we throw an error.
  58 make_grid :: Double -> Values3D -> Grid
  59 make_grid grid_size values
  60     | grid_size <= 0 = error "grid size must be positive"
  61     | otherwise = Grid grid_size values (cubes grid_size values)
  62
  63
  64 -- | Returns a three-dimensional array of cubes centered on the grid
  65 --   points (h*i, h*j, h*k) with the appropriate 'FunctionValues'.
  66 cubes :: Double -> Values3D -> CubeGrid
  67 cubes delta fvs
  68   = array (lbounds, ubounds)
  69            [ ((i,j,k), cube_ijk)
  70                  | i <- [0..xmax],
  71                    j <- [0..ymax],
  72                    k <- [0..zmax],
  73                    let tet_vol = (1/24)*(delta^(3::Int)),
  74                    let cube_ijk =
  75                          Cube delta i j k (make_values fvs i j k) tet_vol]
  76      where
  77        xmax = xsize - 1
  78        ymax = ysize - 1
  79        zmax = zsize - 1
  80        lbounds = (0, 0, 0)
  81        ubounds = (xmax, ymax, zmax)
  82        (xsize, ysize, zsize) = dims fvs
  83
  84
  85 -- | Takes a grid and a position as an argument and returns the cube
  86 --   centered on that position. If there is no cube there (i.e. the
  87 --   position is outside of the grid), it will throw an error.
  88 cube_at :: Grid -> Int -> Int -> Int -> Cube
  89 cube_at g i j k
  90     | i < 0      = error "i < 0 in cube_at"
  91     | i >= xsize = error "i >= xsize in cube_at"
  92     | j < 0      = error "j < 0 in cube_at"
  93     | j >= ysize = error "j >= ysize in cube_at"
  94     | k < 0      = error "k < 0 in cube_at"
  95     | k >= zsize = error "k >= zsize in cube_at"
  96     | otherwise = (cube_grid g) ! (i,j,k)
  97       where
  98         fvs = function_values g
  99         (xsize, ysize, zsize) = dims fvs
 100
 101 --   The first cube along any axis covers (-h/2, h/2). The second
 102 --   covers (h/2, 3h/2).  The third, (3h/2, 5h/2), and so on.
 103 --
 104 --   We translate the (x,y,z) coordinates forward by 'h/2' so that the
 105 --   first covers (0, h), the second covers (h, 2h), etc. This makes
 106 --   it easy to figure out which cube contains the given point.
 107 calculate_containing_cube_coordinate :: Grid -> Double -> Int
 108 calculate_containing_cube_coordinate g coord
 109     -- Don't use a cube on the boundary if we can help it. This
 110     -- returns cube #1 if we would have returned cube #0 and cube #1
 111     -- exists.
 112     | coord < offset = 0
 113     | coord == offset && (xsize > 1 && ysize > 1 && zsize > 1) = 1
 114     | otherwise = (ceiling ( (coord + offset) / cube_width )) - 1
 115     where
 116       (xsize, ysize, zsize) = dims (function_values g)
 117       cube_width = (h g)
 118       offset = cube_width / 2
 119
 120
 121 -- | Takes a 'Grid', and returns a 'Cube' containing the given 'Point'.
 122 --   Since our grid is rectangular, we can figure this out without having
 123 --   to check every cube.
 124 find_containing_cube :: Grid -> Point -> Cube
 125 find_containing_cube g p =
 126     cube_at g i j k
 127     where
 128       (x, y, z) = p
 129       i = calculate_containing_cube_coordinate g x
 130       j = calculate_containing_cube_coordinate g y
 131       k = calculate_containing_cube_coordinate g z
 132
 133
 134 {-# INLINE zoom_lookup #-}
 135 zoom_lookup :: Grid -> ScaleFactor -> a -> (R.DIM3 -> Double)
 136 zoom_lookup g scale_factor _ =
 137     zoom_result g scale_factor
 138
 139
 140 {-# INLINE zoom_result #-}
 141 zoom_result :: Grid -> ScaleFactor -> R.DIM3 -> Double
 142 zoom_result g (sfx, sfy, sfz) (R.Z R.:. m R.:. n R.:. o) =
 143   f p
 144   where
 145     offset = (h g)/2
 146     m' = (fromIntegral m) / (fromIntegral sfx) - offset
 147     n' = (fromIntegral n) / (fromIntegral sfy) - offset
 148     o' = (fromIntegral o) / (fromIntegral sfz) - offset
 149     p  = (m', n', o') :: Point
 150     cube = find_containing_cube g p
 151     t = find_containing_tetrahedron cube p
 152     f = polynomial t
 153
 154 zoom :: Grid -> ScaleFactor -> Values3D
 155 zoom g scale_factor
 156     | xsize == 0 || ysize == 0 || zsize == 0 = empty3d
 157     | otherwise =
 158         R.force $ R.unsafeTraverse arr transExtent (zoom_lookup g scale_factor)
 159           where
 160             arr = function_values g
 161             (xsize, ysize, zsize) = dims arr
 162             transExtent = zoom_shape scale_factor
 163
 164
 165
 166
 167 -- | Check all coefficients of tetrahedron0 belonging to the cube
 168 --   centered on (1,1,1) with a grid constructed from the trilinear
 169 --   values. See example one in the paper.
 170 --
 171 --   We also verify that the four vertices on face0 of the cube are
 172 --   in the correct location.
 173 --
 174 trilinear_c0_t0_tests :: Test.Framework.Test
 175 trilinear_c0_t0_tests =
 176   testGroup "trilinear c0 t0"
 177     [testGroup "coefficients"
 178       [testCase "c0030 is correct" test_trilinear_c0030,
 179        testCase "c0003 is correct" test_trilinear_c0003,
 180        testCase "c0021 is correct" test_trilinear_c0021,
 181        testCase "c0012 is correct" test_trilinear_c0012,
 182        testCase "c0120 is correct" test_trilinear_c0120,
 183        testCase "c0102 is correct" test_trilinear_c0102,
 184        testCase "c0111 is correct" test_trilinear_c0111,
 185        testCase "c0210 is correct" test_trilinear_c0210,
 186        testCase "c0201 is correct" test_trilinear_c0201,
 187        testCase "c0300 is correct" test_trilinear_c0300,
 188        testCase "c1020 is correct" test_trilinear_c1020,
 189        testCase "c1002 is correct" test_trilinear_c1002,
 190        testCase "c1011 is correct" test_trilinear_c1011,
 191        testCase "c1110 is correct" test_trilinear_c1110,
 192        testCase "c1101 is correct" test_trilinear_c1101,
 193        testCase "c1200 is correct" test_trilinear_c1200,
 194        testCase "c2010 is correct" test_trilinear_c2010,
 195        testCase "c2001 is correct" test_trilinear_c2001,
 196        testCase "c2100 is correct" test_trilinear_c2100,
 197        testCase "c3000 is correct" test_trilinear_c3000],
 198
 199     testGroup "face0 vertices"
 200       [testCase "v0 is correct" test_trilinear_f0_t0_v0,
 201        testCase "v1 is correct" test_trilinear_f0_t0_v1,
 202        testCase "v2 is correct" test_trilinear_f0_t0_v2,
 203        testCase "v3 is correct" test_trilinear_f0_t0_v3]
 204     ]
 205   where
 206     g = make_grid 1 trilinear
 207     cube = cube_at g 1 1 1
 208     t = tetrahedron cube 0
 209
 210     test_trilinear_c0030 :: Assertion
 211     test_trilinear_c0030 =
 212       assertAlmostEqual "c0030 is correct" (c t 0 0 3 0) (17/8)
 213
 214     test_trilinear_c0003 :: Assertion
 215     test_trilinear_c0003 =
 216       assertAlmostEqual "c0003 is correct" (c t 0 0 0 3) (27/8)
 217
 218     test_trilinear_c0021 :: Assertion
 219     test_trilinear_c0021 =
 220       assertAlmostEqual "c0021 is correct" (c t 0 0 2 1) (61/24)
 221
 222     test_trilinear_c0012 :: Assertion
 223     test_trilinear_c0012 =
 224       assertAlmostEqual "c0012 is correct" (c t 0 0 1 2) (71/24)
 225
 226     test_trilinear_c0120 :: Assertion
 227     test_trilinear_c0120 =
 228       assertAlmostEqual "c0120 is correct" (c t 0 1 2 0) (55/24)
 229
 230     test_trilinear_c0102 :: Assertion
 231     test_trilinear_c0102 =
 232       assertAlmostEqual "c0102 is correct" (c t 0 1 0 2) (73/24)
 233
 234     test_trilinear_c0111 :: Assertion
 235     test_trilinear_c0111 =
 236       assertAlmostEqual "c0111 is correct" (c t 0 1 1 1) (8/3)
 237
 238     test_trilinear_c0210 :: Assertion
 239     test_trilinear_c0210 =
 240       assertAlmostEqual "c0210 is correct" (c t 0 2 1 0) (29/12)
 241
 242     test_trilinear_c0201 :: Assertion
 243     test_trilinear_c0201 =
 244       assertAlmostEqual "c0201 is correct" (c t 0 2 0 1) (11/4)
 245
 246     test_trilinear_c0300 :: Assertion
 247     test_trilinear_c0300 =
 248       assertAlmostEqual "c0300 is correct" (c t 0 3 0 0) (5/2)
 249
 250     test_trilinear_c1020 :: Assertion
 251     test_trilinear_c1020 =
 252       assertAlmostEqual "c1020 is correct" (c t 1 0 2 0) (8/3)
 253
 254     test_trilinear_c1002 :: Assertion
 255     test_trilinear_c1002 =
 256       assertAlmostEqual "c1002 is correct" (c t 1 0 0 2) (23/6)
 257
 258     test_trilinear_c1011 :: Assertion
 259     test_trilinear_c1011 =
 260       assertAlmostEqual "c1011 is correct" (c t 1 0 1 1) (13/4)
 261
 262     test_trilinear_c1110 :: Assertion
 263     test_trilinear_c1110 =
 264       assertAlmostEqual "c1110 is correct" (c t 1 1 1 0) (23/8)
 265
 266     test_trilinear_c1101 :: Assertion
 267     test_trilinear_c1101 =
 268       assertAlmostEqual "c1101 is correct" (c t 1 1 0 1) (27/8)
 269
 270     test_trilinear_c1200 :: Assertion
 271     test_trilinear_c1200 =
 272       assertAlmostEqual "c1200 is correct" (c t 1 2 0 0) 3
 273
 274     test_trilinear_c2010 :: Assertion
 275     test_trilinear_c2010 =
 276       assertAlmostEqual "c2010 is correct" (c t 2 0 1 0) (10/3)
 277
 278     test_trilinear_c2001 :: Assertion
 279     test_trilinear_c2001 =
 280       assertAlmostEqual "c2001 is correct" (c t 2 0 0 1) 4
 281
 282     test_trilinear_c2100 :: Assertion
 283     test_trilinear_c2100 =
 284       assertAlmostEqual "c2100 is correct" (c t 2 1 0 0) (7/2)
 285
 286     test_trilinear_c3000 :: Assertion
 287     test_trilinear_c3000 =
 288       assertAlmostEqual "c3000 is correct" (c t 3 0 0 0) 4
 289
 290     test_trilinear_f0_t0_v0 :: Assertion
 291     test_trilinear_f0_t0_v0 =
 292       assertEqual "v0 is correct" (v0 t) (1, 1, 1)
 293
 294     test_trilinear_f0_t0_v1 :: Assertion
 295     test_trilinear_f0_t0_v1 =
 296       assertEqual "v1 is correct" (v1 t) (0.5, 1, 1)
 297
 298     test_trilinear_f0_t0_v2 :: Assertion
 299     test_trilinear_f0_t0_v2 =
 300       assertEqual "v2 is correct" (v2 t) (0.5, 0.5, 1.5)
 301
 302     test_trilinear_f0_t0_v3 :: Assertion
 303     test_trilinear_f0_t0_v3 =
 304       assertClose "v3 is correct" (v3 t) (0.5, 1.5, 1.5)
 305
 306
 307 test_trilinear_reproduced :: Assertion
 308 test_trilinear_reproduced =
 309     assertTrue "trilinears are reproduced correctly" $
 310              and [p (i', j', k') ~= value_at trilinear i j k
 311                     | i <- [0..2],
 312                       j <- [0..2],
 313                       k <- [0..2],
 314                       t <- tetrahedra c0,
 315                       let p = polynomial t,
 316                       let i' = fromIntegral i,
 317                       let j' = fromIntegral j,
 318                       let k' = fromIntegral k]
 319     where
 320       g = make_grid 1 trilinear
 321       c0 = cube_at g 1 1 1
 322
 323
 324 test_zeros_reproduced :: Assertion
 325 test_zeros_reproduced =
 326     assertTrue "the zero function is reproduced correctly" $
 327              and [p (i', j', k') ~= value_at zeros i j k
 328                     | i <- [0..2],
 329                       j <- [0..2],
 330                       k <- [0..2],
 331                       let i' = fromIntegral i,
 332                       let j' = fromIntegral j,
 333                       let k' = fromIntegral k]
 334     where
 335       g = make_grid 1 zeros
 336       c0 = cube_at g 1 1 1
 337       t0 = tetrahedron c0 0
 338       p = polynomial t0
 339
 340
 341 -- | Make sure we can reproduce a 9x9x9 trilinear from the 3x3x3 one.
 342 test_trilinear9x9x9_reproduced :: Assertion
 343 test_trilinear9x9x9_reproduced =
 344     assertTrue "trilinear 9x9x9 is reproduced correctly" $
 345       and [p (i', j', k') ~= value_at trilinear9x9x9 i j k
 346             | i <- [0..8],
 347               j <- [0..8],
 348               k <- [0..8],
 349               t <- tetrahedra c0,
 350               let p = polynomial t,
 351               let i' = (fromIntegral i) * 0.5,
 352               let j' = (fromIntegral j) * 0.5,
 353               let k' = (fromIntegral k) * 0.5]
 354     where
 355       g = make_grid 1 trilinear
 356       c0 = cube_at g 1 1 1
 357
 358
 359 -- | The point 'p' in this test lies on the boundary of tetrahedra 12 and 15.
 360 --   However, the 'contains_point' test fails due to some numerical innacuracy.
 361 --   This bug should have been fixed by setting a positive tolerance level.
 362 --
 363 --   Example from before the fix:
 364 --
 365 --   b1 (tetrahedron c 20) (0, 17.5, 0.5)
 366 --   -0.0
 367 --
 368 test_tetrahedra_collision_sensitivity :: Assertion
 369 test_tetrahedra_collision_sensitivity =
 370   assertTrue "tetrahedron collision tests isn't too sensitive" $
 371              contains_point t20 p
 372   where
 373     g = make_grid 1 naturals_1d
 374     cube = cube_at g 0 18 0
 375     p = (0, 17.5, 0.5) :: Point
 376     t20 = tetrahedron cube 20
 377
 378
 379 prop_cube_indices_never_go_out_of_bounds :: Grid -> Gen Bool
 380 prop_cube_indices_never_go_out_of_bounds g =
 381   do
 382     let delta = Grid.h g
 383     let coordmin = negate (delta/2)
 384
 385     let (xsize, ysize, zsize) = dims $ function_values g
 386     let xmax = delta*(fromIntegral xsize) - (delta/2)
 387     let ymax = delta*(fromIntegral ysize) - (delta/2)
 388     let zmax = delta*(fromIntegral zsize) - (delta/2)
 389
 390     x <- choose (coordmin, xmax)
 391     y <- choose (coordmin, ymax)
 392     z <- choose (coordmin, zmax)
 393
 394     let idx_x = calculate_containing_cube_coordinate g x
 395     let idx_y = calculate_containing_cube_coordinate g y
 396     let idx_z = calculate_containing_cube_coordinate g z
 397
 398     return $
 399       idx_x >= 0 &&
 400       idx_x <= xsize - 1 &&
 401       idx_y >= 0 &&
 402       idx_y <= ysize - 1 &&
 403       idx_z >= 0 &&
 404       idx_z <= zsize - 1
 405
 406
 407
 408 grid_tests :: Test.Framework.Test
 409 grid_tests =
 410     testGroup "Grid Tests" [
 411       trilinear_c0_t0_tests,
 412       testCase "tetrahedra collision test isn't too sensitive"
 413          test_tetrahedra_collision_sensitivity,
 414       testCase "trilinear reproduced" test_trilinear_reproduced,
 415       testCase "zeros reproduced" test_zeros_reproduced ]
 416
 417
 418 -- Do the slow tests last so we can stop paying attention.
 419 slow_tests :: Test.Framework.Test
 420 slow_tests =
 421     testGroup "Slow Tests" [
 422       testProperty "cube indices within bounds"
 423                    prop_cube_indices_never_go_out_of_bounds,
 424       testCase "trilinear9x9x9 reproduced" test_trilinear9x9x9_reproduced ]