X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FCube.hs;h=6652e8b6331f039cd90cd9c7791dab61f2fd3f58;hb=d9eed953bd810f6928de536617dc21121a8a645b;hp=11762a5efeb0b1366d0455da5e97198f82ece0fc;hpb=71c69c67074e6eb6ce7520739bc729691525b20b;p=spline3.git diff --git a/src/Cube.hs b/src/Cube.hs index 11762a5..6652e8b 100644 --- a/src/Cube.hs +++ b/src/Cube.hs @@ -24,15 +24,13 @@ import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose) import Cardinal import Comparisons ((~=), (~~=)) -import qualified Face (Face(Face, v0, v1, v2, v3)) +import qualified Face (Face(..), center) import FunctionValues (FunctionValues, eval, rotate) import Misc (all_equal, disjoint) import Point (Point(..), dot) -import Tetrahedron (Tetrahedron(..), c, volume) -import ThreeDimensional +import Tetrahedron (Tetrahedron(..), barycenter, c, volume) -data Cube = Cube { h :: !Double, - i :: !Int, +data Cube = Cube { i :: !Int, j :: !Int, k :: !Int, fv :: !FunctionValues, @@ -42,13 +40,12 @@ 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 @@ -61,7 +58,6 @@ instance Arbitrary Cube where 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 +73,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 = Point 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 (Point 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,7 +130,7 @@ 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 cc = center cube v0' = cc + ( Point delta (-delta) delta ) v1' = cc + ( Point delta delta delta ) @@ -167,7 +143,7 @@ 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 cc = center cube v0' = cc + ( Point delta (-delta) (-delta) ) v1' = cc + ( Point delta delta (-delta) ) @@ -179,7 +155,7 @@ back_face cube = Face.Face v0' v1' v2' v3' down_face :: Cube -> Face.Face down_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = 1/2 cc = center cube v0' = cc + ( Point (-delta) (-delta) (-delta) ) v1' = cc + ( Point (-delta) delta (-delta) ) @@ -192,7 +168,7 @@ 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 cc = center cube v0' = cc + ( Point (-delta) (-delta) delta ) v1' = cc + ( Point (-delta) delta delta ) @@ -203,7 +179,7 @@ front_face cube = Face.Face v0' v1' v2' v3' left_face :: Cube -> Face.Face left_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = 1/2 cc = center cube v0' = cc + ( Point delta (-delta) delta ) v1' = cc + ( Point (-delta) (-delta) delta ) @@ -215,7 +191,7 @@ left_face cube = Face.Face v0' v1' v2' v3' right_face :: Cube -> Face.Face right_face cube = Face.Face v0' v1' v2' v3' where - delta = (1/2)*(h cube) + delta = 1/2 cc = center cube v0' = cc + ( Point (-delta) delta delta) v1' = cc + ( Point delta delta delta ) @@ -230,7 +206,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 +216,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 +227,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 +238,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 +249,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 +260,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 +271,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 +284,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 +295,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 +306,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 +319,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 +334,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 +347,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 +358,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 +369,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 +382,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 +393,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 +404,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 +415,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 +428,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 +440,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 +451,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 +462,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 +475,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. @@ -629,6 +602,7 @@ find_containing_tetrahedron cube p = top_half = in_top_half cube p left_half = in_left_half cube p + candidates :: V.Vector Tetrahedron candidates = if front_half then @@ -659,10 +633,20 @@ find_containing_tetrahedron cube p = -- Use the dot product instead of Euclidean distance here to save -- a sqrt(). So, "distances" below really means "distances -- squared." - distances = V.map ((dot p) . center) candidates + 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 @@ -712,14 +696,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