X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FTetrahedron.hs;h=a9e6c2c6c299a9ec8df42af3a3b137b5c97728b1;hb=cf7ad91b3a45ea05a8459b9461cdc99210183a20;hp=f3b53198362768b8fdbe83085cc55b8e6306b7cc;hpb=1cd0b90dae4b2a0ea35447427e7962b6fe053308;p=spline3.git diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs index f3b5319..a9e6c2c 100644 --- a/src/Tetrahedron.hs +++ b/src/Tetrahedron.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE BangPatterns #-} module Tetrahedron ( Tetrahedron(..), b0, -- Cube test @@ -18,27 +19,26 @@ import qualified Data.Vector as V ( sum ) -import Prelude hiding (LT) import Test.Framework (Test, testGroup) import Test.Framework.Providers.HUnit (testCase) import Test.Framework.Providers.QuickCheck2 (testProperty) -import Test.HUnit +import Test.HUnit (Assertion, assertEqual) import Test.QuickCheck (Arbitrary(..), Gen, Property, (==>)) import Comparisons ((~=), nearly_ge) import FunctionValues (FunctionValues(..), empty_values) import Misc (factorial) -import Point -import RealFunction -import ThreeDimensional +import Point (Point(..), scale) +import RealFunction (RealFunction, cmult, fexp) +import ThreeDimensional (ThreeDimensional(..)) data Tetrahedron = Tetrahedron { function_values :: FunctionValues, - v0 :: Point, - v1 :: Point, - v2 :: Point, - v3 :: Point, - precomputed_volume :: Double + v0 :: !Point, + v1 :: !Point, + v2 :: !Point, + v3 :: !Point, + precomputed_volume :: !Double } deriving (Eq) @@ -71,6 +71,7 @@ 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 && @@ -120,23 +121,6 @@ polynomial t = ((c t 3 0 0 0) `cmult` (beta t 3 0 0 0)) --- | Returns the domain point of t with indices i,j,k,l. --- Simply an alias for the domain_point function. -xi :: Tetrahedron -> Int -> Int -> Int -> Int -> Point -xi = domain_point - --- | Returns the domain point of t with indices i,j,k,l. -domain_point :: Tetrahedron -> Int -> Int -> Int -> Int -> Point -domain_point t i j k l - | i + j + k + l == 3 = weighted_sum `scale` (1/3) - | otherwise = error "domain point index out of bounds" - where - v0' = (v0 t) `scale` (fromIntegral i) - v1' = (v1 t) `scale` (fromIntegral j) - v2' = (v2 t) `scale` (fromIntegral k) - v3' = (v3 t) `scale` (fromIntegral l) - weighted_sum = v0' + v1' + v2' + v3' - -- | The Bernstein polynomial on t with indices i,j,k,l. Denoted by a -- capital 'B' in the Sorokina/Zeilfelder paper. @@ -160,7 +144,7 @@ beta t i j k l -- Zeilfelder, pp. 84-86. If incorrect indices are supplied, the -- function will simply error. c :: Tetrahedron -> Int -> Int -> Int -> Int -> Double -c t i j k l = +c !t !i !j !k !l = coefficient i j k l where fvs = function_values t @@ -310,10 +294,10 @@ det :: Point -> Point -> Point -> Point -> Double det p0 p1 p2 p3 = term5 + term6 where - (x1, y1, z1) = p0 - (x2, y2, z2) = p1 - (x3, y3, z3) = p2 - (x4, y4, z4) = p3 + Point x1 y1 z1 = p0 + Point x2 y2 z2 = p1 + Point x3 y3 z3 = p2 + Point x4 y4 z4 = p3 term1 = ((x2 - x4)*y1 - (x1 - x4)*y2 + (x1 - x2)*y4)*z3 term2 = ((x2 - x3)*y1 - (x1 - x3)*y2 + (x1 - x2)*y3)*z4 term3 = ((x3 - x4)*y2 - (x2 - x4)*y3 + (x2 - x3)*y4)*z1 @@ -324,6 +308,7 @@ det p0 p1 p2 p3 = -- | Computed using the formula from Lai & Schumaker, Definition 15.4, -- page 436. +{-# INLINE volume #-} volume :: Tetrahedron -> Double volume t | v0' == v1' = 0 @@ -341,6 +326,7 @@ volume t -- | The barycentric coordinates of a point with respect to v0. +{-# INLINE b0 #-} b0 :: Tetrahedron -> (RealFunction Point) b0 t point = (volume inner_tetrahedron) / (precomputed_volume t) where @@ -348,6 +334,7 @@ b0 t point = (volume inner_tetrahedron) / (precomputed_volume t) -- | The barycentric coordinates of a point with respect to v1. +{-# INLINE b1 #-} b1 :: Tetrahedron -> (RealFunction Point) b1 t point = (volume inner_tetrahedron) / (precomputed_volume t) where @@ -355,6 +342,7 @@ b1 t point = (volume inner_tetrahedron) / (precomputed_volume t) -- | The barycentric coordinates of a point with respect to v2. +{-# INLINE b2 #-} b2 :: Tetrahedron -> (RealFunction Point) b2 t point = (volume inner_tetrahedron) / (precomputed_volume t) where @@ -362,6 +350,7 @@ b2 t point = (volume inner_tetrahedron) / (precomputed_volume t) -- | The barycentric coordinates of a point with respect to v3. +{-# INLINE b3 #-} b3 :: Tetrahedron -> (RealFunction Point) b3 t point = (volume inner_tetrahedron) / (precomputed_volume t) where @@ -383,10 +372,10 @@ tetrahedron1_geometry_tests = [ testCase "volume1" volume1, testCase "doesn't contain point1" doesnt_contain_point1] where - p0 = (0, -0.5, 0) - p1 = (0, 0.5, 0) - p2 = (2, 0, 0) - p3 = (1, 0, 1) + p0 = Point 0 (-0.5) 0 + p1 = Point 0 0.5 0 + p2 = Point 2 0 0 + p3 = Point 1 0 1 t = Tetrahedron { v0 = p0, v1 = p1, v2 = p2, @@ -404,7 +393,7 @@ tetrahedron1_geometry_tests = doesnt_contain_point1 = assertEqual "doesn't contain an exterior point" False contained where - exterior_point = (5, 2, -9.0212) + exterior_point = Point 5 2 (-9.0212) contained = contains_point t exterior_point @@ -418,10 +407,10 @@ tetrahedron2_geometry_tests = [ testCase "volume1" volume1, testCase "contains point1" contains_point1] where - p0 = (0, -0.5, 0) - p1 = (2, 0, 0) - p2 = (0, 0.5, 0) - p3 = (1, 0, 1) + p0 = Point 0 (-0.5) 0 + p1 = Point 2 0 0 + p2 = Point 0 0.5 0 + p3 = Point 1 0 1 t = Tetrahedron { v0 = p0, v1 = p1, v2 = p2, @@ -437,7 +426,7 @@ tetrahedron2_geometry_tests = contains_point1 :: Assertion contains_point1 = assertEqual "contains an inner point" True contained where - inner_point = (1, 0, 0.5) + inner_point = Point 1 0 0.5 contained = contains_point t inner_point @@ -451,16 +440,16 @@ containment_tests = testCase "doesn't contain point4" doesnt_contain_point4, testCase "doesn't contain point5" doesnt_contain_point5] where - p2 = (0.5, 0.5, 1) - p3 = (0.5, 0.5, 0.5) - exterior_point = (0, 0, 0) + p2 = Point 0.5 0.5 1 + p3 = Point 0.5 0.5 0.5 + exterior_point = Point 0 0 0 doesnt_contain_point2 :: Assertion doesnt_contain_point2 = assertEqual "doesn't contain an exterior point" False contained where - p0 = (0, 1, 1) - p1 = (1, 1, 1) + p0 = Point 0 1 1 + p1 = Point 1 1 1 t = Tetrahedron { v0 = p0, v1 = p1, v2 = p2, @@ -474,8 +463,8 @@ containment_tests = doesnt_contain_point3 = assertEqual "doesn't contain an exterior point" False contained where - p0 = (1, 1, 1) - p1 = (1, 0, 1) + p0 = Point 1 1 1 + p1 = Point 1 0 1 t = Tetrahedron { v0 = p0, v1 = p1, v2 = p2, @@ -489,8 +478,8 @@ containment_tests = doesnt_contain_point4 = assertEqual "doesn't contain an exterior point" False contained where - p0 = (1, 0, 1) - p1 = (0, 0, 1) + p0 = Point 1 0 1 + p1 = Point 0 0 1 t = Tetrahedron { v0 = p0, v1 = p1, v2 = p2, @@ -504,8 +493,8 @@ containment_tests = doesnt_contain_point5 = assertEqual "doesn't contain an exterior point" False contained where - p0 = (0, 0, 1) - p1 = (0, 1, 1) + p0 = Point 0 0 1 + p1 = Point 0 1 1 t = Tetrahedron { v0 = p0, v1 = p1, v2 = p2, @@ -612,39 +601,6 @@ prop_b3_v2_always_zero t = (volume t) > 0 ==> (b3 t) (v2 t) ~= 0 --- | Used for convenience in the next few tests; not a test itself. -p :: Tetrahedron -> Int -> Int -> Int -> Int -> Double -p t i j k l = (polynomial t) (xi t i j k l) - --- | Given in Sorokina and Zeilfelder, p. 78. -prop_c3000_identity :: Tetrahedron -> Property -prop_c3000_identity t = - (volume t) > 0 ==> - c t 3 0 0 0 ~= p t 3 0 0 0 - --- | Given in Sorokina and Zeilfelder, p. 78. -prop_c2100_identity :: Tetrahedron -> Property -prop_c2100_identity t = - (volume t) > 0 ==> - c t 2 1 0 0 ~= (term1 - term2 + term3 - term4) - where - term1 = (1/3)*(p t 0 3 0 0) - term2 = (5/6)*(p t 3 0 0 0) - term3 = 3*(p t 2 1 0 0) - term4 = (3/2)*(p t 1 2 0 0) - --- | Given in Sorokina and Zeilfelder, p. 78. -prop_c1110_identity :: Tetrahedron -> Property -prop_c1110_identity t = - (volume t) > 0 ==> - c t 1 1 1 0 ~= (term1 + term2 - term3 - term4) - where - term1 = (1/3)*((p t 3 0 0 0) + (p t 0 3 0 0) + (p t 0 0 3 0)) - term2 = (9/2)*(p t 1 1 1 0) - term3 = (3/4)*((p t 2 1 0 0) + (p t 1 2 0 0) + (p t 2 0 1 0)) - term4 = (3/4)*((p t 1 0 2 0) + (p t 0 2 1 0) + (p t 0 1 2 0)) - - prop_swapping_vertices_doesnt_affect_coefficients1 :: Tetrahedron -> Bool prop_swapping_vertices_doesnt_affect_coefficients1 t = c t 0 0 1 2 == c t' 0 0 1 2 @@ -683,10 +639,57 @@ tetrahedron_tests = p78_24_properties :: Test.Framework.Test p78_24_properties = - testGroup "p. 78, Section (2.4) Properties" [ - testProperty "c3000 identity" prop_c3000_identity, - testProperty "c2100 identity" prop_c2100_identity, - testProperty "c1110 identity" prop_c1110_identity] + testGroup "p. 78, Section (2.4) Properties" [ + testProperty "c3000 identity" prop_c3000_identity, + testProperty "c2100 identity" prop_c2100_identity, + testProperty "c1110 identity" prop_c1110_identity] + where + -- | Returns the domain point of t with indices i,j,k,l. + domain_point :: Tetrahedron -> Int -> Int -> Int -> Int -> Point + domain_point t i j k l + | i + j + k + l == 3 = weighted_sum `scale` (1/3) + | otherwise = error "domain point index out of bounds" + where + v0' = (v0 t) `scale` (fromIntegral i) + v1' = (v1 t) `scale` (fromIntegral j) + v2' = (v2 t) `scale` (fromIntegral k) + v3' = (v3 t) `scale` (fromIntegral l) + weighted_sum = v0' + v1' + v2' + v3' + + + -- | Used for convenience in the next few tests. + p :: Tetrahedron -> Int -> Int -> Int -> Int -> Double + p t i j k l = (polynomial t) (domain_point t i j k l) + + + -- | Given in Sorokina and Zeilfelder, p. 78. + prop_c3000_identity :: Tetrahedron -> Property + prop_c3000_identity t = + (volume t) > 0 ==> + c t 3 0 0 0 ~= p t 3 0 0 0 + + -- | Given in Sorokina and Zeilfelder, p. 78. + prop_c2100_identity :: Tetrahedron -> Property + prop_c2100_identity t = + (volume t) > 0 ==> + c t 2 1 0 0 ~= (term1 - term2 + term3 - term4) + where + term1 = (1/3)*(p t 0 3 0 0) + term2 = (5/6)*(p t 3 0 0 0) + term3 = 3*(p t 2 1 0 0) + term4 = (3/2)*(p t 1 2 0 0) + + -- | Given in Sorokina and Zeilfelder, p. 78. + prop_c1110_identity :: Tetrahedron -> Property + prop_c1110_identity t = + (volume t) > 0 ==> + c t 1 1 1 0 ~= (term1 + term2 - term3 - term4) + where + term1 = (1/3)*((p t 3 0 0 0) + (p t 0 3 0 0) + (p t 0 0 3 0)) + term2 = (9/2)*(p t 1 1 1 0) + term3 = (3/4)*((p t 2 1 0 0) + (p t 1 2 0 0) + (p t 2 0 1 0)) + term4 = (3/4)*((p t 1 0 2 0) + (p t 0 2 1 0) + (p t 0 1 2 0)) + tetrahedron_properties :: Test.Framework.Test