From 6fb9ab6b6068870323e996da931fc04c7710e3e4 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Tue, 3 May 2011 20:57:37 -0400 Subject: [PATCH] Begin overhauling the program to handle other tetrahedra. Main is pretty much disabled at this point, waiting for me to implement the tetrahedron[0-23] functions. --- src/Cardinal.hs | 37 ++- src/Cube.hs | 204 +++++++++------- src/Face.hs | 205 +++++----------- src/FunctionValues.hs | 60 ++++- src/Grid.hs | 27 ++- src/Main.hs | 91 +++---- src/Tetrahedron.hs | 542 ++++++++++-------------------------------- 7 files changed, 469 insertions(+), 697 deletions(-) diff --git a/src/Cardinal.hs b/src/Cardinal.hs index e4e9464..f62896f 100644 --- a/src/Cardinal.hs +++ b/src/Cardinal.hs @@ -5,19 +5,46 @@ data Cardinal = F | B | L | R - | T | D + | T + | FL + | FR + | FD + | FT + | BL + | BR + | BD + | BT + | LD + | LT + | RD + | RT + | FLD + | FLT + | FRD + | FRT + | BLD + | BLT + | BRD + | BRT + | I + | Scalar Double | Sum Cardinal Cardinal | Difference Cardinal Cardinal | Product Cardinal Cardinal - | ScalarProduct Double Cardinal + | Quotient Cardinal Cardinal deriving (Show, Eq) instance Num Cardinal where x + y = Sum x y x - y = Difference x y - x * y = Sum x y - negate x = ScalarProduct (-1) x + x * y = Product x y + negate x = Product (Scalar (-1)) x abs x = x signum x = x - fromInteger _ = F -- Whatever. + fromInteger x = Scalar (fromIntegral x) + +instance Fractional Cardinal where + x / y = Quotient x y + recip x = Quotient (Scalar 1) x + fromRational x = Scalar (fromRational x) diff --git a/src/Cube.hs b/src/Cube.hs index 3011281..23cbe4d 100644 --- a/src/Cube.hs +++ b/src/Cube.hs @@ -1,59 +1,42 @@ module Cube where -import Grid +import Face +import FunctionValues +--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, +data Cube = Cube { h :: Double, i :: Int, j :: Int, k :: Int, - datum :: Double } + fv :: FunctionValues } 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" + " (zmax: " ++ (show (zmax c)) ++ ")" empty_cube :: Cube -empty_cube = Cube empty_grid 0 0 0 0 +empty_cube = Cube 0 0 0 0 empty_values -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) + delta = h c -- | The right-side boundary of the cube. See Sorokina and Zeilfelder, -- p. 76. @@ -61,7 +44,7 @@ xmax :: Cube -> Double xmax c = (2*i' + 1)*delta / 2 where i' = fromIntegral (i c) :: Double - delta = h (grid c) + delta = h c -- | The front boundary of the cube. See Sorokina and Zeilfelder, -- p. 76. @@ -69,7 +52,7 @@ ymin :: Cube -> Double ymin c = (2*j' - 1)*delta / 2 where j' = fromIntegral (j c) :: Double - delta = h (grid c) + delta = h c -- | The back boundary of the cube. See Sorokina and Zeilfelder, -- p. 76. @@ -77,7 +60,7 @@ ymax :: Cube -> Double ymax c = (2*j' + 1)*delta / 2 where j' = fromIntegral (j c) :: Double - delta = h (grid c) + delta = h c -- | The bottom boundary of the cube. See Sorokina and Zeilfelder, -- p. 76. @@ -85,7 +68,7 @@ zmin :: Cube -> Double zmin c = (2*k' - 1)*delta / 2 where k' = fromIntegral (k c) :: Double - delta = h (grid c) + delta = h c -- | The top boundary of the cube. See Sorokina and Zeilfelder, -- p. 76. @@ -93,14 +76,14 @@ zmax :: Cube -> Double zmax c = (2*k' + 1)*delta / 2 where k' = fromIntegral (k c) :: Double - delta = h (grid c) + delta = h 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) + delta = h c i' = fromIntegral (i c) :: Double j' = fromIntegral (j c) :: Double k' = fromIntegral (k c) :: Double @@ -118,73 +101,120 @@ instance ThreeDimensional Cube where | 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) +-- 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) - (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) +-- 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) +-- signum (Cube g1 i1 j1 k1 d1) = +-- Cube g1 (signum i1) (signum j1) (signum k1) (signum d1) - fromInteger x = empty_cube { datum = (fromIntegral x) } +-- 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) +-- 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) +-- recip (Cube g1 i1 j1 k1 d1) = +-- Cube g1 i1 j1 k1 (recip d1) - fromRational q = empty_cube { datum = fromRational q } +-- fromRational q = empty_cube { datum = fromRational q } --- | Constructs a cube, switching the x and z axes. -reverse_cube :: Grid -> Int -> Int -> Int -> Double -> Cube -reverse_cube g k' j' i' = Cube g i' j' k' -- | Return the cube corresponding to the grid point i,j,k. The list -- of cubes is stored as [z][y][x] but we'll be requesting it by -- [x][y][z] so we flip the indices in the last line. -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 = - (((cubes g) !! k') !! j') !! i' - - --- These next three functions basically form a 'for' loop, looping --- through the xs, ys, and zs in that order. - --- | The cubes_from_values function will return a function that takes --- a list of values and returns a list of cubes. It could just as --- well be written to take the values as a parameter; the omission --- of the last parameter is known as an eta reduce. -cubes_from_values :: Grid -> Int -> Int -> ([Double] -> [Cube]) -cubes_from_values g i' j' = - zipWith (reverse_cube g i' j') [0..] - --- | The cubes_from_planes function will return a function that takes --- a list of planes and returns a list of cubes. It could just as --- well be written to take the planes as a parameter; the omission --- of the last parameter is known as an eta reduce. -cubes_from_planes :: Grid -> Int -> ([[Double]] -> [[Cube]]) -cubes_from_planes g i' = - zipWith (cubes_from_values g i') [0..] - --- | Takes a grid as an argument, and returns a three-dimensional list --- of cubes centered on its grid points. -cubes :: Grid -> [[[Cube]]] -cubes g = zipWith (cubes_from_planes g) [0..] (function_values g) +-- 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 = +-- (((cubes g) !! k') !! j') !! i' + + + + + + +-- Face stuff. + +-- | The top (in the direction of z) face of the cube. +top_face :: Cube -> Face +top_face c = Face v0' v1' v2' v3' + where + delta = (1/2)*(h c) + 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 (in the direction of x) face of the cube. +back_face :: Cube -> Face +back_face c = Face v0' v1' v2' v3' + where + delta = (1/2)*(h c) + 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 (in the direction of -z) of the cube. +down_face :: Cube -> Face +down_face c = Face v0' v1' v2' v3' + where + delta = (1/2)*(h c) + 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 (in the direction of -x) face of the cube. +front_face :: Cube -> Face +front_face c = Face v0' v1' v2' v3' + where + delta = (1/2)*(h c) + 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 (in the direction of -y) face of the cube. +left_face :: Cube -> Face +left_face c = Face v0' v1' v2' v3' + where + delta = (1/2)*(h c) + 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 (in the direction of y) face of the cube. +right_face :: Cube -> Face +right_face c = Face v0' v1' v2' v3' + where + delta = (1/2)*(h c) + v0' = (center c) + (-delta, delta, -delta) + v1' = (center c) + (delta, delta, -delta) + v2' = (center c) + (delta, delta, delta) + v3' = (center c) + (-delta, delta, delta) + diff --git a/src/Face.hs b/src/Face.hs index 6f7e99b..bc1dfee 100644 --- a/src/Face.hs +++ b/src/Face.hs @@ -1,160 +1,83 @@ 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, +data Face = Face { 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) + show f = "Face:\n" ++ + " v0: " ++ (show (v0 f)) ++ "\n" ++ + " v1: " ++ (show (v1 f)) ++ "\n" ++ + " v2: " ++ (show (v2 f)) ++ "\n" ++ + " v3: " ++ (show (v3 f)) ++ "\n" 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 (in the direction of z) 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 back (in the direction of x) 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 (in the direction of -z) 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 front (in the direction of -x) 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 left (in the direction of -y) 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 right (in the direction of y) 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 +-- 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/FunctionValues.hs b/src/FunctionValues.hs index d477a38..e9ff064 100644 --- a/src/FunctionValues.hs +++ b/src/FunctionValues.hs @@ -1,6 +1,8 @@ module FunctionValues where +import Prelude hiding (LT) + import Cardinal data FunctionValues = @@ -44,7 +46,63 @@ eval f L = left f eval f R = right f eval f T = top f eval f D = down f +eval f FL = front_left f +eval f FR = front_right f +eval f FD = front_down f +eval f FT = front_top f +eval f BL = back_left f +eval f BR = back_right f +eval f BD = back_down f +eval f BT = back_top f +eval f LD = left_down f +eval f LT = left_top f +eval f RD = right_down f +eval f RT = right_top f +eval f FLD = front_left_down f +eval f FLT = front_left_top f +eval f FRD = front_right_down f +eval f FRT = front_right_top f +eval f BLD = back_left_down f +eval f BLT = back_left_top f +eval f BRD = back_right_down f +eval f BRT = back_right_top f +eval f I = interior f +eval _ (Scalar x) = x eval f (Sum x y) = (eval f x) + (eval f y) eval f (Difference x y) = (eval f x) - (eval f y) eval f (Product x y) = (eval f x) * (eval f y) -eval f (ScalarProduct x y) = x * (eval f y) +eval f (Quotient x y) = (eval f x) / (eval f y) + +value_at :: [[[Double]]] -> Int -> Int -> Int -> Double +value_at values i j k = + ((values !! k) !! j) !! i + +make_values :: [[[Double]]] -> Int -> Int -> Int -> FunctionValues +make_values values i j k = + empty_values { front = value_at values (i-1) j k, + back = value_at values (i+1) j k, + left = value_at values i (j-1) k, + right = value_at values i (j+1) k, + down = value_at values i j (k-1), + top = value_at values i j (k+1), + front_left = value_at values (i-1) (j-1) k, + front_right = value_at values (i-1) (j+1) k, + front_down =value_at values (i-1) j (k-1), + front_top = value_at values (i-1) j (k+1), + back_left = value_at values (i+1) (j-1) k, + back_right = value_at values (i+1) (j+1) k, + back_down = value_at values (i+1) j (k-1), + back_top = value_at values (i+1) j (k+1), + left_down = value_at values i (j-1) (k-1), + left_top = value_at values i (j-1) (k+1), + right_top = value_at values i (j+1) (k+1), + right_down = value_at values i (j+1) (k-1), + front_left_down = value_at values (i-1) (j-1) (k-1), + front_left_top = value_at values (i-1) (j-1) (k+1), + front_right_down = value_at values (i-1) (j+1) (k-1), + front_right_top = value_at values (i-1) (j+1) (k+1), + back_left_down = value_at values (i-1) (j-1) (k-1), + back_left_top = value_at values (i+1) (j-1) (k+1), + back_right_down = value_at values (i+1) (j+1) (k-1), + back_right_top = value_at values (i+1) (j+1) (k+1), + interior = value_at values i j k } diff --git a/src/Grid.hs b/src/Grid.hs index e67ac76..6d360e2 100644 --- a/src/Grid.hs +++ b/src/Grid.hs @@ -1,13 +1,12 @@ -- | 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) +module Grid where +import Cube (Cube(Cube)) +import FunctionValues + -- | 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 @@ -28,3 +27,21 @@ make_grid grid_size values -- | Creates an empty grid with grid size 1. empty_grid :: Grid empty_grid = Grid 1 [[[]]] + + + +-- This is how we do a 'for' loop in Haskell. +-- No, seriously. +cubes :: Grid -> [[[Cube]]] +cubes g + | fvs == [[[]]] = [[[]]] + | head fvs == [[]] = [[[]]] + | otherwise = + [[[ Cube (h g) i j k (make_values fvs i j k) | i <- [0..xsize]] + | j <- [0..ysize]] + | k <- [0..zsize]] + where + fvs = function_values g + zsize = (length fvs) - 1 + ysize = (length $ head fvs) - 1 + xsize = (length $ head $ head fvs) - 1 diff --git a/src/Main.hs b/src/Main.hs index 4467d17..9ad6199 100644 --- a/src/Main.hs +++ b/src/Main.hs @@ -1,14 +1,14 @@ module Main where -import Cube -import Face -import Grid -import Misc (flatten) -import Point -import RealFunction -import Tetrahedron -import ThreeDimensional +--import Cube +--import Face +--import Grid +--import Misc (flatten) +--import Point +--import RealFunction +--import Tetrahedron +--import ThreeDimensional trilinear :: [[[Double]]] trilinear = [ [ [ 1, 2, 3 ], @@ -48,43 +48,44 @@ dummy = [ [ [ 0, 1, 2 ], [ 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 = head good_cubes - good_tets = filter ((flip contains_point) p) (tetrahedrons target_cube) - target_tetrahedron = head good_tets - poly = polynomial target_tetrahedron +--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 = head good_cubes +-- good_tets = filter ((flip contains_point) p) (tetrahedrons target_cube) +-- target_tetrahedron = head good_tets +-- poly = polynomial target_tetrahedron main :: IO () main = do - print $ find_point_value (0,0,0) - print $ find_point_value (1,0,0) - print $ find_point_value (2,0,0) - print $ find_point_value (0,1,0) - print $ find_point_value (1,1,0) - print $ find_point_value (2,1,0) - print $ find_point_value (0,2,0) - print $ find_point_value (1,2,0) - print $ find_point_value (2,2,0) - print $ find_point_value (0,0,1) - print $ find_point_value (1,0,1) - print $ find_point_value (2,0,1) - print $ find_point_value (0,1,1) - print $ find_point_value (1,1,1) - print $ find_point_value (2,1,1) - print $ find_point_value (0,2,1) - print $ find_point_value (1,2,1) - print $ find_point_value (2,2,1) - print $ find_point_value (0,0,2) - print $ find_point_value (1,0,2) - print $ find_point_value (2,0,2) - print $ find_point_value (0,1,2) - print $ find_point_value (1,1,2) - print $ find_point_value (2,1,2) - print $ find_point_value (0,2,2) - print $ find_point_value (1,2,2) - print $ find_point_value (2,2,2) + putStrLn "Hello, World." + -- print $ find_point_value (0,0,0) + -- print $ find_point_value (1,0,0) + -- print $ find_point_value (2,0,0) + -- print $ find_point_value (0,1,0) + -- print $ find_point_value (1,1,0) + -- print $ find_point_value (2,1,0) + -- print $ find_point_value (0,2,0) + -- print $ find_point_value (1,2,0) + -- print $ find_point_value (2,2,0) + -- print $ find_point_value (0,0,1) + -- print $ find_point_value (1,0,1) + -- print $ find_point_value (2,0,1) + -- print $ find_point_value (0,1,1) + -- print $ find_point_value (1,1,1) + -- print $ find_point_value (2,1,1) + -- print $ find_point_value (0,2,1) + -- print $ find_point_value (1,2,1) + -- print $ find_point_value (2,2,1) + -- print $ find_point_value (0,0,2) + -- print $ find_point_value (1,0,2) + -- print $ find_point_value (2,0,2) + -- print $ find_point_value (0,1,2) + -- print $ find_point_value (1,1,2) + -- print $ find_point_value (2,1,2) + -- print $ find_point_value (0,2,2) + -- print $ find_point_value (1,2,2) + -- print $ find_point_value (2,2,2) diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs index ad164d2..81fe91c 100644 --- a/src/Tetrahedron.hs +++ b/src/Tetrahedron.hs @@ -2,40 +2,29 @@ module Tetrahedron where import Numeric.LinearAlgebra hiding (i, scale) +import Prelude hiding (LT) -import Cube (back, - Cube(datum), - down, - front, - Gridded, - left, - right, - top) +import Cardinal +import FunctionValues import Misc (factorial) import Point import RealFunction import ThreeDimensional -data Tetrahedron = Tetrahedron { cube :: Cube, - v0 :: Point, - v1 :: Point, - v2 :: Point, - v3 :: Point } +data Tetrahedron = Tetrahedron { fv :: FunctionValues, + 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) + show t = "Tetrahedron:\n" ++ + " fv: " ++ (show (fv t)) ++ "\n" ++ + " v0: " ++ (show (v0 t)) ++ "\n" ++ + " v1: " ++ (show (v1 t)) ++ "\n" ++ + " v2: " ++ (show (v2 t)) ++ "\n" ++ + " v3: " ++ (show (v3 t)) ++ "\n" instance ThreeDimensional Tetrahedron where @@ -88,394 +77,121 @@ beta t i j k 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 t 0 0 3 0 = eval (fv t) $ + (1/8) * (I + F + L + T + LT + FL + FT + FLT) + +c t 0 0 0 3 = eval (fv t) $ + (1/8) * (I + F + R + T + RT + FR + FT + FRT) + +c t 0 0 2 1 = eval (fv t) $ + (5/24)*(I + F + T + FT) + + (1/24)*(L + FL + LT + FLT) + +c t 0 0 1 2 = eval (fv t) $ + (5/24)*(I + F + T + FT) + + (1/24)*(R + FR + RT + FRT) + +c t 0 1 2 0 = eval (fv t) $ + (5/24)*(I + F) + + (1/8)*(L + T + FL + FT) + + (1/24)*(LT + FLT) + +c t 0 1 0 2 = eval (fv t) $ + (5/24)*(I + F) + + (1/8)*(R + T + FR + FT) + + (1/24)*(RT + FRT) + +c t 0 1 1 1 = eval (fv t) $ + (13/48)*(I + F) + + (7/48)*(T + FT) + + (1/32)*(L + R + FL + FR) + + (1/96)*(LT + RT + FLT + FRT) + +c t 0 2 1 0 = eval (fv t) $ + (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) + +c t 0 2 0 1 = eval (fv t) $ + (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) + +c t 0 3 0 0 = eval (fv t) $ + (13/48)*(I + F) + + (5/96)*(L + R + T + D + FL + FR + FT + FD) + + (1/192)*(RT + RD + LT + LD + FRT + FRD + FLT + FLD) + +c t 1 0 2 0 = eval (fv t) $ + (1/4)*I + + (1/6)*(F + L + T) + + (1/12)*(LT + FL + FT) + +c t 1 0 0 2 = eval (fv t) $ + (1/4)*I + + (1/6)*(F + R + T) + + (1/12)*(RT + FR + FT) + +c t 1 0 1 1 = eval (fv t) $ + (1/3)*I + + (5/24)*(F + T) + + (1/12)*FT + + (1/24)*(L + R) + + (1/48)*(LT + RT + FL + FR) + +c t 1 1 1 0 = eval (fv t) $ + (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) + +c t 1 1 0 1 = eval (fv t) $ + (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) + +c t 1 2 0 0 = eval (fv t) $ + (1/3)*I + + (5/24)*F + + (7/96)*(L + R + T + D) + + (1/32)*(FL + FR + FT + FD) + + (1/96)*(RT + RD + LT + LD) + +c t 2 0 1 0 = eval (fv t) $ + (3/8)*I + + (7/48)*(F + T + L) + + (1/48)*(R + D + B + LT + FL + FT) + + (1/96)*(RT + BT + FR + FD + LD + BL) + +c t 2 0 0 1 = eval (fv t) $ + (3/8)*I + + (7/48)*(F + T + R) + + (1/48)*(L + D + B + RT + FR + FT) + + (1/96)*(LT + BT + FL + FD + RD + BR) + +c t 2 1 0 0 = eval (fv t) $ + (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) + +c t 3 0 0 0 = eval (fv t) $ + (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) c _ _ _ _ _ = error "coefficient index out of bounds" -- 2.43.2