| 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)
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.
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.
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.
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.
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.
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
| 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)
+
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
module FunctionValues
where
+import Prelude hiding (LT)
+
import Cardinal
data FunctionValues =
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 }
-- | 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
-- | 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
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 ],
[ 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)
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
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"
projects / spline3.git / commitdiff