--- /dev/null
+[^.]*
\ No newline at end of file
--- /dev/null
+GHC_WARNINGS := -Wall
+GHC_WARNINGS += -fwarn-hi-shadowing
+GHC_WARNINGS += -fwarn-missing-signatures
+GHC_WARNINGS += -fwarn-name-shadowing
+GHC_WARNINGS += -fwarn-orphans
+GHC_WARNINGS += -fwarn-type-defaults
+
+BIN := spline
+
+.PHONY : test doc src_html
+
+$(BIN): src/*.hs
+ ghc -O2 $(GHC_WARNINGS) --make -o bin/${BIN} src/*.hs
+
+all: $(BIN) test_src
+
+test_src: src/Tests/*.hs
+ ghc -O2 $(GHC_WARNINGS) --make -o bin/${BIN} src/*.hs src/Tests/*.hs
+
+profile: src/*.hs
+ ghc -O2 $(GHC_WARNINGS) -prof -auto-all --make -o bin/$(BIN) src/*.hs
+
+clean:
+ rm -f bin/$(BIN)
+ rm -f src/*.hi
+ rm -f src/*.o
+ rm -f src/Tests/*.hi
+ rm -f src/Tests/*.o
+ rm -f *.prof
+ rm -rf doc/*
+
+test:
+ runghc -i"src" test/TestSuite.hs
+
+
+src_html:
+ util/hscolour_srcs
+
+
+DOCS := -i /usr/share/doc/ghc-6.12.3/html/libraries/base-4.2.0.2/base.haddock
+DOCS += -i /usr/share/doc/storable-complex-0.2.1/html/storable-complex.haddock
+DOCS += -i /usr/share/doc/hmatrix-0.10.0.1/html/hmatrix.haddock
+
+doc: src_html
+ haddock $(DOCS) --html --use-unicode \
+ --odir=doc/ --title="spline3" \
+ --source-module="src/%{MODULE/.//}.html" \
+ --source-entity="src/%{MODULE/.//}.html#%{NAME}" \
+ src/*.hs src/Tests/*.hs
--- /dev/null
+module Cube
+where
+
+import Grid
+import Point
+import ThreeDimensional
+
+class Gridded a where
+ back :: a -> Cube
+ down :: a -> Cube
+ front :: a -> Cube
+ left :: a -> Cube
+ right :: a -> Cube
+ top :: a -> Cube
+
+
+data Cube = Cube { grid :: Grid,
+ i :: Int,
+ j :: Int,
+ k :: Int,
+ datum :: Double }
+ deriving (Eq)
+
+
+instance Show Cube where
+ show c =
+ "Cube_" ++ (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k c)) ++
+ " (Grid: " ++ (show (grid c)) ++ ")" ++
+ " (Center: " ++ (show (center c)) ++ ")" ++
+ " (xmin: " ++ (show (xmin c)) ++ ")" ++
+ " (xmax: " ++ (show (xmax c)) ++ ")" ++
+ " (ymin: " ++ (show (ymin c)) ++ ")" ++
+ " (ymax: " ++ (show (ymax c)) ++ ")" ++
+ " (zmin: " ++ (show (zmin c)) ++ ")" ++
+ " (zmax: " ++ (show (zmax c)) ++ ")" ++
+ " (datum: " ++ (show (datum c)) ++ ")\n\n"
+
+empty_cube :: Cube
+empty_cube = Cube empty_grid 0 0 0 0
+
+-- TODO: The paper considers 'i' to be the front/back direction,
+-- whereas I have it in the left/right direction.
+instance Gridded Cube where
+ back c = cube_at (grid c) ((i c) + 1) (j c) (k c)
+ down c = cube_at (grid c) (i c) (j c) ((k c) - 1)
+ front c = cube_at (grid c) ((i c) - 1) (j c) (k c)
+ left c = cube_at (grid c) (i c) ((j c) - 1) (k c)
+ right c = cube_at (grid c) (i c) ((j c) + 1) (k c)
+ top c = cube_at (grid c) (i c) (j c) ((k c) + 1)
+
+-- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
+-- p. 76.
+xmin :: Cube -> Double
+xmin c = (2*i' - 1)*delta / 2
+ where
+ i' = fromIntegral (i c) :: Double
+ delta = h (grid c)
+
+-- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
+-- p. 76.
+xmax :: Cube -> Double
+xmax c = (2*i' + 1)*delta / 2
+ where
+ i' = fromIntegral (i c) :: Double
+ delta = h (grid c)
+
+-- | The front boundary of the cube. See Sorokina and Zeilfelder,
+-- p. 76.
+ymin :: Cube -> Double
+ymin c = (2*j' - 1)*delta / 2
+ where
+ j' = fromIntegral (j c) :: Double
+ delta = h (grid c)
+
+-- | The back boundary of the cube. See Sorokina and Zeilfelder,
+-- p. 76.
+ymax :: Cube -> Double
+ymax c = (2*j' + 1)*delta / 2
+ where
+ j' = fromIntegral (j c) :: Double
+ delta = h (grid c)
+
+-- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
+-- p. 76.
+zmin :: Cube -> Double
+zmin c = (2*k' - 1)*delta / 2
+ where
+ k' = fromIntegral (k c) :: Double
+ delta = h (grid c)
+
+-- | The top boundary of the cube. See Sorokina and Zeilfelder,
+-- p. 76.
+zmax :: Cube -> Double
+zmax c = (2*k' + 1)*delta / 2
+ where
+ k' = fromIntegral (k c) :: Double
+ delta = h (grid c)
+
+instance ThreeDimensional Cube where
+ -- | The center of Cube_ijk coincides with v_ijk at
+ -- (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
+ center c = (x, y, z)
+ where
+ delta = h (grid c)
+ i' = fromIntegral (i c) :: Double
+ j' = fromIntegral (j c) :: Double
+ k' = fromIntegral (k c) :: Double
+ x = (delta * i')
+ y = (delta * j')
+ z = (delta * k')
+
+ contains_point c p
+ | (x_coord p) < (xmin c) = False
+ | (x_coord p) > (xmax c) = False
+ | (y_coord p) < (ymin c) = False
+ | (y_coord p) > (ymax c) = False
+ | (z_coord p) < (zmin c) = False
+ | (z_coord p) > (zmax c) = False
+ | otherwise = True
+
+
+instance Num Cube where
+ (Cube g1 i1 j1 k1 d1) + (Cube _ i2 j2 k2 d2) =
+ Cube g1 (i1 + i2) (j1 + j2) (k1 + k2) (d1 + d2)
+
+ (Cube g1 i1 j1 k1 d1) - (Cube _ i2 j2 k2 d2) =
+ Cube g1 (i1 - i2) (j1 - j2) (k1 - k2) (d1 - d2)
+
+ (Cube g1 i1 j1 k1 d1) * (Cube _ i2 j2 k2 d2) =
+ Cube g1 (i1 * i2) (j1 * j2) (k1 * k2) (d1 * d2)
+
+ abs (Cube g1 i1 j1 k1 d1) =
+ Cube g1 (abs i1) (abs j1) (abs k1) (abs d1)
+
+ signum (Cube g1 i1 j1 k1 d1) =
+ Cube g1 (signum i1) (signum j1) (signum k1) (signum d1)
+
+ fromInteger x = empty_cube { datum = (fromIntegral x) }
+
+instance Fractional Cube where
+ (Cube g1 i1 j1 k1 d1) / (Cube _ _ _ _ d2) =
+ Cube g1 i1 j1 k1 (d1 / d2)
+
+ recip (Cube g1 i1 j1 k1 d1) =
+ Cube g1 i1 j1 k1 (recip d1)
+
+ fromRational q = empty_cube { datum = fromRational q }
+
+reverse_cube :: Grid -> Int -> Int -> Int -> Double -> Cube
+reverse_cube g k' j' i' datum' = Cube g i' j' k' datum'
+
+
+cube_at :: Grid -> Int -> Int -> Int -> Cube
+cube_at g i' j' k'
+ | i' >= length (function_values g) = Cube g i' j' k' 0
+ | i' < 0 = Cube g i' j' k' 0
+ | j' >= length ((function_values g) !! i') = Cube g i' j' k' 0
+ | j' < 0 = Cube g i' j' k' 0
+ | k' >= length (((function_values g) !! i') !! j') = Cube g i' j' k' 0
+ | k' < 0 = Cube g i' j' k' 0
+ | otherwise =
+ Cube g i' j' k' ((((function_values g) !! i') !! j') !! k')
+
+-- These next three functions basically form a 'for' loop, looping
+-- through the xs, ys, and zs in that order.
+cubes_from_values :: Grid -> Int -> Int -> [Double] -> [Cube]
+cubes_from_values g i' j' vals =
+ zipWith (reverse_cube g i' j') [0..] vals
+
+cubes_from_planes :: Grid -> Int -> [[Double]] -> [[Cube]]
+cubes_from_planes g i' planes =
+ zipWith (cubes_from_values g i') [0..] planes
+
+cubes :: Grid -> [[[Cube]]]
+cubes g = zipWith (cubes_from_planes g) [0..] (function_values g)
--- /dev/null
+module Face
+where
+
+import Cube
+import Grid
+import Point
+import Tetrahedron hiding (c, cube, v0, v1, v2, v3)
+import ThreeDimensional
+
+data Face = Face { cube :: Cube,
+ v0 :: Point,
+ v1 :: Point,
+ v2 :: Point,
+ v3 :: Point }
+ deriving (Eq)
+
+instance Show Face where
+ show f = "Face (Cube_" ++ (show i') ++ "," ++ (show j') ++ "," ++
+ (show k') ++ ") " ++ "(v0: " ++ (show (v0 f)) ++ ") (v1: " ++
+ (show (v1 f)) ++ ") (v2: " ++ (show (v2 f)) ++ ") (v3: " ++
+ (show (v3 f)) ++ ")\n\n"
+ where
+ i' = i (cube f)
+ j' = j (cube f)
+ k' = k (cube f)
+
+instance ThreeDimensional Face where
+ center f = ((v0 f) + (v1 f) + (v2 f) + (v3 f)) `scale` (1/4)
+ -- Too lazy to implement this right now.
+ contains_point _ _ = False
+
+-- The top face of the cube.
+face0 :: Cube -> Face
+face0 c = Face c v0' v1' v2' v3'
+ where
+ g = grid c
+ delta = (1/2)*(h g)
+ v0' = (center c) + (-delta, delta, delta)
+ v1' = (center c) + (delta, delta, delta)
+ v2' = (center c) + (delta, -delta, delta)
+ v3' = (center c) + (-delta, -delta, delta)
+
+-- The right face of the cube.
+face1 :: Cube -> Face
+face1 c = Face c v0' v1' v2' v3'
+ where
+ g = grid c
+ delta = (1/2)*(h g)
+ v0' = (center c) + (delta, delta, delta)
+ v1' = (center c) + (delta, delta, -delta)
+ v2' = (center c) + (delta, -delta, -delta)
+ v3' = (center c) + (delta, -delta, delta)
+
+
+-- The bottom face of the cube.
+face2 :: Cube -> Face
+face2 c = Face c v0' v1' v2' v3'
+ where
+ g = grid c
+ delta = (1/2)*(h g)
+ v0' = (center c) + (delta, delta, -delta)
+ v1' = (center c) + (-delta, delta, -delta)
+ v2' = (center c) + (-delta, -delta, -delta)
+ v3' = (center c) + (delta, -delta, -delta)
+
+
+-- The left face of the cube.
+face3 :: Cube -> Face
+face3 c = Face c v0' v1' v2' v3'
+ where
+ g = grid c
+ delta = (1/2)*(h g)
+ v0' = (center c) + (-delta, delta, -delta)
+ v1' = (center c) + (-delta, delta, delta)
+ v2' = (center c) + (-delta, -delta, delta)
+ v3' = (center c) + (-delta, -delta, -delta)
+
+
+-- The front face of the cube.
+face4 :: Cube -> Face
+face4 c = Face c v0' v1' v2' v3'
+ where
+ g = grid c
+ delta = (1/2)*(h g)
+ v0' = (center c) + (-delta, -delta, delta)
+ v1' = (center c) + (delta, -delta, delta)
+ v2' = (center c) + (delta, -delta, -delta)
+ v3' = (center c) + (-delta, -delta, -delta)
+
+
+-- The back face of the cube.
+face5 :: Cube -> Face
+face5 c = Face c v0' v1' v2' v3'
+ where
+ g = grid c
+ delta = (1/2)*(h g)
+ v0' = (center c) + (-delta, delta, -delta)
+ v1' = (center c) + (delta, delta, -delta)
+ v2' = (center c) + (delta, delta, delta)
+ v3' = (center c) + (-delta, delta, delta)
+
+
+tetrahedron0 :: Face -> Tetrahedron
+tetrahedron0 f =
+ Tetrahedron c v0' v1' v2' v3'
+ where
+ c = cube f
+ v0' = v0 f
+ v1' = v1 f
+ v2' = center f
+ v3' = center c
+
+tetrahedron1 :: Face -> Tetrahedron
+tetrahedron1 f =
+ Tetrahedron c v0' v1' v2' v3'
+ where
+ c = cube f
+ v0' = v1 f
+ v1' = v2 f
+ v2' = center f
+ v3' = center c
+
+
+tetrahedron2 :: Face -> Tetrahedron
+tetrahedron2 f =
+ Tetrahedron c v0' v1' v2' v3'
+ where
+ c = cube f
+ v0' = v2 f
+ v1' = v3 f
+ v2' = center f
+ v3' = center c
+
+
+tetrahedron3 :: Face -> Tetrahedron
+tetrahedron3 f =
+ Tetrahedron c v0' v1' v2' v3'
+ where
+ c = cube f
+ v0' = v3 f
+ v1' = v0 f
+ v2' = center f
+ v3' = center c
+
+tetrahedrons :: Cube -> [Tetrahedron]
+tetrahedrons c =
+ concat [
+ [tetrahedron0 f0, tetrahedron1 f0, tetrahedron2 f0, tetrahedron3 f0],
+ [tetrahedron0 f1, tetrahedron1 f1, tetrahedron2 f1, tetrahedron3 f2],
+ [tetrahedron0 f2, tetrahedron1 f2, tetrahedron2 f2, tetrahedron3 f2],
+ [tetrahedron0 f3, tetrahedron1 f3, tetrahedron2 f3, tetrahedron3 f3],
+ [tetrahedron0 f4, tetrahedron1 f4, tetrahedron2 f4, tetrahedron3 f4],
+ [tetrahedron0 f5, tetrahedron1 f5, tetrahedron2 f5, tetrahedron3 f5] ]
+ where
+ f0 = face0 c
+ f1 = face1 c
+ f2 = face2 c
+ f3 = face3 c
+ f4 = face4 c
+ f5 = face5 c
--- /dev/null
+-- | 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.
+module Grid (empty_grid,
+ function_values,
+ Grid,
+ h,
+ make_grid)
+where
+
+-- | 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 :: [[[Double]]] }
+ deriving (Eq, Show)
+
+
+-- | The constructor that we want people to use. If we're passed a
+-- non-positive grid size, we throw an error.
+make_grid :: Double -> [[[Double]]] -> Grid
+make_grid grid_size values
+ | grid_size <= 0 = error "grid size must be positive"
+ | otherwise = Grid grid_size values
+
+
+-- | Creates an empty grid with grid size 1.
+empty_grid :: Grid
+empty_grid = Grid 1 [[[]]]
--- /dev/null
+module Main
+where
+
+import Cube
+import Face
+import Grid
+import Misc (flatten)
+import Point
+import RealFunction
+import Tetrahedron
+import ThreeDimensional
+
+trilinear :: [[[Double]]]
+trilinear = [ [ [ 1, 2, 3 ],
+ [ 1, 3, 5 ],
+ [ 1, 4, 7 ] ],
+ [ [ 1, 2, 3 ],
+ [ 1, 4, 7 ],
+ [ 1, 6, 11 ] ],
+ [ [ 1, 2, 3 ],
+ [ 1, 5, 9 ],
+ [ 1, 8, 15 ]]]
+
+zeros :: [[[Double]]]
+zeros = [ [ [ 0, 0, 0 ],
+ [ 0, 0, 0 ],
+ [ 0, 0, 0 ] ],
+ --
+ [ [ 0, 0, 0 ],
+ [ 0, 0, 0 ],
+ [ 0, 0, 0 ] ],
+ --
+ [ [ 0, 0, 0 ],
+ [ 0, 0, 0 ],
+ [ 0, 0, 0 ]]]
+
+dummy :: [[[Double]]]
+dummy = [ [ [ 0, 1, 2 ],
+ [ 3, 4, 5 ],
+ [ 6, 7, 8 ] ],
+ --
+ [ [ 9, 10, 11 ],
+ [ 12, 13, 14 ],
+ [ 15, 16, 17 ] ],
+ --
+ [ [ 18, 19, 20 ],
+ [ 21, 22, 23 ],
+ [ 24, 25, 26 ]]]
+
+
+find_point_value :: RealFunction Point
+find_point_value p = poly p
+ where
+ g0 = make_grid 1 trilinear
+ the_cubes = flatten (cubes g0)
+ good_cubes = filter ((flip contains_point) p) the_cubes
+ target_cube = good_cubes !! 0
+ good_tets = filter ((flip contains_point) p) (tetrahedrons target_cube)
+ target_tetrahedron = good_tets !! 0
+ poly = polynomial target_tetrahedron
+
+main :: IO ()
+main = do
+ putStrLn $ show $ find_point_value (0,0,0)
+ putStrLn $ show $ find_point_value (1,0,0)
+ putStrLn $ show $ find_point_value (2,0,0)
+ putStrLn $ show $ find_point_value (0,1,0)
+ putStrLn $ show $ find_point_value (1,1,0)
+ putStrLn $ show $ find_point_value (2,1,0)
+ putStrLn $ show $ find_point_value (0,2,0)
+ putStrLn $ show $ find_point_value (1,2,0)
+ putStrLn $ show $ find_point_value (2,2,0)
+ putStrLn $ show $ find_point_value (0,0,1)
+ putStrLn $ show $ find_point_value (1,0,1)
+ putStrLn $ show $ find_point_value (2,0,1)
+ putStrLn $ show $ find_point_value (0,1,1)
+ putStrLn $ show $ find_point_value (1,1,1)
+ putStrLn $ show $ find_point_value (2,1,1)
+ putStrLn $ show $ find_point_value (0,2,1)
+ putStrLn $ show $ find_point_value (1,2,1)
+ putStrLn $ show $ find_point_value (2,2,1)
+ putStrLn $ show $ find_point_value (0,0,2)
+ putStrLn $ show $ find_point_value (1,0,2)
+ putStrLn $ show $ find_point_value (2,0,2)
+ putStrLn $ show $ find_point_value (0,1,2)
+ putStrLn $ show $ find_point_value (1,1,2)
+ putStrLn $ show $ find_point_value (2,1,2)
+ putStrLn $ show $ find_point_value (0,2,2)
+ putStrLn $ show $ find_point_value (1,2,2)
+ putStrLn $ show $ find_point_value (2,2,2)
+ -- let g0 = make_grid 1 trilinear
+ -- let the_cubes = flatten (cubes g0)
+ -- putStrLn $ show $ the_cubes
+ -- let p = (2, 0, 0)
+ -- let target_cubes = filter ((flip contains_point) p) the_cubes
+ -- putStrLn $ show $ target_cubes
+ -- let target_cube = (take 1 target_cubes) !! 0
+ -- putStrLn $ show $ target_cube
+ -- let target_tetrahedra = filter ((flip contains_point) p) (tetrahedrons target_cube)
+ -- let target_tetrahedron = (take 1 target_tetrahedra) !! 0
+ -- putStrLn $ show $ target_tetrahedron
+ -- let poly = polynomial target_tetrahedron
+ -- putStrLn $ show $ poly
+ -- putStrLn $ show $ poly p
--- /dev/null
+-- | The Misc module contains helper functions that seem out of place
+-- anywhere else.
+module Misc
+where
+
+
+-- | The standard factorial function. See
+-- <http://www.willamette.edu/~fruehr/haskell/evolution.html> for
+-- possible improvements.
+factorial :: Int -> Int
+factorial n
+ | n <= 1 = 1
+ | n > 20 = error "integer overflow in factorial function"
+ | otherwise = product [1..n]
+
+
+-- | Takes a three-dimensional list, and flattens it into a
+-- one-dimensional one.
+flatten :: [[[a]]] -> [a]
+flatten xs = concat $ concat xs
--- /dev/null
+{-# LANGUAGE TypeSynonymInstances #-}
+
+module Point
+where
+
+type Point = (Double, Double, Double)
+
+x_coord :: Point -> Double
+x_coord (x, _, _) = x
+
+y_coord :: Point -> Double
+y_coord (_, y, _) = y
+
+z_coord :: Point -> Double
+z_coord (_, _, z) = z
+
+instance Num Point where
+ p1 + p2 = (x1+x2, y1+y2, z1+z2)
+ where
+ x1 = x_coord p1
+ x2 = x_coord p2
+ y1 = y_coord p1
+ y2 = y_coord p2
+ z1 = z_coord p1
+ z2 = z_coord p2
+
+ p1 - p2 = (x1-x2, y1-y2, z1-z2)
+ where
+ x1 = x_coord p1
+ x2 = x_coord p2
+ y1 = y_coord p1
+ y2 = y_coord p2
+ z1 = z_coord p1
+ z2 = z_coord p2
+
+ p1 * p2 = (x1*x2, y1*y2, z1*z2)
+ where
+ x1 = x_coord p1
+ x2 = x_coord p2
+ y1 = y_coord p1
+ y2 = y_coord p2
+ z1 = z_coord p1
+ z2 = z_coord p2
+
+ abs (x, y, z) = (abs x, abs y, abs z)
+ signum (x, y, z) = (signum x, signum y, signum z)
+ fromInteger n = (fromInteger n, fromInteger n, fromInteger n)
+
+
+scale :: Point -> Double -> Point
+scale (x, y, z) d = (x*d, y*d, z*d)
--- /dev/null
+{-# LANGUAGE TypeSynonymInstances #-}
+
+module RealFunction
+where
+
+
+type RealFunction a = (a -> Double)
+
+instance Show (RealFunction a) where
+ show _ = "Real Function"
+
+instance Eq (RealFunction a) where
+ _ == _ = False
+
+instance Num (RealFunction a) where
+ f1 + f2 = \x -> (f1 x) + (f2 x)
+ f1 - f2 = \x -> (f1 x) - (f2 x)
+ f1 * f2 = \x -> (f1 x) * (f2 x)
+ negate f = \x -> -1 * (f x)
+ abs f = \x -> abs (f x)
+ signum f = \x -> signum (f x)
+ fromInteger i = \_ -> (fromInteger i)
+
+
+-- Takes a constant, and a function as arguments. Returns a new
+-- function representing the original function times the constant.
+cmult :: Double -> (RealFunction a) -> (RealFunction a)
+cmult coeff f = (*coeff) . f
+
+-- Takes a function f and an exponent n. Returns a new function, f^n.
+fexp :: (RealFunction a) -> Int -> (RealFunction a)
+fexp f n
+ | n == 0 = (\_ -> 1)
+ | otherwise = \x -> (f x)^n
--- /dev/null
+module Tests.Cube
+where
+
+import Test.QuickCheck
+
+import Cube
+import Grid (Grid)
+import Tests.Grid ()
+
+instance Arbitrary Cube where
+ arbitrary = do
+ g' <- arbitrary :: Gen Grid
+ i' <- choose (coordmin, coordmax)
+ j' <- choose (coordmin, coordmax)
+ k' <- choose (coordmin, coordmax)
+ d' <- arbitrary :: Gen Double
+ return (Cube g' i' j' k' d')
+ where
+ coordmin = -268435456 -- -(2^29 / 2)
+ coordmax = 268435456 -- +(2^29 / 2)
--- /dev/null
+module Tests.Face
+where
+
+import Test.QuickCheck
+
+import Cube (Cube(grid))
+import Face (tetrahedrons)
+import Grid (Grid(h))
+import Tetrahedron (volume)
+
+-- QuickCheck Tests.
+prop_all_volumes_nonnegative :: Cube -> Property
+prop_all_volumes_nonnegative c =
+ (delta > 0) ==> (null negative_volumes)
+ where
+ delta = h (grid c)
+ ts = tetrahedrons c
+ volumes = map volume ts
+ negative_volumes = filter (< 0) volumes
--- /dev/null
+module Tests.Grid
+where
+
+import Test.QuickCheck
+import Grid
+
+instance Arbitrary Grid where
+ arbitrary = do
+ (Positive h') <- arbitrary :: Gen (Positive Double)
+ fv <- arbitrary :: Gen [[[Double]]]
+ return (make_grid h' fv)
--- /dev/null
+module Tests.Misc
+where
+
+import Test.HUnit
+import Test.QuickCheck
+
+import Misc
+
+prop_factorial_greater :: Int -> Property
+prop_factorial_greater n =
+ n <= 20 ==> factorial n >= n
+
+
+test_flatten1 :: Test
+test_flatten1 =
+ TestCase $ assertEqual "flatten actually works" expected_list actual_list
+ where
+ target = [[[1::Int]], [[2, 3]]]
+ expected_list = [1, 2, 3]
+ actual_list = flatten target
+
+misc_tests :: [Test]
+misc_tests = [ test_flatten1 ]
--- /dev/null
+module Tests.Tetrahedron
+where
+
+import Test.HUnit
+import Test.QuickCheck
+
+import Cube
+import Point
+import Tests.Cube()
+import Tetrahedron
+import ThreeDimensional
+
+instance Arbitrary Tetrahedron where
+ arbitrary = do
+ rnd_c0 <- arbitrary :: Gen Cube
+ rnd_v0 <- arbitrary :: Gen Point
+ rnd_v1 <- arbitrary :: Gen Point
+ rnd_v2 <- arbitrary :: Gen Point
+ rnd_v3 <- arbitrary :: Gen Point
+ return (Tetrahedron rnd_c0 rnd_v0 rnd_v1 rnd_v2 rnd_v3)
+
+almost_equals :: Double -> Double -> Bool
+almost_equals x y = (abs (x - y)) < 0.0001
+
+(~=) :: Double -> Double -> Bool
+(~=) = almost_equals
+
+
+-- HUnit Tests
+
+-- Since p0, p1, p2 are in clockwise order, we expect the volume here
+-- to be negative.
+test_volume1 :: Test
+test_volume1 =
+ TestCase $ assertEqual "volume is correct" True (vol ~= (-1/3))
+ where
+ p0 = (0, -0.5, 0)
+ p1 = (0, 0.5, 0)
+ p2 = (2, 0, 0)
+ p3 = (1, 0, 1)
+ t = Tetrahedron { cube = empty_cube,
+ v0 = p0,
+ v1 = p1,
+ v2 = p2,
+ v3 = p3 }
+ vol = volume t
+
+
+-- Now, p0, p1, and p2 are in counter-clockwise order. The volume
+-- should therefore be positive.
+test_volume2 :: Test
+test_volume2 =
+ TestCase $ assertEqual "volume is correct" True (vol ~= (1/3))
+ where
+ p0 = (0, -0.5, 0)
+ p1 = (2, 0, 0)
+ p2 = (0, 0.5, 0)
+ p3 = (1, 0, 1)
+ t = Tetrahedron { cube = empty_cube,
+ v0 = p0,
+ v1 = p1,
+ v2 = p2,
+ v3 = p3 }
+ vol = volume t
+
+test_contains_point1 :: Test
+test_contains_point1 =
+ TestCase $ assertEqual "contains an inner point" True (contains_point t inner_point)
+ where
+ p0 = (0, -0.5, 0)
+ p1 = (0, 0.5, 0)
+ p2 = (2, 0, 0)
+ p3 = (1, 0, 1)
+ inner_point = (1, 0, 0.5)
+ t = Tetrahedron { cube = empty_cube,
+ v0 = p0,
+ v1 = p1,
+ v2 = p2,
+ v3 = p3 }
+
+
+test_doesnt_contain_point1 :: Test
+test_doesnt_contain_point1 =
+ TestCase $ assertEqual "doesn't contain an exterior point" False (contains_point t exterior_point)
+ where
+ p0 = (0, -0.5, 0)
+ p1 = (0, 0.5, 0)
+ p2 = (2, 0, 0)
+ p3 = (1, 0, 1)
+ exterior_point = (5, 2, -9.0212)
+ c_empty = empty_cube
+ t = Tetrahedron { cube = c_empty,
+ v0 = p0,
+ v1 = p1,
+ v2 = p2,
+ v3 = p3 }
+
+
+test_doesnt_contain_point2 :: Test
+test_doesnt_contain_point2 =
+ TestCase $ assertEqual "doesn't contain an exterior point" False (contains_point t exterior_point)
+ where
+ p0 = (0, 1, 1)
+ p1 = (1, 1, 1)
+ p2 = (0.5, 0.5, 1)
+ p3 = (0.5, 0.5, 0.5)
+ exterior_point = (0, 0, 0)
+ c_empty = empty_cube
+ t = Tetrahedron { cube = c_empty,
+ v0 = p0,
+ v1 = p1,
+ v2 = p2,
+ v3 = p3 }
+
+test_doesnt_contain_point3 :: Test
+test_doesnt_contain_point3 =
+ TestCase $ assertEqual "doesn't contain an exterior point" False (contains_point t exterior_point)
+ where
+ p0 = (1, 1, 1)
+ p1 = (1, 0, 1)
+ p2 = (0.5, 0.5, 1)
+ p3 = (0.5, 0.5, 0.5)
+ exterior_point = (0, 0, 0)
+ c_empty = empty_cube
+ t = Tetrahedron { cube = c_empty,
+ v0 = p0,
+ v1 = p1,
+ v2 = p2,
+ v3 = p3 }
+
+test_doesnt_contain_point4 :: Test
+test_doesnt_contain_point4 =
+ TestCase $ assertEqual "doesn't contain an exterior point" False (contains_point t exterior_point)
+ where
+ p0 = (1, 0, 1)
+ p1 = (0, 0, 1)
+ p2 = (0.5, 0.5, 1)
+ p3 = (0.5, 0.5, 0.5)
+ exterior_point = (0, 0, 0)
+ c_empty = empty_cube
+ t = Tetrahedron { cube = c_empty,
+ v0 = p0,
+ v1 = p1,
+ v2 = p2,
+ v3 = p3 }
+
+test_doesnt_contain_point5 :: Test
+test_doesnt_contain_point5 =
+ TestCase $ assertEqual "doesn't contain an exterior point" False (contains_point t exterior_point)
+ where
+ p0 = (0, 0, 1)
+ p1 = (0, 1, 1)
+ p2 = (0.5, 0.5, 1)
+ p3 = (0.5, 0.5, 0.5)
+ exterior_point = (0, 0, 0)
+ c_empty = empty_cube
+ t = Tetrahedron { cube = c_empty,
+ v0 = p0,
+ v1 = p1,
+ v2 = p2,
+ v3 = p3 }
+
+tetrahedron_tests :: [Test]
+tetrahedron_tests = [test_volume1,
+ test_volume2,
+ test_contains_point1,
+ test_doesnt_contain_point1,
+ test_doesnt_contain_point2,
+ test_doesnt_contain_point3,
+ test_doesnt_contain_point4,
+ test_doesnt_contain_point5 ]
+
+prop_b0_v0_always_unity :: Tetrahedron -> Property
+prop_b0_v0_always_unity t =
+ (volume t) > 0 ==> (b0 t) (v0 t) ~= 1.0
+
+prop_b0_v1_always_zero :: Tetrahedron -> Property
+prop_b0_v1_always_zero t =
+ (volume t) > 0 ==> (b0 t) (v1 t) ~= 0
+
+prop_b0_v2_always_zero :: Tetrahedron -> Property
+prop_b0_v2_always_zero t =
+ (volume t) > 0 ==> (b0 t) (v2 t) ~= 0
+
+prop_b0_v3_always_zero :: Tetrahedron -> Property
+prop_b0_v3_always_zero t =
+ (volume t) > 0 ==> (b0 t) (v3 t) ~= 0
+
+prop_b1_v1_always_unity :: Tetrahedron -> Property
+prop_b1_v1_always_unity t =
+ (volume t) > 0 ==> (b1 t) (v1 t) ~= 1.0
+
+prop_b1_v0_always_zero :: Tetrahedron -> Property
+prop_b1_v0_always_zero t =
+ (volume t) > 0 ==> (b1 t) (v0 t) ~= 0
+
+prop_b1_v2_always_zero :: Tetrahedron -> Property
+prop_b1_v2_always_zero t =
+ (volume t) > 0 ==> (b1 t) (v2 t) ~= 0
+
+prop_b1_v3_always_zero :: Tetrahedron -> Property
+prop_b1_v3_always_zero t =
+ (volume t) > 0 ==> (b1 t) (v3 t) ~= 0
+
+prop_b2_v2_always_unity :: Tetrahedron -> Property
+prop_b2_v2_always_unity t =
+ (volume t) > 0 ==> (b2 t) (v2 t) ~= 1.0
+
+prop_b2_v0_always_zero :: Tetrahedron -> Property
+prop_b2_v0_always_zero t =
+ (volume t) > 0 ==> (b2 t) (v0 t) ~= 0
+
+prop_b2_v1_always_zero :: Tetrahedron -> Property
+prop_b2_v1_always_zero t =
+ (volume t) > 0 ==> (b2 t) (v1 t) ~= 0
+
+prop_b2_v3_always_zero :: Tetrahedron -> Property
+prop_b2_v3_always_zero t =
+ (volume t) > 0 ==> (b2 t) (v3 t) ~= 0
+
+prop_b3_v3_always_unity :: Tetrahedron -> Property
+prop_b3_v3_always_unity t =
+ (volume t) > 0 ==> (b3 t) (v3 t) ~= 1.0
+
+prop_b3_v0_always_zero :: Tetrahedron -> Property
+prop_b3_v0_always_zero t =
+ (volume t) > 0 ==> (b3 t) (v0 t) ~= 0
+
+prop_b3_v1_always_zero :: Tetrahedron -> Property
+prop_b3_v1_always_zero t =
+ (volume t) > 0 ==> (b3 t) (v1 t) ~= 0
+
+prop_b3_v2_always_zero :: Tetrahedron -> Property
+prop_b3_v2_always_zero t =
+ (volume t) > 0 ==> (b3 t) (v2 t) ~= 0
+
+
+-- Used for convenience in the next few tests.
+p :: Tetrahedron -> Int -> Int -> Int -> Int -> Double
+p t i j k l = (polynomial t) (xi t i j k l)
+
+-- | Given in Sorokina and Zeilfelder, p. 78.
+prop_c3000_identity :: Tetrahedron -> Property
+prop_c3000_identity t =
+ (volume t) > 0 ==>
+ c t 3 0 0 0 ~= p t 3 0 0 0
+
+-- | Given in Sorokina and Zeilfelder, p. 78.
+prop_c2100_identity :: Tetrahedron -> Property
+prop_c2100_identity t =
+ (volume t) > 0 ==>
+ c t 2 1 0 0 ~= (term1 - term2 + term3 - term4)
+ where
+ term1 = (1/3)*(p t 0 3 0 0)
+ term2 = (5/6)*(p t 3 0 0 0)
+ term3 = 3*(p t 2 1 0 0)
+ term4 = (3/2)*(p t 1 2 0 0)
+
+-- | Given in Sorokina and Zeilfelder, p. 78.
+prop_c1110_identity :: Tetrahedron -> Property
+prop_c1110_identity t =
+ (volume t) > 0 ==>
+ c t 1 1 1 0 ~= (term1 + term2 - term3 - term4)
+ where
+ term1 = (1/3)*((p t 3 0 0 0) + (p t 0 3 0 0) + (p t 0 0 3 0))
+ term2 = (9/2)*(p t 1 1 1 0)
+ term3 = (3/4)*((p t 2 1 0 0) + (p t 1 2 0 0) + (p t 2 0 1 0))
+ term4 = (3/4)*((p t 1 0 2 0) + (p t 0 2 1 0) + (p t 0 1 2 0))
--- /dev/null
+module Tetrahedron
+where
+
+import Numeric.LinearAlgebra hiding (i, scale)
+
+import Cube (back,
+ Cube(datum),
+ down,
+ front,
+ Gridded,
+ left,
+ right,
+ top)
+import Misc (factorial)
+import Point
+import RealFunction
+import ThreeDimensional
+
+data Tetrahedron = Tetrahedron { cube :: Cube,
+ v0 :: Point,
+ v1 :: Point,
+ v2 :: Point,
+ v3 :: Point }
+ deriving (Eq)
+
+instance Show Tetrahedron where
+ show t = "Tetrahedron (Cube: " ++ (show (cube t)) ++ ") " ++
+ "(v0: " ++ (show (v0 t)) ++ ") (v1: " ++ (show (v1 t)) ++
+ ") (v2: " ++ (show (v2 t)) ++ ") (v3: " ++ (show (v3 t)) ++
+ ")\n\n"
+
+instance Gridded Tetrahedron where
+ back t = back (cube t)
+ down t = down (cube t)
+ front t = front (cube t)
+ left t = left (cube t)
+ right t = right (cube t)
+ top t = top (cube t)
+
+
+instance ThreeDimensional Tetrahedron where
+ center t = ((v0 t) + (v1 t) + (v2 t) + (v3 t)) `scale` (1/4)
+ contains_point t p =
+ (b0 t p) >= 0 && (b1 t p) >= 0 && (b2 t p) >= 0 && (b3 t p) >= 0
+
+
+polynomial :: Tetrahedron -> (RealFunction Point)
+polynomial t =
+ sum [ (c t i j k l) `cmult` (beta t i j k l) | i <- [0..3],
+ j <- [0..3],
+ k <- [0..3],
+ l <- [0..3],
+ i + j + k + l == 3]
+
+
+-- | Returns the domain point of t with indices i,j,k,l.
+xi :: Tetrahedron -> Int -> Int -> Int -> Int -> Point
+xi t i j k l
+ | i + j + k + l == 3 = weighted_sum `scale` (1/3)
+ | otherwise = error "xi index out of bounds"
+ where
+ v0' = (v0 t) `scale` (fromIntegral i)
+ v1' = (v1 t) `scale` (fromIntegral j)
+ v2' = (v2 t) `scale` (fromIntegral k)
+ v3' = (v3 t) `scale` (fromIntegral l)
+ weighted_sum = v0' + v1' + v2' + v3'
+
+
+
+-- | The Bernstein polynomial on t with indices i,j,k,l. Denoted by a
+-- capital 'B' in the Sorokina/Zeilfelder paper.
+beta :: Tetrahedron -> Int -> Int -> Int -> Int -> (RealFunction Point)
+beta t i j k l
+ | (i + j + k + l == 3) =
+ coefficient `cmult` (b0_term * b1_term * b2_term * b3_term)
+ | otherwise = error "basis function index out of bounds"
+ where
+ denominator = (factorial i)*(factorial j)*(factorial k)*(factorial l)
+ coefficient = 6 / (fromIntegral denominator)
+ b0_term = (b0 t) `fexp` i
+ b1_term = (b1 t) `fexp` j
+ b2_term = (b2 t) `fexp` k
+ b3_term = (b3 t) `fexp` l
+
+
+c :: Tetrahedron -> Int -> Int -> Int -> Int -> Double
+c x 0 0 3 0 = datum $ (1/8) * (i + f + l + t + lt + fl + ft + flt)
+ where
+ f = front x
+ flt = front (left (top x))
+ fl = front (left x)
+ ft = front (top x)
+ i = cube x
+ l = left x
+ lt = left (top x)
+ t = top x
+
+
+c x 0 0 0 3 = datum $ (1/8) * (i + f + r + t + rt + fr + ft + frt)
+ where
+ f = front x
+ fr = front (right x)
+ frt = front (right (top x))
+ ft = front (top x)
+ i = cube x
+ r = right x
+ rt = right (top x)
+ t = top x
+
+
+c x 0 0 2 1 = datum $ (5/24)*(i + f + t + ft) + (1/24)*(l + fl + lt + flt)
+ where
+ f = front x
+ flt = front (left (top x))
+ fl = front (left x)
+ ft = front (top x)
+ i = cube x
+ l = left x
+ lt = left (top x)
+ t = top x
+
+
+c x 0 0 1 2 = datum $ (5/24)*(i + f + t + ft) + (1/24)*(r + fr + rt + frt)
+ where
+ f = front x
+ frt = front (right (top x))
+ fr = front (right x)
+ ft = front (top x)
+ i = cube x
+ r = right x
+ rt = right (top x)
+ t = top x
+
+
+c x 0 1 2 0 = datum $
+ (5/24)*(i + f) +
+ (1/8)*(l + t + fl + ft) +
+ (1/24)*(lt + flt)
+ where
+ f = front x
+ flt = front (left (top x))
+ fl = front (left x)
+ ft = front (top x)
+ i = cube x
+ l = left x
+ lt = left (top x)
+ t = top x
+
+
+c x 0 1 0 2 = datum $
+ (5/24)*(i + f) +
+ (1/8)*(r + t + fr + ft) +
+ (1/24)*(rt + frt)
+ where
+ f = front x
+ fr = front (right x)
+ frt = front (right (top x))
+ ft = front (top x)
+ i = cube x
+ r = right x
+ rt = right (top x)
+ t = top x
+
+
+c x 0 1 1 1 = datum $
+ (13/48)*(i + f) +
+ (7/48)*(t + ft) +
+ (1/32)*(l + r + fl + fr) +
+ (1/96)*(lt + rt + flt + frt)
+ where
+ f = front x
+ flt = front (left (top x))
+ fl = front (left x)
+ fr = front (right x)
+ frt = front (right (top x))
+ ft = front (top x)
+ i = cube x
+ l = left x
+ lt = left (top x)
+ r = right x
+ rt = right (top x)
+ t = top x
+
+
+c x 0 2 1 0 = datum $
+ (13/48)*(i + f) +
+ (17/192)*(l + t + fl + ft) +
+ (1/96)*(lt + flt) +
+ (1/64)*(r + d + fr + fd) +
+ (1/192)*(rt + ld + frt + fld)
+ where
+ d = down x
+ f = front x
+ fd = front (down x)
+ fld = front (left (down x))
+ flt = front (left (top x))
+ fl = front (left x)
+ fr = front (right x)
+ frt = front (right (top x))
+ ft = front (top x)
+ i = cube x
+ l = left x
+ ld = left (down x)
+ lt = left (top x)
+ r = right x
+ rt = right (top x)
+ t = top x
+
+
+c x 0 2 0 1 = datum $
+ (13/48)*(i + f) +
+ (17/192)*(r + t + fr + ft) +
+ (1/96)*(rt + frt) +
+ (1/64)*(l + d + fl + fd) +
+ (1/192)*(rd + lt + flt + frd)
+ where
+ d = down x
+ f = front x
+ fd = front (down x)
+ flt = front (left (top x))
+ fl = front (left x)
+ frd = front (right (down x))
+ fr = front (right x)
+ frt = front (right (top x))
+ ft = front (top x)
+ i = cube x
+ l = left x
+ lt = left (top x)
+ r = right x
+ rd = right (down x)
+ rt = right (top x)
+ t = top x
+
+
+c x 0 3 0 0 = datum $
+ (13/48)*(i + f) +
+ (5/96)*(l + r + t + d + fl + fr + ft + fd) +
+ (1/192)*(rt + rd + lt + ld + frt + frd + flt + fld)
+ where
+ d = down x
+ f = front x
+ fd = front (down x)
+ fld = front (left (down x))
+ flt = front (left (top x))
+ fl = front (left x)
+ frd = front (right (down x))
+ fr = front (right x)
+ frt = front (right (top x))
+ ft = front (top x)
+ i = cube x
+ l = left x
+ ld = left (down x)
+ lt = left (top x)
+ r = right x
+ rd = right (down x)
+ rt = right (top x)
+ t = top x
+
+c x 1 0 2 0 = datum $ (1/4)*i + (1/6)*(f + l + t) + (1/12)*(lt + fl + ft)
+ where
+ f = front x
+ fl = front (left x)
+ ft = front (top x)
+ i = cube x
+ l = left x
+ lt = left (top x)
+ t = top x
+
+
+c x 1 0 0 2 = datum $ (1/4)*i + (1/6)*(f + r + t) + (1/12)*(rt + fr + ft)
+ where
+ f = front x
+ fr = front (right x)
+ ft = front (top x)
+ i = cube x
+ r = right x
+ rt = right (top x)
+ t = top x
+
+
+c x 1 0 1 1 = datum $
+ (1/3)*i +
+ (5/24)*(f + t) +
+ (1/12)*ft +
+ (1/24)*(l + r) +
+ (1/48)*(lt + rt + fl + fr)
+ where
+ f = front x
+ fl = front (left x)
+ fr = front (right x)
+ ft = front (top x)
+ i = cube x
+ l = left x
+ lt = left (top x)
+ r = right x
+ rt = right (top x)
+ t = top x
+
+
+c x 1 1 1 0 = datum $
+ (1/3)*i +
+ (5/24)*f +
+ (1/8)*(l + t) +
+ (5/96)*(fl + ft) +
+ (1/48)*(d + r + lt) +
+ (1/96)*(fd + ld + rt + fr)
+ where
+ d = down x
+ f = front x
+ fd = front (down x)
+ fl = front (left x)
+ fr = front (right x)
+ ft = front (top x)
+ i = cube x
+ l = left x
+ ld = left (down x)
+ lt = left (top x)
+ r = right x
+ rt = right (top x)
+ t = top x
+
+
+c x 1 1 0 1 = datum $
+ (1/3)*i +
+ (5/24)*f +
+ (1/8)*(r + t) +
+ (5/96)*(fr + ft) +
+ (1/48)*(d + l + rt) +
+ (1/96)*(fd + lt + rd + fl)
+ where
+ d = down x
+ f = front x
+ fd = front (down x)
+ fl = front (left x)
+ fr = front (right x)
+ ft = front (top x)
+ i = cube x
+ l = left x
+ lt = left (top x)
+ r = right x
+ rd = right (down x)
+ rt = right (top x)
+ t = top x
+
+
+
+c x 1 2 0 0 = datum $
+ (1/3)*i +
+ (5/24)*f +
+ (7/96)*(l + r + t + d) +
+ (1/32)*(fl + fr + ft + fd) +
+ (1/96)*(rt + rd + lt + ld)
+ where
+ d = down x
+ f = front x
+ fd = front (down x)
+ fl = front (left x)
+ fr = front (right x)
+ ft = front (top x)
+ i = cube x
+ l = left x
+ ld = left (down x)
+ lt = left (top x)
+ r = right x
+ rd = right (down x)
+ rt = right (top x)
+ t = top x
+
+
+c x 2 0 1 0 = datum $
+ (3/8)*i +
+ (7/48)*(f + t + l) +
+ (1/48)*(r + d + b + lt + fl + ft) +
+ (1/96)*(rt + bt + fr + fd + ld + bl)
+ where
+ b = back x
+ bl = back (left x)
+ bt = back (top x)
+ d = down x
+ f = front x
+ fd = front (down x)
+ fl = front (left x)
+ fr = front (right x)
+ ft = front (top x)
+ i = cube x
+ l = left x
+ ld = left (down x)
+ lt = left (top x)
+ r = right x
+ rt = right (top x)
+ t = top x
+
+
+c x 2 0 0 1 = datum $
+ (3/8)*i +
+ (7/48)*(f + t + r) +
+ (1/48)*(l + d + b + rt + fr + ft) +
+ (1/96)*(lt + bt + fl + fd + rd + br)
+ where
+ b = back x
+ br = back (right x)
+ bt = back (top x)
+ d = down x
+ f = front x
+ fd = front (down x)
+ fl = front (left x)
+ fr = front (right x)
+ ft = front (top x)
+ i = cube x
+ l = left x
+ lt = left (top x)
+ r = right x
+ rd = right (down x)
+ rt = right (top x)
+ t = top x
+
+
+
+c x 2 1 0 0 = datum $
+ (3/8)*i +
+ (1/12)*(t + r + l + d) +
+ (1/64)*(ft + fr + fl + fd) +
+ (7/48)*f +
+ (1/48)*b +
+ (1/96)*(rt + ld + lt + rd) +
+ (1/192)*(bt + br + bl + bd)
+ where
+ b = back x
+ bd = back (down x)
+ bl = back (left x)
+ br = back (right x)
+ bt = back (top x)
+ d = down x
+ f = front x
+ fd = front (down x)
+ fl = front (left x)
+ fr = front (right x)
+ ft = front (top x)
+ i = cube x
+ l = left x
+ ld = left (down x)
+ lt = left (top x)
+ r = right x
+ rd = right (down x)
+ rt = right (top x)
+ t = top x
+
+
+c x 3 0 0 0 = datum $
+ (3/8)*i +
+ (1/12)*(t + f + l + r + d + b) +
+ (1/96)*(lt + fl + ft + rt + bt + fr) +
+ (1/96)*(fd + ld + bd + br + rd + bl)
+ where
+ b = back x
+ bd = back (down x)
+ bl = back (left x)
+ br = back (right x)
+ bt = back (top x)
+ d = down x
+ f = front x
+ fd = front (down x)
+ fl = front (left x)
+ fr = front (right x)
+ ft = front (top x)
+ i = cube x
+ l = left x
+ ld = left (down x)
+ lt = left (top x)
+ r = right x
+ rd = right (down x)
+ rt = right (top x)
+ t = top x
+
+
+c _ _ _ _ _ = error "coefficient index out of bounds"
+
+
+
+vol_matrix :: Tetrahedron -> Matrix Double
+vol_matrix t = (4><4) $
+ [1, 1, 1, 1,
+ x1, x2, x3, x4,
+ y1, y2, y3, y4,
+ z1, z2, z3, z4 ]
+ where
+ x1 = x_coord (v0 t)
+ x2 = x_coord (v1 t)
+ x3 = x_coord (v2 t)
+ x4 = x_coord (v3 t)
+ y1 = y_coord (v0 t)
+ y2 = y_coord (v1 t)
+ y3 = y_coord (v2 t)
+ y4 = y_coord (v3 t)
+ z1 = z_coord (v0 t)
+ z2 = z_coord (v1 t)
+ z3 = z_coord (v2 t)
+ z4 = z_coord (v3 t)
+
+-- Computed using the formula from Lai & Schumaker, Definition 15.4,
+-- page 436.
+volume :: Tetrahedron -> Double
+volume t
+ | (v0 t) == (v1 t) = 0
+ | (v0 t) == (v2 t) = 0
+ | (v0 t) == (v3 t) = 0
+ | (v1 t) == (v2 t) = 0
+ | (v1 t) == (v3 t) = 0
+ | (v2 t) == (v3 t) = 0
+ | otherwise = (1/6)*(det (vol_matrix t))
+
+
+b0 :: Tetrahedron -> (RealFunction Point)
+b0 t point = (volume inner_tetrahedron) / (volume t)
+ where
+ inner_tetrahedron = t { v0 = point }
+
+b1 :: Tetrahedron -> (RealFunction Point)
+b1 t point = (volume inner_tetrahedron) / (volume t)
+ where
+ inner_tetrahedron = t { v1 = point }
+
+b2 :: Tetrahedron -> (RealFunction Point)
+b2 t point = (volume inner_tetrahedron) / (volume t)
+ where
+ inner_tetrahedron = t { v2 = point }
+
+b3 :: Tetrahedron -> (RealFunction Point)
+b3 t point = (volume inner_tetrahedron) / (volume t)
+ where
+ inner_tetrahedron = t { v3 = point }
--- /dev/null
+module ThreeDimensional
+where
+
+import Point(Point)
+
+class ThreeDimensional a where
+ center :: a -> Point
+ contains_point :: a -> Point -> Bool
--- /dev/null
+import Test.HUnit
+import Test.QuickCheck
+
+import Tests.Face
+import Tests.Misc
+import Tests.Tetrahedron
+
+-- The list of HUnit tests.
+test_suite = TestList (concat [misc_tests,
+ tetrahedron_tests])
+
+main :: IO ()
+main = do
+ putStrLn "HUnit"
+ putStrLn "-----"
+ runTestTT test_suite
+ putStrLn ""
+ putStrLn "QuickCheck"
+ putStrLn "----------"
+ let qc_args = stdArgs { maxSuccess = 1000,
+ maxDiscard = 5000,
+ maxSize = 1000 }
+
+ putStr "prop_all_volumes_nonnegative... "
+ quickCheckWith qc_args prop_all_volumes_nonnegative
+
+ putStr "prop_factorial_greater... "
+ quickCheckWith qc_args prop_factorial_greater
+
+ putStr "prop_b0_v0_always_unity... "
+ quickCheckWith qc_args prop_b0_v0_always_unity
+
+ putStr "prop_b0_v1_always_zero... "
+ quickCheckWith qc_args prop_b0_v1_always_zero
+
+ putStr "prop_b0_v2_always_zero... "
+ quickCheckWith qc_args prop_b0_v2_always_zero
+
+ putStr "prop_b0_v3_always_zero... "
+ quickCheckWith qc_args prop_b0_v3_always_zero
+
+ putStr "prop_b1_v1_always_unity... "
+ quickCheckWith qc_args prop_b1_v1_always_unity
+
+ putStr "prop_b1_v0_always_zero... "
+ quickCheckWith qc_args prop_b1_v0_always_zero
+
+ putStr "prop_b1_v2_always_zero... "
+ quickCheckWith qc_args prop_b1_v2_always_zero
+
+ putStr "prop_b1_v3_always_zero... "
+ quickCheckWith qc_args prop_b1_v3_always_zero
+
+ putStr "prop_b2_v2_always_unity... "
+ quickCheckWith qc_args prop_b2_v2_always_unity
+
+ putStr "prop_b2_v0_always_zero... "
+ quickCheckWith qc_args prop_b2_v0_always_zero
+
+ putStr "prop_b2_v1_always_zero... "
+ quickCheckWith qc_args prop_b2_v1_always_zero
+
+ putStr "prop_b2_v3_always_zero... "
+ quickCheckWith qc_args prop_b2_v3_always_zero
+
+ putStr "prop_b3_v3_always_unity... "
+ quickCheckWith qc_args prop_b3_v3_always_unity
+
+ putStr "prop_b3_v0_always_zero... "
+ quickCheckWith qc_args prop_b3_v0_always_zero
+
+ putStr "prop_b3_v1_always_zero... "
+ quickCheckWith qc_args prop_b3_v1_always_zero
+
+ putStr "prop_b3_v2_always_zero... "
+ quickCheckWith qc_args prop_b3_v2_always_zero
+
+ putStr "prop_c3000_identity... "
+ quickCheckWith qc_args prop_c3000_identity
+
+ putStr "prop_c2100_identity... "
+ quickCheckWith qc_args prop_c2100_identity
+
+ putStr "prop_c1110_identity... "
+ quickCheckWith qc_args prop_c1110_identity
+
+ return ()
--- /dev/null
+#!/bin/bash
+#
+# hscolour_srcs
+#
+# Called from the makefile to generate colourised source HTML which is
+# then used by Haddock.
+#
+
+# Find all source files. The way we do this is important, because of
+# the dirname and basename invocations below. The paths here must
+# coincide with the ones passed to 'haddock' in the makefile.
+SRCS=`find ./src -iname '*.hs'`
+
+# We have to create the output directories before we start.
+mkdir -p doc/src/Tests
+
+for file in $SRCS; do
+ DIRNAME=`dirname $file`
+ BASENAME=`basename $file .hs`
+ OUTFILE="doc/${DIRNAME}/${BASENAME}.html"
+ HsColour -html -anchor $file > $OUTFILE
+done
projects / spline3.git / commitdiff