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