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