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