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 PolynomialArray (PolynomialArray)
  31 import ScaleFactor
  32 import Tetrahedron (Tetrahedron, c, number, 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 cube_ijk =
  76                          Cube delta i j k (make_values fvs i j k) 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 -> PolynomialArray -> ScaleFactor -> a -> (R.DIM3 -> Double)
 137 zoom_lookup g polynomials scale_factor _ =
 138     zoom_result g polynomials scale_factor
 139
 140
 141 {-# INLINE zoom_result #-}
 142 zoom_result :: Grid -> PolynomialArray -> ScaleFactor -> R.DIM3 -> Double
 143 zoom_result g polynomials (sfx, sfy, sfz) (R.Z R.:. m R.:. n R.:. o) =
 144   (polynomials ! (i, j, k, (number t))) 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     -- Figure out i,j,k without importing those functions.
 153     Cube _ i j k _ _ = cube
 154     t = find_containing_tetrahedron cube p
 155
 156
 157 zoom :: Grid -> PolynomialArray -> ScaleFactor -> Values3D
 158 zoom g polynomials scale_factor
 159     | xsize == 0 || ysize == 0 || zsize == 0 = empty3d
 160     | otherwise =
 161         R.force $ R.traverse arr transExtent (zoom_lookup g polynomials 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 -- | Check all coefficients of tetrahedron0 belonging to the cube
 171 --   centered on (1,1,1) with a grid constructed from the trilinear
 172 --   values. See example one in the paper.
 173 --
 174 --   We also verify that the four vertices on face0 of the cube are
 175 --   in the correct location.
 176 --
 177 trilinear_c0_t0_tests :: Test.Framework.Test
 178 trilinear_c0_t0_tests =
 179   testGroup "trilinear c0 t0"
 180     [testGroup "coefficients"
 181       [testCase "c0030 is correct" test_trilinear_c0030,
 182        testCase "c0003 is correct" test_trilinear_c0003,
 183        testCase "c0021 is correct" test_trilinear_c0021,
 184        testCase "c0012 is correct" test_trilinear_c0012,
 185        testCase "c0120 is correct" test_trilinear_c0120,
 186        testCase "c0102 is correct" test_trilinear_c0102,
 187        testCase "c0111 is correct" test_trilinear_c0111,
 188        testCase "c0210 is correct" test_trilinear_c0210,
 189        testCase "c0201 is correct" test_trilinear_c0201,
 190        testCase "c0300 is correct" test_trilinear_c0300,
 191        testCase "c1020 is correct" test_trilinear_c1020,
 192        testCase "c1002 is correct" test_trilinear_c1002,
 193        testCase "c1011 is correct" test_trilinear_c1011,
 194        testCase "c1110 is correct" test_trilinear_c1110,
 195        testCase "c1101 is correct" test_trilinear_c1101,
 196        testCase "c1200 is correct" test_trilinear_c1200,
 197        testCase "c2010 is correct" test_trilinear_c2010,
 198        testCase "c2001 is correct" test_trilinear_c2001,
 199        testCase "c2100 is correct" test_trilinear_c2100,
 200        testCase "c3000 is correct" test_trilinear_c3000],
 201
 202     testGroup "face0 vertices"
 203       [testCase "v0 is correct" test_trilinear_f0_t0_v0,
 204        testCase "v1 is correct" test_trilinear_f0_t0_v1,
 205        testCase "v2 is correct" test_trilinear_f0_t0_v2,
 206        testCase "v3 is correct" test_trilinear_f0_t0_v3]
 207     ]
 208   where
 209     g = make_grid 1 trilinear
 210     cube = cube_at g 1 1 1
 211     t = tetrahedron cube 0
 212
 213     test_trilinear_c0030 :: Assertion
 214     test_trilinear_c0030 =
 215       assertAlmostEqual "c0030 is correct" (c t 0 0 3 0) (17/8)
 216
 217     test_trilinear_c0003 :: Assertion
 218     test_trilinear_c0003 =
 219       assertAlmostEqual "c0003 is correct" (c t 0 0 0 3) (27/8)
 220
 221     test_trilinear_c0021 :: Assertion
 222     test_trilinear_c0021 =
 223       assertAlmostEqual "c0021 is correct" (c t 0 0 2 1) (61/24)
 224
 225     test_trilinear_c0012 :: Assertion
 226     test_trilinear_c0012 =
 227       assertAlmostEqual "c0012 is correct" (c t 0 0 1 2) (71/24)
 228
 229     test_trilinear_c0120 :: Assertion
 230     test_trilinear_c0120 =
 231       assertAlmostEqual "c0120 is correct" (c t 0 1 2 0) (55/24)
 232
 233     test_trilinear_c0102 :: Assertion
 234     test_trilinear_c0102 =
 235       assertAlmostEqual "c0102 is correct" (c t 0 1 0 2) (73/24)
 236
 237     test_trilinear_c0111 :: Assertion
 238     test_trilinear_c0111 =
 239       assertAlmostEqual "c0111 is correct" (c t 0 1 1 1) (8/3)
 240
 241     test_trilinear_c0210 :: Assertion
 242     test_trilinear_c0210 =
 243       assertAlmostEqual "c0210 is correct" (c t 0 2 1 0) (29/12)
 244
 245     test_trilinear_c0201 :: Assertion
 246     test_trilinear_c0201 =
 247       assertAlmostEqual "c0201 is correct" (c t 0 2 0 1) (11/4)
 248
 249     test_trilinear_c0300 :: Assertion
 250     test_trilinear_c0300 =
 251       assertAlmostEqual "c0300 is correct" (c t 0 3 0 0) (5/2)
 252
 253     test_trilinear_c1020 :: Assertion
 254     test_trilinear_c1020 =
 255       assertAlmostEqual "c1020 is correct" (c t 1 0 2 0) (8/3)
 256
 257     test_trilinear_c1002 :: Assertion
 258     test_trilinear_c1002 =
 259       assertAlmostEqual "c1002 is correct" (c t 1 0 0 2) (23/6)
 260
 261     test_trilinear_c1011 :: Assertion
 262     test_trilinear_c1011 =
 263       assertAlmostEqual "c1011 is correct" (c t 1 0 1 1) (13/4)
 264
 265     test_trilinear_c1110 :: Assertion
 266     test_trilinear_c1110 =
 267       assertAlmostEqual "c1110 is correct" (c t 1 1 1 0) (23/8)
 268
 269     test_trilinear_c1101 :: Assertion
 270     test_trilinear_c1101 =
 271       assertAlmostEqual "c1101 is correct" (c t 1 1 0 1) (27/8)
 272
 273     test_trilinear_c1200 :: Assertion
 274     test_trilinear_c1200 =
 275       assertAlmostEqual "c1200 is correct" (c t 1 2 0 0) 3
 276
 277     test_trilinear_c2010 :: Assertion
 278     test_trilinear_c2010 =
 279       assertAlmostEqual "c2010 is correct" (c t 2 0 1 0) (10/3)
 280
 281     test_trilinear_c2001 :: Assertion
 282     test_trilinear_c2001 =
 283       assertAlmostEqual "c2001 is correct" (c t 2 0 0 1) 4
 284
 285     test_trilinear_c2100 :: Assertion
 286     test_trilinear_c2100 =
 287       assertAlmostEqual "c2100 is correct" (c t 2 1 0 0) (7/2)
 288
 289     test_trilinear_c3000 :: Assertion
 290     test_trilinear_c3000 =
 291       assertAlmostEqual "c3000 is correct" (c t 3 0 0 0) 4
 292
 293     test_trilinear_f0_t0_v0 :: Assertion
 294     test_trilinear_f0_t0_v0 =
 295       assertEqual "v0 is correct" (v0 t) (1, 1, 1)
 296
 297     test_trilinear_f0_t0_v1 :: Assertion
 298     test_trilinear_f0_t0_v1 =
 299       assertEqual "v1 is correct" (v1 t) (0.5, 1, 1)
 300
 301     test_trilinear_f0_t0_v2 :: Assertion
 302     test_trilinear_f0_t0_v2 =
 303       assertEqual "v2 is correct" (v2 t) (0.5, 0.5, 1.5)
 304
 305     test_trilinear_f0_t0_v3 :: Assertion
 306     test_trilinear_f0_t0_v3 =
 307       assertClose "v3 is correct" (v3 t) (0.5, 1.5, 1.5)
 308
 309
 310 test_trilinear_reproduced :: Assertion
 311 test_trilinear_reproduced =
 312     assertTrue "trilinears are reproduced correctly" $
 313              and [p (i', j', k') ~= value_at trilinear i j k
 314                     | i <- [0..2],
 315                       j <- [0..2],
 316                       k <- [0..2],
 317                       t <- tetrahedra c0,
 318                       let p = polynomial t,
 319                       let i' = fromIntegral i,
 320                       let j' = fromIntegral j,
 321                       let k' = fromIntegral k]
 322     where
 323       g = make_grid 1 trilinear
 324       c0 = cube_at g 1 1 1
 325
 326
 327 test_zeros_reproduced :: Assertion
 328 test_zeros_reproduced =
 329     assertTrue "the zero function is reproduced correctly" $
 330              and [p (i', j', k') ~= value_at zeros i j k
 331                     | i <- [0..2],
 332                       j <- [0..2],
 333                       k <- [0..2],
 334                       let i' = fromIntegral i,
 335                       let j' = fromIntegral j,
 336                       let k' = fromIntegral k]
 337     where
 338       g = make_grid 1 zeros
 339       c0 = cube_at g 1 1 1
 340       t0 = tetrahedron c0 0
 341       p = polynomial t0
 342
 343
 344 -- | Make sure we can reproduce a 9x9x9 trilinear from the 3x3x3 one.
 345 test_trilinear9x9x9_reproduced :: Assertion
 346 test_trilinear9x9x9_reproduced =
 347     assertTrue "trilinear 9x9x9 is reproduced correctly" $
 348       and [p (i', j', k') ~= value_at trilinear9x9x9 i j k
 349             | i <- [0..8],
 350               j <- [0..8],
 351               k <- [0..8],
 352               t <- tetrahedra c0,
 353               let p = polynomial t,
 354               let i' = (fromIntegral i) * 0.5,
 355               let j' = (fromIntegral j) * 0.5,
 356               let k' = (fromIntegral k) * 0.5]
 357     where
 358       g = make_grid 1 trilinear
 359       c0 = cube_at g 1 1 1
 360
 361
 362 -- | The point 'p' in this test lies on the boundary of tetrahedra 12 and 15.
 363 --   However, the 'contains_point' test fails due to some numerical innacuracy.
 364 --   This bug should have been fixed by setting a positive tolerance level.
 365 --
 366 --   Example from before the fix:
 367 --
 368 --   > b0 (tetrahedron c 15) p
 369 --   -3.4694469519536365e-18
 370 --
 371 test_tetrahedra_collision_sensitivity :: Assertion
 372 test_tetrahedra_collision_sensitivity =
 373   assertTrue "tetrahedron collision tests isn't too sensitive" $
 374              contains_point t15 p
 375   where
 376     g = make_grid 1 naturals_1d
 377     cube = cube_at g 0 17 1
 378     p = (0, 16.75, 0.5) :: Point
 379     t15 = tetrahedron cube 15
 380
 381
 382 prop_cube_indices_never_go_out_of_bounds :: Grid -> Gen Bool
 383 prop_cube_indices_never_go_out_of_bounds g =
 384   do
 385     let delta = Grid.h g
 386     let coordmin = negate (delta/2)
 387
 388     let (xsize, ysize, zsize) = dims $ function_values g
 389     let xmax = delta*(fromIntegral xsize) - (delta/2)
 390     let ymax = delta*(fromIntegral ysize) - (delta/2)
 391     let zmax = delta*(fromIntegral zsize) - (delta/2)
 392
 393     x <- choose (coordmin, xmax)
 394     y <- choose (coordmin, ymax)
 395     z <- choose (coordmin, zmax)
 396
 397     let idx_x = calculate_containing_cube_coordinate g x
 398     let idx_y = calculate_containing_cube_coordinate g y
 399     let idx_z = calculate_containing_cube_coordinate g z
 400
 401     return $
 402       idx_x >= 0 &&
 403       idx_x <= xsize - 1 &&
 404       idx_y >= 0 &&
 405       idx_y <= ysize - 1 &&
 406       idx_z >= 0 &&
 407       idx_z <= zsize - 1
 408
 409
 410
 411 grid_tests :: Test.Framework.Test
 412 grid_tests =
 413     testGroup "Grid Tests" [
 414       trilinear_c0_t0_tests,
 415       testCase "tetrahedra collision test isn't too sensitive"
 416          test_tetrahedra_collision_sensitivity,
 417       testCase "trilinear reproduced" test_trilinear_reproduced,
 418       testCase "zeros reproduced" test_zeros_reproduced ]
 419
 420
 421 -- Do the slow tests last so we can stop paying attention.
 422 slow_tests :: Test.Framework.Test
 423 slow_tests =
 424     testGroup "Slow Tests" [
 425       testProperty "cube indices within bounds"
 426                    prop_cube_indices_never_go_out_of_bounds,
 427       testCase "trilinear9x9x9 reproduced" test_trilinear9x9x9_reproduced ]