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