From: Michael Orlitzky Date: Sun, 6 Nov 2011 23:16:54 +0000 (-0500) Subject: Drop the no-longer-useful ThreeDimensional class. X-Git-Tag: 0.0.1~71 X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=5f01596d42cca3ec2b8236d697adb468cfcdb055;p=spline3.git Drop the no-longer-useful ThreeDimensional class. --- diff --git a/src/Cube.hs b/src/Cube.hs index 9f93143..d863c29 100644 --- a/src/Cube.hs +++ b/src/Cube.hs @@ -24,12 +24,11 @@ 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, @@ -122,30 +121,20 @@ zmax cube = (k' + 1/2)*delta 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 +-- (ih, jh, kh). See Sorokina and Zeilfelder, p. 76. +center :: Cube -> Point +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' -- Face stuff. @@ -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,7 +488,7 @@ 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 @@ -626,6 +615,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 @@ -656,10 +646,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 @@ -716,7 +716,8 @@ prop_all_volumes_exact cube = where delta = h 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 diff --git a/src/Face.hs b/src/Face.hs index bc316a7..785acaf 100644 --- a/src/Face.hs +++ b/src/Face.hs @@ -1,11 +1,12 @@ -- | The Face module just contains the definition of the 'Face' data -- type and its two typeclass instances. -module Face - ( Face(..) ) +module Face ( + Face(..), + center + ) where import Point -import ThreeDimensional data Face = Face { v0 :: !Point, v1 :: !Point, @@ -14,22 +15,17 @@ data Face = Face { v0 :: !Point, deriving (Eq) instance Show Face where - show f = "Face:\n" ++ - " v0: " ++ (show (v0 f)) ++ "\n" ++ - " v1: " ++ (show (v1 f)) ++ "\n" ++ - " v2: " ++ (show (v2 f)) ++ "\n" ++ - " v3: " ++ (show (v3 f)) ++ "\n" + show (Face v0' v1' v2' v3') = + "Face:\n" ++ + " v0: " ++ (show v0') ++ "\n" ++ + " v1: " ++ (show v1') ++ "\n" ++ + " v2: " ++ (show v2') ++ "\n" ++ + " v3: " ++ (show v3') ++ "\n" --- | The 'Face' type is an instance of 'ThreeDimensional' so that we --- can call the 'center' function on it. This is useful because the --- center of a face is always a vertex of a tetrahedron. -instance ThreeDimensional Face where - -- | Since a face is square, we can just average the four vertices - -- to find the center. - center f = ((v0 f) + (v1 f) + (v2 f) + (v3 f)) `scale` (1/4) - - -- | It's possible to implement this, but it hasn't been done - -- yet. A face will contain a point if the point lies in the same - -- plane as the vertices of the face, and if it falls on the - -- correct side of the four sides of the face. - contains_point _ _ = False +-- | Returns the center of the given face. Since a face is just +-- square, we can average the four vertices to find its center. This +-- is useful because the center of a face is always a vertex of a +-- tetrahedron. +center :: Face -> Point +center (Face v0' v1' v2' v3') = + (v0' + v1' + v2' + v3') `scale` (1/4) diff --git a/src/Grid.hs b/src/Grid.hs index a1058d0..94b9e1d 100644 --- a/src/Grid.hs +++ b/src/Grid.hs @@ -32,8 +32,12 @@ import Examples (trilinear, trilinear9x9x9, zeros, naturals_1d) import FunctionValues (make_values, value_at) import Point (Point(..)) import ScaleFactor (ScaleFactor) -import Tetrahedron (Tetrahedron, c, polynomial, v0, v1, v2, v3) -import ThreeDimensional (ThreeDimensional(..)) +import Tetrahedron ( + Tetrahedron(v0,v1,v2,v3), + c, + contains_point, + polynomial, + ) import Values (Values3D, dims, empty3d, zoom_shape) diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs index a6d69a2..7595729 100644 --- a/src/Tetrahedron.hs +++ b/src/Tetrahedron.hs @@ -5,7 +5,9 @@ module Tetrahedron ( b1, -- Cube test b2, -- Cube test b3, -- Cube test + barycenter, c, + contains_point, polynomial, tetrahedron_properties, tetrahedron_tests, @@ -30,7 +32,6 @@ import FunctionValues (FunctionValues(..), empty_values) import Misc (factorial) import Point (Point(..), scale) import RealFunction (RealFunction, cmult, fexp) -import ThreeDimensional (ThreeDimensional(..)) data Tetrahedron = Tetrahedron { function_values :: FunctionValues, @@ -67,34 +68,40 @@ instance Show Tetrahedron where " v3: " ++ (show (v3 t)) ++ "\n" -instance ThreeDimensional Tetrahedron where - center (Tetrahedron _ v0' v1' v2' v3' _) = - (v0' + v1' + v2' + v3') `scale` (1/4) - - -- contains_point is only used in tests. - contains_point t p0 = - b0_unscaled `nearly_ge` 0 && - b1_unscaled `nearly_ge` 0 && - b2_unscaled `nearly_ge` 0 && - b3_unscaled `nearly_ge` 0 +-- | Find the barycenter of the given tetrahedron. +-- We just average the four vertices. +barycenter :: Tetrahedron -> Point +barycenter (Tetrahedron _ v0' v1' v2' v3' _) = + (v0' + v1' + v2' + v3') `scale` (1/4) + +-- | A point is internal to a tetrahedron if all of its barycentric +-- coordinates with respect to that tetrahedron are non-negative. +contains_point :: Tetrahedron -> Point -> Bool +contains_point t p0 = + b0_unscaled `nearly_ge` 0 && + b1_unscaled `nearly_ge` 0 && + b2_unscaled `nearly_ge` 0 && + b3_unscaled `nearly_ge` 0 + where + -- Drop the useless division and volume calculation that we + -- would do if we used the regular b0,..b3 functions. + b0_unscaled :: Double + b0_unscaled = volume inner_tetrahedron where - -- Drop the useless division and volume calculation that we - -- would do if we used the regular b0,..b3 functions. - b0_unscaled :: Double - b0_unscaled = volume inner_tetrahedron - where inner_tetrahedron = t { v0 = p0 } - - b1_unscaled :: Double - b1_unscaled = volume inner_tetrahedron - where inner_tetrahedron = t { v1 = p0 } - - b2_unscaled :: Double - b2_unscaled = volume inner_tetrahedron - where inner_tetrahedron = t { v2 = p0 } - - b3_unscaled :: Double - b3_unscaled = volume inner_tetrahedron - where inner_tetrahedron = t { v3 = p0 } + inner_tetrahedron = t { v0 = p0 } + + b1_unscaled :: Double + b1_unscaled = volume inner_tetrahedron + where inner_tetrahedron = t { v1 = p0 } + + b2_unscaled :: Double + b2_unscaled = volume inner_tetrahedron + where inner_tetrahedron = t { v2 = p0 } + + b3_unscaled :: Double + b3_unscaled = volume inner_tetrahedron + where inner_tetrahedron = t { v3 = p0 } + {-# INLINE polynomial #-} polynomial :: Tetrahedron -> (RealFunction Point) diff --git a/src/ThreeDimensional.hs b/src/ThreeDimensional.hs deleted file mode 100644 index 4b74f18..0000000 --- a/src/ThreeDimensional.hs +++ /dev/null @@ -1,9 +0,0 @@ -module ThreeDimensional - ( ThreeDimensional(..) ) -where - -import Point(Point) - -class ThreeDimensional a where - center :: a -> Point - contains_point :: a -> Point -> Bool