import Cardinal
import Comparisons ((~=), (~~=))
import qualified Face (Face(Face, v0, v1, v2, v3))
-import FunctionValues
+import FunctionValues (FunctionValues, eval, rotate)
import Misc (all_equal, disjoint)
import Point
-import Tetrahedron (
- Tetrahedron(..),
- c,
- b0,
- b1,
- b2,
- b3,
- volume
- )
+import Tetrahedron (Tetrahedron(..), c, volume)
import ThreeDimensional
-data Cube = Cube { h :: Double,
- i :: Int,
- j :: Int,
- k :: Int,
- fv :: FunctionValues,
- tetrahedra_volume :: Double }
+data Cube = Cube { h :: !Double,
+ i :: !Int,
+ j :: !Int,
+ k :: !Int,
+ fv :: !FunctionValues,
+ tetrahedra_volume :: !Double }
deriving (Eq)
fv' <- arbitrary :: Gen FunctionValues
(Positive tet_vol) <- arbitrary :: Gen (Positive Double)
return (Cube h' i' j' k' fv' tet_vol)
- where
- coordmin = -268435456 -- -(2^29 / 2)
- coordmax = 268435456 -- +(2^29 / 2)
+ where
+ -- The idea here is that, when cubed in the volume formula,
+ -- these numbers don't overflow 64 bits. This number is not
+ -- magic in any other sense than that it does not cause test
+ -- failures, while 2^23 does.
+ coordmax = 4194304 -- 2^22
+ coordmin = -coordmax
instance Show Cube where
" ymin: " ++ (show (ymin cube)) ++ "\n" ++
" ymax: " ++ (show (ymax cube)) ++ "\n" ++
" zmin: " ++ (show (zmin cube)) ++ "\n" ++
- " zmax: " ++ (show (zmax cube)) ++ "\n" ++
- " fv: " ++ (show (Cube.fv cube)) ++ "\n"
+ " zmax: " ++ (show (zmax cube)) ++ "\n"
where
subscript =
(show (i cube)) ++ "," ++ (show (j cube)) ++ "," ++ (show (k cube))
-- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
-- p. 76.
xmin :: Cube -> Double
-xmin cube = (2*i' - 1)*delta / 2
+xmin cube = (i' - 1/2)*delta
where
i' = fromIntegral (i cube) :: Double
delta = h cube
-- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
-- p. 76.
xmax :: Cube -> Double
-xmax cube = (2*i' + 1)*delta / 2
+xmax cube = (i' + 1/2)*delta
where
i' = fromIntegral (i cube) :: Double
delta = h cube
-- | The front boundary of the cube. See Sorokina and Zeilfelder,
-- p. 76.
ymin :: Cube -> Double
-ymin cube = (2*j' - 1)*delta / 2
+ymin cube = (j' - 1/2)*delta
where
j' = fromIntegral (j cube) :: Double
delta = h cube
-- | The back boundary of the cube. See Sorokina and Zeilfelder,
-- p. 76.
ymax :: Cube -> Double
-ymax cube = (2*j' + 1)*delta / 2
+ymax cube = (j' + 1/2)*delta
where
j' = fromIntegral (j cube) :: Double
delta = h cube
-- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
-- p. 76.
zmin :: Cube -> Double
-zmin cube = (2*k' - 1)*delta / 2
+zmin cube = (k' - 1/2)*delta
where
k' = fromIntegral (k cube) :: Double
delta = h cube
-- | The top boundary of the cube. See Sorokina and Zeilfelder,
-- p. 76.
zmax :: Cube -> Double
-zmax cube = (2*k' + 1)*delta / 2
+zmax cube = (k' + 1/2)*delta
where
k' = fromIntegral (k cube) :: Double
delta = h cube
-- | Since the grid size is necessarily positive, all tetrahedra
--- (which comprise cubes of positive volume) must have positive volume
--- as well.
+-- (which comprise cubes of positive volume) must have positive
+-- volume as well.
prop_all_volumes_positive :: Cube -> Bool
prop_all_volumes_positive cube =
- null nonpositive_volumes
+ all (>= 0) volumes
where
ts = tetrahedra cube
volumes = map volume ts
- nonpositive_volumes = filter (<= 0) volumes
+
-- | In fact, since all of the tetrahedra are identical, we should
-- already know their volumes. There's 24 tetrahedra to a cube, so
t6 = tetrahedron cube 6
-
--- | Given in Sorokina and Zeilfelder, p. 78.
-prop_cijk1_identity :: Cube -> Bool
-prop_cijk1_identity cube =
- and [ c t0 i j k 1 ~=
- (c t1 (i+1) j k 0) * ((b0 t0) (v3 t1)) +
- (c t1 i (j+1) k 0) * ((b1 t0) (v3 t1)) +
- (c t1 i j (k+1) 0) * ((b2 t0) (v3 t1)) +
- (c t1 i j k 1) * ((b3 t0) (v3 t1)) | i <- [0..2],
- j <- [0..2],
- k <- [0..2],
- i + j + k == 2]
- where
- t0 = tetrahedron cube 0
- t1 = tetrahedron cube 1
-
-
-- | The function values at the interior should be the same for all
-- tetrahedra.
prop_interior_values_all_identical :: Cube -> Bool
t20 = tetrahedron cube 20
-
-
-
-p78_25_properties :: Test.Framework.Test
-p78_25_properties =
- testGroup "p. 78, Section (2.5) Properties" [
- testProperty "c_ijk1 identity" prop_cijk1_identity ]
-
p79_26_properties :: Test.Framework.Test
p79_26_properties =
testGroup "p. 79, Section (2.6) Properties" [
cube_properties :: Test.Framework.Test
cube_properties =
testGroup "Cube Properties" [
- p78_25_properties,
p79_26_properties,
p79_27_properties,
p79_28_properties,