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