From d9eed953bd810f6928de536617dc21121a8a645b Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sun, 6 Nov 2011 21:32:55 -0500 Subject: [PATCH] Remove the 'h' parameter from the model entirely by defining h=1. --- src/Cube.hs | 55 ++++++++++++++++---------------------------- src/Everything.hs | 1 - src/Grid.hs | 58 +++++++++++++++++++++-------------------------- 3 files changed, 46 insertions(+), 68 deletions(-) diff --git a/src/Cube.hs b/src/Cube.hs index d863c29..6652e8b 100644 --- a/src/Cube.hs +++ b/src/Cube.hs @@ -30,8 +30,7 @@ import Misc (all_equal, disjoint) import Point (Point(..), dot) import Tetrahedron (Tetrahedron(..), barycenter, c, volume) -data Cube = Cube { h :: !Double, - i :: !Int, +data Cube = Cube { i :: !Int, j :: !Int, k :: !Int, fv :: !FunctionValues, @@ -41,13 +40,12 @@ data Cube = Cube { h :: !Double, instance Arbitrary Cube where arbitrary = do - (Positive h') <- arbitrary :: Gen (Positive Double) i' <- choose (coordmin, coordmax) j' <- choose (coordmin, coordmax) k' <- choose (coordmin, coordmax) fv' <- arbitrary :: Gen FunctionValues (Positive tet_vol) <- arbitrary :: Gen (Positive Double) - return (Cube h' i' j' k' fv' tet_vol) + return (Cube i' j' k' fv' tet_vol) where -- The idea here is that, when cubed in the volume formula, -- these numbers don't overflow 64 bits. This number is not @@ -60,7 +58,6 @@ instance Arbitrary Cube where instance Show Cube where show cube = "Cube_" ++ subscript ++ "\n" ++ - " h: " ++ (show (h cube)) ++ "\n" ++ " Center: " ++ (show (center cube)) ++ "\n" ++ " xmin: " ++ (show (xmin cube)) ++ "\n" ++ " xmax: " ++ (show (xmax cube)) ++ "\n" ++ @@ -76,65 +73,55 @@ instance Show Cube where -- | The left-side boundary of the cube. See Sorokina and Zeilfelder, -- p. 76. xmin :: Cube -> Double -xmin cube = (i' - 1/2)*delta +xmin cube = (i' - 1/2) where i' = fromIntegral (i cube) :: Double - delta = h cube -- | The right-side boundary of the cube. See Sorokina and Zeilfelder, -- p. 76. xmax :: Cube -> Double -xmax cube = (i' + 1/2)*delta +xmax cube = (i' + 1/2) where i' = fromIntegral (i cube) :: Double - delta = h cube -- | The front boundary of the cube. See Sorokina and Zeilfelder, -- p. 76. ymin :: Cube -> Double -ymin cube = (j' - 1/2)*delta +ymin cube = (j' - 1/2) where j' = fromIntegral (j cube) :: Double - delta = h cube -- | The back boundary of the cube. See Sorokina and Zeilfelder, -- p. 76. ymax :: Cube -> Double -ymax cube = (j' + 1/2)*delta +ymax cube = (j' + 1/2) where j' = fromIntegral (j cube) :: Double - delta = h cube -- | The bottom boundary of the cube. See Sorokina and Zeilfelder, -- p. 76. zmin :: Cube -> Double -zmin cube = (k' - 1/2)*delta +zmin cube = (k' - 1/2) where k' = fromIntegral (k cube) :: Double - delta = h cube -- | The top boundary of the cube. See Sorokina and Zeilfelder, -- p. 76. zmax :: Cube -> Double -zmax cube = (k' + 1/2)*delta +zmax cube = (k' + 1/2) where k' = fromIntegral (k cube) :: Double - delta = h cube -- | The center of Cube_ijk coincides with v_ijk at --- (ih, jh, kh). See Sorokina and Zeilfelder, p. 76. +-- (i, j, k). See Sorokina and Zeilfelder, p. 76. center :: Cube -> Point center cube = Point x y z where - delta = h cube - i' = fromIntegral (i cube) :: Double - j' = fromIntegral (j cube) :: Double - k' = fromIntegral (k cube) :: Double - x = delta * i' - y = delta * j' - z = delta * k' + x = fromIntegral (i cube) :: Double + y = fromIntegral (j cube) :: Double + z = fromIntegral (k cube) :: Double -- Face stuff. @@ -143,7 +130,7 @@ center cube = top_face :: Cube -> Face.Face top_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = 1/2 cc = center cube v0' = cc + ( Point delta (-delta) delta ) v1' = cc + ( Point delta delta delta ) @@ -156,7 +143,7 @@ top_face cube = Face.Face v0' v1' v2' v3' back_face :: Cube -> Face.Face back_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = 1/2 cc = center cube v0' = cc + ( Point delta (-delta) (-delta) ) v1' = cc + ( Point delta delta (-delta) ) @@ -168,7 +155,7 @@ back_face cube = Face.Face v0' v1' v2' v3' down_face :: Cube -> Face.Face down_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = 1/2 cc = center cube v0' = cc + ( Point (-delta) (-delta) (-delta) ) v1' = cc + ( Point (-delta) delta (-delta) ) @@ -181,7 +168,7 @@ down_face cube = Face.Face v0' v1' v2' v3' front_face :: Cube -> Face.Face front_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = 1/2 cc = center cube v0' = cc + ( Point (-delta) (-delta) delta ) v1' = cc + ( Point (-delta) delta delta ) @@ -192,7 +179,7 @@ front_face cube = Face.Face v0' v1' v2' v3' left_face :: Cube -> Face.Face left_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = 1/2 cc = center cube v0' = cc + ( Point delta (-delta) delta ) v1' = cc + ( Point (-delta) (-delta) delta ) @@ -204,7 +191,7 @@ left_face cube = Face.Face v0' v1' v2' v3' right_face :: Cube -> Face.Face right_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = 1/2 cc = center cube v0' = cc + ( Point (-delta) delta delta) v1' = cc + ( Point delta delta delta ) @@ -709,12 +696,10 @@ prop_all_volumes_positive cube = -- | In fact, since all of the tetrahedra are identical, we should -- already know their volumes. There's 24 tetrahedra to a cube, so --- we'd expect the volume of each one to be (1/24)*h^3. +-- we'd expect the volume of each one to be 1/24. prop_all_volumes_exact :: Cube -> Bool prop_all_volumes_exact cube = - and [volume t ~~= (1/24)*(delta^(3::Int)) | t <- tetrahedra cube] - where - delta = h cube + and [volume t ~~= 1/24 | t <- tetrahedra cube] -- | All tetrahedron should have their v0 located at the center of the -- cube. diff --git a/src/Everything.hs b/src/Everything.hs index af74566..5225c74 100644 --- a/src/Everything.hs +++ b/src/Everything.hs @@ -8,7 +8,6 @@ where import Cardinal as X import Comparisons as X -import Cube as X hiding (h) import Examples as X import Face as X hiding (v0,v1,v2,v3) import FunctionValues as X diff --git a/src/Grid.hs b/src/Grid.hs index 6627549..26f4425 100644 --- a/src/Grid.hs +++ b/src/Grid.hs @@ -1,7 +1,6 @@ {-# LANGUAGE BangPatterns #-} --- | The Grid module just contains the Grid type and two constructors --- for it. We hide the main Grid constructor because we don't want --- to allow instantiation of a grid with h <= 0. +-- | The Grid module contains the Grid type, its tests, and the 'zoom' +-- function used to build the interpolation. module Grid ( cube_at, grid_tests, @@ -18,7 +17,6 @@ import Test.Framework.Providers.QuickCheck2 (testProperty) import Test.QuickCheck ((==>), Arbitrary(..), Gen, - Positive(..), Property, choose) import Assertions (assertAlmostEqual, assertTrue) @@ -40,19 +38,18 @@ import Values (Values3D, dims, empty3d, zoom_shape) -- | Our problem is defined on a Grid. The grid size is given by the --- positive number h. The function values are the values of the --- function at the grid points, which are distance h from one --- another in each direction (x,y,z). -data Grid = Grid { h :: Double, -- MUST BE GREATER THAN ZERO! - function_values :: Values3D } +-- positive number h, which we have defined to always be 1 for +-- performance reasons (and simplicity). The function values are the +-- values of the function at the grid points, which are distance h=1 +-- from one another in each direction (x,y,z). +data Grid = Grid { function_values :: Values3D } deriving (Show) instance Arbitrary Grid where arbitrary = do - (Positive h') <- arbitrary :: Gen (Positive Double) fvs <- arbitrary :: Gen Values3D - return $ Grid h' fvs + return $ Grid fvs @@ -62,19 +59,18 @@ instance Arbitrary Grid where -- does improve performance. cube_at :: Grid -> Int -> Int -> Int -> Cube cube_at !g !i !j !k = - Cube delta i j k fvs' tet_vol + Cube i j k fvs' tet_vol where fvs = function_values g fvs' = make_values fvs i j k - delta = h g - tet_vol = (1/24)*(delta^(3::Int)) + tet_vol = 1/24 --- The first cube along any axis covers (-h/2, h/2). The second --- covers (h/2, 3h/2). The third, (3h/2, 5h/2), and so on. +-- The first cube along any axis covers (-1/2, 1/2). The second +-- covers (1/2, 3/2). The third, (3/2, 5/2), and so on. -- --- We translate the (x,y,z) coordinates forward by 'h/2' so that the --- first covers (0, h), the second covers (h, 2h), etc. This makes +-- We translate the (x,y,z) coordinates forward by 1/2 so that the +-- first covers (0, 1), the second covers (1, 2), etc. This makes -- it easy to figure out which cube contains the given point. calculate_containing_cube_coordinate :: Grid -> Double -> Int calculate_containing_cube_coordinate g coord @@ -83,11 +79,10 @@ calculate_containing_cube_coordinate g coord -- exists. | coord < offset = 0 | coord == offset && (xsize > 1 && ysize > 1 && zsize > 1) = 1 - | otherwise = (ceiling ( (coord + offset) / cube_width )) - 1 + | otherwise = (ceiling (coord + offset)) - 1 where (xsize, ysize, zsize) = dims (function_values g) - cube_width = (h g) - offset = cube_width / 2 + offset = 1/2 -- | Takes a 'Grid', and returns a 'Cube' containing the given 'Point'. @@ -111,8 +106,8 @@ zoom_result :: Values3D -> ScaleFactor -> R.DIM3 -> Double zoom_result v3d (sfx, sfy, sfz) (R.Z R.:. m R.:. n R.:. o) = f p where - g = Grid 1 v3d - offset = (h g)/2 + g = Grid v3d + offset = 1/2 m' = (fromIntegral m) / (fromIntegral sfx) - offset n' = (fromIntegral n) / (fromIntegral sfy) - offset o' = (fromIntegral o) / (fromIntegral sfz) - offset @@ -172,7 +167,7 @@ trilinear_c0_t0_tests = testCase "v3 is correct" test_trilinear_f0_t0_v3] ] where - g = Grid 1 trilinear + g = Grid trilinear cube = cube_at g 1 1 1 t = tetrahedron cube 0 @@ -287,7 +282,7 @@ test_trilinear_reproduced = let j' = fromIntegral j, let k' = fromIntegral k] where - g = Grid 1 trilinear + g = Grid trilinear cs = [ cube_at g ci cj ck | ci <- [0..2], cj <- [0..2], ck <- [0..2] ] @@ -305,7 +300,7 @@ test_zeros_reproduced = t0 <- tetrahedra c0, let p = polynomial t0 ] where - g = Grid 1 zeros + g = Grid zeros cs = [ cube_at g ci cj ck | ci <- [0..2], cj <- [0..2], ck <- [0..2] ] @@ -323,7 +318,7 @@ test_trilinear9x9x9_reproduced = let j' = (fromIntegral j) * 0.5, let k' = (fromIntegral k) * 0.5] where - g = Grid 1 trilinear + g = Grid trilinear c0 = cube_at g 1 1 1 @@ -331,13 +326,12 @@ test_trilinear9x9x9_reproduced = prop_cube_indices_never_go_out_of_bounds :: Grid -> Gen Bool prop_cube_indices_never_go_out_of_bounds g = do - let delta = Grid.h g - let coordmin = negate (delta/2) + let coordmin = negate (1/2) let (xsize, ysize, zsize) = dims $ function_values g - let xmax = delta*(fromIntegral xsize) - (delta/2) - let ymax = delta*(fromIntegral ysize) - (delta/2) - let zmax = delta*(fromIntegral zsize) - (delta/2) + let xmax = (fromIntegral xsize) - (1/2) + let ymax = (fromIntegral ysize) - (1/2) + let zmax = (fromIntegral zsize) - (1/2) x <- choose (coordmin, xmax) y <- choose (coordmin, ymax) -- 2.43.2