X-Git-Url: https://gitweb.michael.orlitzky.com/?p=spline3.git;a=blobdiff_plain;f=src%2FCube.hs;h=8596b11846838b91f0b4db3667d130ac4152b7ed;hp=32b4b652c4ceebf76f4837fd5ed6d1b2a4723cc0;hb=83ef0aaeae074756e4ee90d72d3e27e74e136061;hpb=e46f39ac6687ba24e5f3ff3d8fef4efe309c40cd diff --git a/src/Cube.hs b/src/Cube.hs index 32b4b65..8596b11 100644 --- a/src/Cube.hs +++ b/src/Cube.hs @@ -3,11 +3,10 @@ module Cube ( cube_properties, find_containing_tetrahedron, tetrahedra, - tetrahedron - ) + tetrahedron ) where -import Data.Maybe (fromJust) +import Data.Maybe ( fromJust ) import qualified Data.Vector as V ( Vector, findIndex, @@ -15,24 +14,36 @@ import qualified Data.Vector as V ( minimum, singleton, snoc, - unsafeIndex - ) -import Prelude hiding (LT) -import Test.Framework (Test, testGroup) -import Test.Framework.Providers.QuickCheck2 (testProperty) -import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose) - -import Cardinal -import Comparisons ((~=), (~~=)) -import qualified Face (Face(Face, v0, v1, v2, v3)) -import FunctionValues (FunctionValues, eval, rotate) -import Misc (all_equal, disjoint) -import Point -import Tetrahedron (Tetrahedron(..), c, volume) -import ThreeDimensional - -data Cube = Cube { h :: !Double, - i :: !Int, + unsafeIndex) +import Prelude hiding ( LT ) +import Test.Tasty ( TestTree, testGroup ) +import Test.Tasty.QuickCheck ( + Arbitrary( arbitrary ), + Gen, + Positive( Positive ), + choose, + testProperty ) +import Cardinal ( + Cardinal(F, B, L, R, D, T, FL, FR, FD, FT, + BL, BR, BD, BT, LD, LT, RD, RT, I), + ccwx, + ccwy, + ccwz, + cwx, + cwy, + cwz ) +import Comparisons ( (~=), (~~=) ) +import qualified Face ( Face(..), center ) +import FunctionValues ( FunctionValues, eval, rotate ) +import Misc ( all_equal, disjoint ) +import Point ( Point( Point ), dot ) +import Tetrahedron ( + Tetrahedron(Tetrahedron, function_values, v0, v1, v2, v3), + barycenter, + c, + volume ) + +data Cube = Cube { i :: !Int, j :: !Int, k :: !Int, fv :: !FunctionValues, @@ -42,26 +53,24 @@ data Cube = Cube { h :: !Double, instance Arbitrary Cube where arbitrary = do - (Positive h') <- arbitrary :: Gen (Positive Double) i' <- choose (coordmin, coordmax) j' <- choose (coordmin, coordmax) k' <- choose (coordmin, coordmax) fv' <- arbitrary :: Gen FunctionValues (Positive tet_vol) <- arbitrary :: Gen (Positive Double) - return (Cube h' i' j' k' fv' tet_vol) + return (Cube i' j' k' fv' tet_vol) 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 + coordmax = 4194304 :: Int -- 2^22 coordmin = -coordmax instance Show Cube where show cube = "Cube_" ++ subscript ++ "\n" ++ - " h: " ++ (show (h cube)) ++ "\n" ++ " Center: " ++ (show (center cube)) ++ "\n" ++ " xmin: " ++ (show (xmin cube)) ++ "\n" ++ " xmax: " ++ (show (xmax cube)) ++ "\n" ++ @@ -77,76 +86,56 @@ instance Show Cube where -- | The left-side boundary of the cube. See Sorokina and Zeilfelder, -- p. 76. xmin :: Cube -> Double -xmin cube = (i' - 1/2)*delta +xmin cube = (i' - 1/2) 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 = (i' + 1/2)*delta +xmax cube = (i' + 1/2) 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 = (j' - 1/2)*delta +ymin cube = (j' - 1/2) 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 = (j' + 1/2)*delta +ymax cube = (j' + 1/2) 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 = (k' - 1/2)*delta +zmin cube = (k' - 1/2) 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 = (k' + 1/2)*delta +zmax cube = (k' + 1/2) where k' = fromIntegral (k cube) :: Double - delta = h cube - -instance ThreeDimensional Cube where - -- | The center of Cube_ijk coincides with v_ijk at - -- (ih, jh, kh). See Sorokina and Zeilfelder, p. 76. - center cube = (x, y, z) - where - delta = h cube - i' = fromIntegral (i cube) :: Double - j' = fromIntegral (j cube) :: Double - k' = fromIntegral (k cube) :: Double - x = delta * i' - y = delta * j' - z = delta * k' - - -- | It's easy to tell if a point is within a cube; just make sure - -- that it falls on the proper side of each of the cube's faces. - contains_point cube (x, y, z) - | x < (xmin cube) = False - | x > (xmax cube) = False - | y < (ymin cube) = False - | y > (ymax cube) = False - | z < (zmin cube) = False - | z > (zmax cube) = False - | otherwise = True +-- | The center of Cube_ijk coincides with v_ijk at +-- (i, j, k). See Sorokina and Zeilfelder, p. 76. +center :: Cube -> Point +center cube = + Point x y z + where + x = fromIntegral (i cube) :: Double + y = fromIntegral (j cube) :: Double + z = fromIntegral (k cube) :: Double + -- Face stuff. @@ -154,12 +143,12 @@ instance ThreeDimensional Cube where top_face :: Cube -> Face.Face top_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = (1/2) :: Double cc = center cube - v0' = cc + (delta, -delta, delta) - v1' = cc + (delta, delta, delta) - v2' = cc + (-delta, delta, delta) - v3' = cc + (-delta, -delta, delta) + v0' = cc + ( Point delta (-delta) delta ) + v1' = cc + ( Point delta delta delta ) + v2' = cc + ( Point (-delta) delta delta ) + v3' = cc + ( Point (-delta) (-delta) delta ) @@ -167,24 +156,24 @@ top_face cube = Face.Face v0' v1' v2' v3' back_face :: Cube -> Face.Face back_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = (1/2) :: Double cc = center cube - v0' = cc + (delta, -delta, -delta) - v1' = cc + (delta, delta, -delta) - v2' = cc + (delta, delta, delta) - v3' = cc + (delta, -delta, delta) + v0' = cc + ( Point delta (-delta) (-delta) ) + v1' = cc + ( Point delta delta (-delta) ) + v2' = cc + ( Point delta delta delta ) + v3' = cc + ( Point delta (-delta) delta ) -- The bottom face (in the direction of -z) of the cube. down_face :: Cube -> Face.Face down_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = (1/2) :: Double cc = center cube - v0' = cc + (-delta, -delta, -delta) - v1' = cc + (-delta, delta, -delta) - v2' = cc + (delta, delta, -delta) - v3' = cc + (delta, -delta, -delta) + v0' = cc + ( Point (-delta) (-delta) (-delta) ) + v1' = cc + ( Point (-delta) delta (-delta) ) + v2' = cc + ( Point delta delta (-delta) ) + v3' = cc + ( Point delta (-delta) (-delta) ) @@ -192,35 +181,35 @@ down_face cube = Face.Face v0' v1' v2' v3' front_face :: Cube -> Face.Face front_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = (1/2) :: Double cc = center cube - v0' = cc + (-delta, -delta, delta) - v1' = cc + (-delta, delta, delta) - v2' = cc + (-delta, delta, -delta) - v3' = cc + (-delta, -delta, -delta) + v0' = cc + ( Point (-delta) (-delta) delta ) + v1' = cc + ( Point (-delta) delta delta ) + v2' = cc + ( Point (-delta) delta (-delta) ) + v3' = cc + ( Point (-delta) (-delta) (-delta) ) -- | The left (in the direction of -y) face of the cube. left_face :: Cube -> Face.Face left_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = (1/2) :: Double cc = center cube - v0' = cc + (delta, -delta, delta) - v1' = cc + (-delta, -delta, delta) - v2' = cc + (-delta, -delta, -delta) - v3' = cc + (delta, -delta, -delta) + v0' = cc + ( Point delta (-delta) delta ) + v1' = cc + ( Point (-delta) (-delta) delta ) + v2' = cc + ( Point (-delta) (-delta) (-delta) ) + v3' = cc + ( Point delta (-delta) (-delta) ) -- | The right (in the direction of y) face of the cube. right_face :: Cube -> Face.Face right_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = (1/2) :: Double cc = center cube - v0' = cc + (-delta, delta, delta) - v1' = cc + (delta, delta, delta) - v2' = cc + (delta, delta, -delta) - v3' = cc + (-delta, delta, -delta) + v0' = cc + ( Point (-delta) delta delta) + v1' = cc + ( Point delta delta delta ) + v2' = cc + ( Point delta delta (-delta) ) + v3' = cc + ( Point (-delta) delta (-delta) ) tetrahedron :: Cube -> Int -> Tetrahedron @@ -230,7 +219,7 @@ tetrahedron cube 0 = where v0' = center cube ff = front_face cube - v1' = center ff + v1' = Face.center ff v2' = Face.v0 ff v3' = Face.v1 ff vol = tetrahedra_volume cube @@ -240,7 +229,7 @@ tetrahedron cube 1 = where v0' = center cube ff = front_face cube - v1' = center ff + v1' = Face.center ff v2' = Face.v1 ff v3' = Face.v2 ff fv' = rotate ccwx (fv cube) @@ -251,7 +240,7 @@ tetrahedron cube 2 = where v0' = center cube ff = front_face cube - v1' = center ff + v1' = Face.center ff v2' = Face.v2 ff v3' = Face.v3 ff fv' = rotate ccwx $ rotate ccwx $ fv cube @@ -262,7 +251,7 @@ tetrahedron cube 3 = where v0' = center cube ff = front_face cube - v1' = center ff + v1' = Face.center ff v2' = Face.v3 ff v3' = Face.v0 ff fv' = rotate cwx (fv cube) @@ -273,7 +262,7 @@ tetrahedron cube 4 = where v0' = center cube tf = top_face cube - v1' = center tf + v1' = Face.center tf v2' = Face.v0 tf v3' = Face.v1 tf fv' = rotate cwy (fv cube) @@ -284,7 +273,7 @@ tetrahedron cube 5 = where v0' = center cube tf = top_face cube - v1' = center tf + v1' = Face.center tf v2' = Face.v1 tf v3' = Face.v2 tf fv' = rotate cwy $ rotate cwz $ fv cube @@ -295,7 +284,7 @@ tetrahedron cube 6 = where v0' = center cube tf = top_face cube - v1' = center tf + v1' = Face.center tf v2' = Face.v2 tf v3' = Face.v3 tf fv' = rotate cwy $ rotate cwz @@ -308,7 +297,7 @@ tetrahedron cube 7 = where v0' = center cube tf = top_face cube - v1' = center tf + v1' = Face.center tf v2' = Face.v3 tf v3' = Face.v0 tf fv' = rotate cwy $ rotate ccwz $ fv cube @@ -319,7 +308,7 @@ tetrahedron cube 8 = where v0' = center cube bf = back_face cube - v1' = center bf + v1' = Face.center bf v2' = Face.v0 bf v3' = Face.v1 bf fv' = rotate cwy $ rotate cwy $ fv cube @@ -330,7 +319,7 @@ tetrahedron cube 9 = where v0' = center cube bf = back_face cube - v1' = center bf + v1' = Face.center bf v2' = Face.v1 bf v3' = Face.v2 bf fv' = rotate cwy $ rotate cwy @@ -343,7 +332,7 @@ tetrahedron cube 10 = where v0' = center cube bf = back_face cube - v1' = center bf + v1' = Face.center bf v2' = Face.v2 bf v3' = Face.v3 bf fv' = rotate cwy $ rotate cwy @@ -358,7 +347,7 @@ tetrahedron cube 11 = where v0' = center cube bf = back_face cube - v1' = center bf + v1' = Face.center bf v2' = Face.v3 bf v3' = Face.v0 bf fv' = rotate cwy $ rotate cwy @@ -371,7 +360,7 @@ tetrahedron cube 12 = where v0' = center cube df = down_face cube - v1' = center df + v1' = Face.center df v2' = Face.v0 df v3' = Face.v1 df fv' = rotate ccwy $ fv cube @@ -382,7 +371,7 @@ tetrahedron cube 13 = where v0' = center cube df = down_face cube - v1' = center df + v1' = Face.center df v2' = Face.v1 df v3' = Face.v2 df fv' = rotate ccwy $ rotate ccwz $ fv cube @@ -393,7 +382,7 @@ tetrahedron cube 14 = where v0' = center cube df = down_face cube - v1' = center df + v1' = Face.center df v2' = Face.v2 df v3' = Face.v3 df fv' = rotate ccwy $ rotate ccwz @@ -406,7 +395,7 @@ tetrahedron cube 15 = where v0' = center cube df = down_face cube - v1' = center df + v1' = Face.center df v2' = Face.v3 df v3' = Face.v0 df fv' = rotate ccwy $ rotate cwz $ fv cube @@ -417,7 +406,7 @@ tetrahedron cube 16 = where v0' = center cube rf = right_face cube - v1' = center rf + v1' = Face.center rf v2' = Face.v0 rf v3' = Face.v1 rf fv' = rotate ccwz $ fv cube @@ -428,7 +417,7 @@ tetrahedron cube 17 = where v0' = center cube rf = right_face cube - v1' = center rf + v1' = Face.center rf v2' = Face.v1 rf v3' = Face.v2 rf fv' = rotate ccwz $ rotate cwy $ fv cube @@ -439,7 +428,7 @@ tetrahedron cube 18 = where v0' = center cube rf = right_face cube - v1' = center rf + v1' = Face.center rf v2' = Face.v2 rf v3' = Face.v3 rf fv' = rotate ccwz $ rotate cwy @@ -452,7 +441,7 @@ tetrahedron cube 19 = where v0' = center cube rf = right_face cube - v1' = center rf + v1' = Face.center rf v2' = Face.v3 rf v3' = Face.v0 rf fv' = rotate ccwz $ rotate ccwy @@ -464,7 +453,7 @@ tetrahedron cube 20 = where v0' = center cube lf = left_face cube - v1' = center lf + v1' = Face.center lf v2' = Face.v0 lf v3' = Face.v1 lf fv' = rotate cwz $ fv cube @@ -475,7 +464,7 @@ tetrahedron cube 21 = where v0' = center cube lf = left_face cube - v1' = center lf + v1' = Face.center lf v2' = Face.v1 lf v3' = Face.v2 lf fv' = rotate cwz $ rotate ccwy $ fv cube @@ -486,7 +475,7 @@ tetrahedron cube 22 = where v0' = center cube lf = left_face cube - v1' = center lf + v1' = Face.center lf v2' = Face.v2 lf v3' = Face.v3 lf fv' = rotate cwz $ rotate ccwy @@ -499,16 +488,13 @@ tetrahedron cube 23 = where v0' = center cube lf = left_face cube - v1' = center lf + v1' = Face.center lf v2' = Face.v3 lf v3' = Face.v0 lf fv' = rotate cwz $ rotate cwy $ fv cube vol = tetrahedra_volume cube --- Feels dirty, but whatever. -tetrahedron _ _ = error "asked for a nonexistent tetrahedron" - -- Only used in tests, so we don't need the added speed -- of Data.Vector. @@ -588,14 +574,14 @@ back_right_down_tetrahedra cube = (tetrahedron cube 18) in_top_half :: Cube -> Point -> Bool -in_top_half cube (_,_,z) = +in_top_half cube (Point _ _ z) = distance_from_top <= distance_from_bottom where distance_from_top = abs $ (zmax cube) - z distance_from_bottom = abs $ (zmin cube) - z in_front_half :: Cube -> Point -> Bool -in_front_half cube (x,_,_) = +in_front_half cube (Point x _ _) = distance_from_front <= distance_from_back where distance_from_front = abs $ (xmin cube) - x @@ -603,7 +589,7 @@ in_front_half cube (x,_,_) = in_left_half :: Cube -> Point -> Bool -in_left_half cube (_,y,_) = +in_left_half cube (Point _ y _) = distance_from_left <= distance_from_right where distance_from_left = abs $ (ymin cube) - y @@ -620,6 +606,7 @@ in_left_half cube (_,y,_) = -- This can throw an exception, but the use of 'head' might -- save us some unnecessary computations. -- +{-# INLINE find_containing_tetrahedron #-} find_containing_tetrahedron :: Cube -> Point -> Tetrahedron find_containing_tetrahedron cube p = candidates `V.unsafeIndex` (fromJust lucky_idx) @@ -628,9 +615,9 @@ find_containing_tetrahedron cube p = top_half = in_top_half cube p left_half = in_left_half cube p - candidates = - if front_half then - + candidates :: V.Vector Tetrahedron + candidates + | front_half = if left_half then if top_half then front_left_top_tetrahedra cube @@ -642,8 +629,7 @@ find_containing_tetrahedron cube p = else front_right_down_tetrahedra cube - else -- bottom half - + | otherwise = -- back half if left_half then if top_half then back_left_top_tetrahedra cube @@ -655,12 +641,23 @@ find_containing_tetrahedron cube p = else back_right_down_tetrahedra cube - -- Use the dot product instead of 'distance' here to save a - -- sqrt(). So, "distances" below really means "distances squared." - distances = V.map ((dot p) . center) candidates + -- Use the dot product instead of Euclidean distance here to save + -- a sqrt(). So, "distances" below really means "distances + -- squared." + distances :: V.Vector Double + distances = V.map ((dot p) . barycenter) candidates + + shortest_distance :: Double shortest_distance = V.minimum distances + + -- Compute the index of the tetrahedron with the center closest to + -- p. This is a bad algorithm, but don't change it! If you make it + -- smarter by finding the index of shortest_distance in distances + -- (this should give the same answer and avoids recomputing the + -- dot product), the program gets slower. Seriously! + lucky_idx :: Maybe Int lucky_idx = V.findIndex - (\t -> (center t) `dot` p == shortest_distance) + (\t -> (barycenter t) `dot` p == shortest_distance) candidates @@ -668,9 +665,7 @@ find_containing_tetrahedron cube p = --- Tests - --- Quickcheck tests. +-- * Tests prop_opposite_octant_tetrahedra_disjoint1 :: Cube -> Bool prop_opposite_octant_tetrahedra_disjoint1 cube = @@ -710,14 +705,13 @@ prop_all_volumes_positive cube = -- | In fact, since all of the tetrahedra are identical, we should -- already know their volumes. There's 24 tetrahedra to a cube, so --- we'd expect the volume of each one to be (1/24)*h^3. +-- we'd expect the volume of each one to be 1/24. prop_all_volumes_exact :: Cube -> Bool prop_all_volumes_exact cube = - and [volume t ~~= (1/24)*(delta^(3::Int)) | t <- tetrahedra cube] - where - delta = h cube + and [volume t ~~= 1/24 | t <- tetrahedra cube] --- | All tetrahedron should have their v0 located at the center of the cube. +-- | All tetrahedron should have their v0 located at the center of the +-- cube. prop_v0_all_equal :: Cube -> Bool prop_v0_all_equal cube = (v0 t0) == (v0 t1) where @@ -991,7 +985,7 @@ prop_interior_values_all_identical cube = -- This test checks the rotation works as expected. prop_c_tilde_2100_rotation_correct :: Cube -> Bool prop_c_tilde_2100_rotation_correct cube = - expr1 == expr2 + expr1 ~= expr2 where t0 = tetrahedron cube 0 t6 = tetrahedron cube 6 @@ -1025,7 +1019,7 @@ prop_c_tilde_2100_rotation_correct cube = -- even meaningful! prop_c_tilde_2100_correct :: Cube -> Bool prop_c_tilde_2100_correct cube = - c t6 2 1 0 0 == expected + c t6 2 1 0 0 ~= expected where t0 = tetrahedron cube 0 t6 = tetrahedron cube 6 @@ -1147,9 +1141,9 @@ prop_t7_shares_edge_with_t20 cube = t20 = tetrahedron cube 20 -p79_26_properties :: Test.Framework.Test +p79_26_properties :: TestTree p79_26_properties = - testGroup "p. 79, Section (2.6) Properties" [ + testGroup "p. 79, Section (2.6) properties" [ testProperty "c0120 identity1" prop_c0120_identity1, testProperty "c0120 identity2" prop_c0120_identity2, testProperty "c0120 identity3" prop_c0120_identity3, @@ -1163,9 +1157,9 @@ p79_26_properties = testProperty "c1200 identity1" prop_c1200_identity1, testProperty "c2100 identity1" prop_c2100_identity1] -p79_27_properties :: Test.Framework.Test +p79_27_properties :: TestTree p79_27_properties = - testGroup "p. 79, Section (2.7) Properties" [ + testGroup "p. 79, Section (2.7) properties" [ testProperty "c0102 identity1" prop_c0102_identity1, testProperty "c0201 identity1" prop_c0201_identity1, testProperty "c0300 identity2" prop_c0300_identity2, @@ -1174,9 +1168,9 @@ p79_27_properties = testProperty "c2100 identity2" prop_c2100_identity2 ] -p79_28_properties :: Test.Framework.Test +p79_28_properties :: TestTree p79_28_properties = - testGroup "p. 79, Section (2.8) Properties" [ + testGroup "p. 79, Section (2.8) properties" [ testProperty "c3000 identity" prop_c3000_identity, testProperty "c2010 identity" prop_c2010_identity, testProperty "c2001 identity" prop_c2001_identity, @@ -1185,9 +1179,9 @@ p79_28_properties = testProperty "c1011 identity" prop_c1011_identity ] -edge_incidence_tests :: Test.Framework.Test +edge_incidence_tests :: TestTree edge_incidence_tests = - testGroup "Edge Incidence Tests" [ + testGroup "Edge incidence tests" [ testProperty "t0 shares edge with t6" prop_t0_shares_edge_with_t6, testProperty "t0 shares edge with t1" prop_t0_shares_edge_with_t1, testProperty "t0 shares edge with t3" prop_t0_shares_edge_with_t3, @@ -1204,9 +1198,9 @@ edge_incidence_tests = testProperty "t6 shares edge with t7" prop_t6_shares_edge_with_t7, testProperty "t7 shares edge with t20" prop_t7_shares_edge_with_t20 ] -cube_properties :: Test.Framework.Test +cube_properties :: TestTree cube_properties = - testGroup "Cube Properties" [ + testGroup "Cube properties" [ p79_26_properties, p79_27_properties, p79_28_properties,