From 89b8b6e94fcc944a1f4611811265f3c6217af850 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Tue, 26 Apr 2011 01:06:06 -0400 Subject: [PATCH] Rename the spline project to spline3. --- bin/.gitignore | 1 + doc/.gitignore | 1 + makefile | 49 ++++ src/Cube.hs | 175 +++++++++++++ src/Face.hs | 160 ++++++++++++ src/Grid.hs | 30 +++ src/Main.hs | 104 ++++++++ src/Misc.hs | 20 ++ src/Point.hs | 51 ++++ src/RealFunction.hs | 34 +++ src/Tests/Cube.hs | 20 ++ src/Tests/Face.hs | 19 ++ src/Tests/Grid.hs | 11 + src/Tests/Misc.hs | 23 ++ src/Tests/Tetrahedron.hs | 268 ++++++++++++++++++++ src/Tetrahedron.hs | 531 +++++++++++++++++++++++++++++++++++++++ src/ThreeDimensional.hs | 8 + test/TestSuite.hs | 87 +++++++ util/hscolour_srcs | 22 ++ 19 files changed, 1614 insertions(+) create mode 100644 bin/.gitignore create mode 100644 doc/.gitignore create mode 100644 makefile create mode 100644 src/Cube.hs create mode 100644 src/Face.hs create mode 100644 src/Grid.hs create mode 100644 src/Main.hs create mode 100644 src/Misc.hs create mode 100644 src/Point.hs create mode 100644 src/RealFunction.hs create mode 100644 src/Tests/Cube.hs create mode 100644 src/Tests/Face.hs create mode 100644 src/Tests/Grid.hs create mode 100644 src/Tests/Misc.hs create mode 100644 src/Tests/Tetrahedron.hs create mode 100644 src/Tetrahedron.hs create mode 100644 src/ThreeDimensional.hs create mode 100644 test/TestSuite.hs create mode 100755 util/hscolour_srcs diff --git a/bin/.gitignore b/bin/.gitignore new file mode 100644 index 0000000..80e609a --- /dev/null +++ b/bin/.gitignore @@ -0,0 +1 @@ +[^.]* \ No newline at end of file diff --git a/doc/.gitignore b/doc/.gitignore new file mode 100644 index 0000000..13e4d83 --- /dev/null +++ b/doc/.gitignore @@ -0,0 +1 @@ +[^.]* diff --git a/makefile b/makefile new file mode 100644 index 0000000..cb7e81f --- /dev/null +++ b/makefile @@ -0,0 +1,49 @@ +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 diff --git a/src/Cube.hs b/src/Cube.hs new file mode 100644 index 0000000..d0276c8 --- /dev/null +++ b/src/Cube.hs @@ -0,0 +1,175 @@ +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) diff --git a/src/Face.hs b/src/Face.hs new file mode 100644 index 0000000..d93ebff --- /dev/null +++ b/src/Face.hs @@ -0,0 +1,160 @@ +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 diff --git a/src/Grid.hs b/src/Grid.hs new file mode 100644 index 0000000..e67ac76 --- /dev/null +++ b/src/Grid.hs @@ -0,0 +1,30 @@ +-- | 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 [[[]]] diff --git a/src/Main.hs b/src/Main.hs new file mode 100644 index 0000000..8376ce9 --- /dev/null +++ b/src/Main.hs @@ -0,0 +1,104 @@ +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 diff --git a/src/Misc.hs b/src/Misc.hs new file mode 100644 index 0000000..6364e95 --- /dev/null +++ b/src/Misc.hs @@ -0,0 +1,20 @@ +-- | The Misc module contains helper functions that seem out of place +-- anywhere else. +module Misc +where + + +-- | The standard factorial function. See +-- 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 diff --git a/src/Point.hs b/src/Point.hs new file mode 100644 index 0000000..d0859bf --- /dev/null +++ b/src/Point.hs @@ -0,0 +1,51 @@ +{-# 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) diff --git a/src/RealFunction.hs b/src/RealFunction.hs new file mode 100644 index 0000000..5e0832d --- /dev/null +++ b/src/RealFunction.hs @@ -0,0 +1,34 @@ +{-# 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 diff --git a/src/Tests/Cube.hs b/src/Tests/Cube.hs new file mode 100644 index 0000000..eee5444 --- /dev/null +++ b/src/Tests/Cube.hs @@ -0,0 +1,20 @@ +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) diff --git a/src/Tests/Face.hs b/src/Tests/Face.hs new file mode 100644 index 0000000..14d000b --- /dev/null +++ b/src/Tests/Face.hs @@ -0,0 +1,19 @@ +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 diff --git a/src/Tests/Grid.hs b/src/Tests/Grid.hs new file mode 100644 index 0000000..e22557b --- /dev/null +++ b/src/Tests/Grid.hs @@ -0,0 +1,11 @@ +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) diff --git a/src/Tests/Misc.hs b/src/Tests/Misc.hs new file mode 100644 index 0000000..f282434 --- /dev/null +++ b/src/Tests/Misc.hs @@ -0,0 +1,23 @@ +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 ] diff --git a/src/Tests/Tetrahedron.hs b/src/Tests/Tetrahedron.hs new file mode 100644 index 0000000..98edf63 --- /dev/null +++ b/src/Tests/Tetrahedron.hs @@ -0,0 +1,268 @@ +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)) diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs new file mode 100644 index 0000000..032ce1e --- /dev/null +++ b/src/Tetrahedron.hs @@ -0,0 +1,531 @@ +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 } diff --git a/src/ThreeDimensional.hs b/src/ThreeDimensional.hs new file mode 100644 index 0000000..2775ee9 --- /dev/null +++ b/src/ThreeDimensional.hs @@ -0,0 +1,8 @@ +module ThreeDimensional +where + +import Point(Point) + +class ThreeDimensional a where + center :: a -> Point + contains_point :: a -> Point -> Bool diff --git a/test/TestSuite.hs b/test/TestSuite.hs new file mode 100644 index 0000000..26c6bc9 --- /dev/null +++ b/test/TestSuite.hs @@ -0,0 +1,87 @@ +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 () diff --git a/util/hscolour_srcs b/util/hscolour_srcs new file mode 100755 index 0000000..8bc595f --- /dev/null +++ b/util/hscolour_srcs @@ -0,0 +1,22 @@ +#!/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 -- 2.43.2