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                       c0 <- cs,
 295                       t <- tetrahedra c0,
 296                       let p = polynomial t,
 297                       let i' = fromIntegral i,
 298                       let j' = fromIntegral j,
 299                       let k' = fromIntegral k]
 300     where
 301       g = make_grid 1 trilinear
 302       cs = [ cube_at g ci cj ck | ci <- [0..2], cj <- [0..2], ck <- [0..2] ]
 303
 304
 305 test_zeros_reproduced :: Assertion
 306 test_zeros_reproduced =
 307     assertTrue "the zero function is reproduced correctly" $
 308              and [p (i', j', k') ~= value_at zeros i j k
 309                     | i <- [0..2],
 310                       j <- [0..2],
 311                       k <- [0..2],
 312                       let i' = fromIntegral i,
 313                       let j' = fromIntegral j,
 314                       let k' = fromIntegral k,
 315                       c0 <- cs,
 316                       t0 <- tetrahedra c0,
 317                       let p = polynomial t0 ]
 318     where
 319       g = make_grid 1 zeros
 320       cs = [ cube_at g ci cj ck | ci <- [0..2], cj <- [0..2], ck <- [0..2] ]
 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 ]