src/Grid.hs

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