X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FTetrahedron.hs;h=9a3e4950a1a984dd527e1fe034e1038b38f1ac33;hb=HEAD;hp=0e2f9029dd80f9c42b32a8779cbcf40cf50c6df2;hpb=4d82669d840c49e162f1101ddd9a25c5f3234f92;p=spline3.git diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs index 0e2f902..9a3e495 100644 --- a/src/Tetrahedron.hs +++ b/src/Tetrahedron.hs @@ -1,4 +1,9 @@ +-- The local "coefficient" function defined within the "c" function +-- pattern matches on a bunch of integers, but doesn't handle the +-- "otherwise" case for performance reasons. +{-# OPTIONS_GHC -Wno-incomplete-patterns #-} {-# LANGUAGE BangPatterns #-} + module Tetrahedron ( Tetrahedron(..), b0, -- Cube test @@ -10,27 +15,31 @@ module Tetrahedron ( polynomial, tetrahedron_properties, tetrahedron_tests, - volume -- Cube test - ) + volume ) -- Cube test where -import qualified Data.Vector as V ( - singleton, - snoc, - sum - ) - -import Test.Framework (Test, testGroup) -import Test.Framework.Providers.HUnit (testCase) -import Test.Framework.Providers.QuickCheck2 (testProperty) -import Test.HUnit (Assertion, assertEqual) -import Test.QuickCheck (Arbitrary(..), Gen, Property, (==>)) - -import Comparisons ((~=)) -import FunctionValues (FunctionValues(..), empty_values) -import Misc (factorial) -import Point (Point(..), scale) -import RealFunction (RealFunction, cmult, fexp) +import Data.Vector ( singleton, snoc ) +import qualified Data.Vector as V ( sum ) +import Test.Tasty ( TestTree, testGroup ) +import Test.Tasty.HUnit ( Assertion, assertEqual, testCase ) +import Test.Tasty.QuickCheck ( + Arbitrary( arbitrary ), + Gen, + Property, + (==>), + testProperty ) + +import Comparisons ( (~=) ) +import FunctionValues ( + FunctionValues( front, back, left, right, top, down, front_left, + front_right, front_down, front_top, back_left, back_right, + back_down, back_top, left_down, left_top, right_down, + right_top, front_left_down, front_left_top, + front_right_down, front_right_top, interior ), + empty_values ) +import Misc ( factorial ) +import Point ( Point(Point), scale ) +import RealFunction ( RealFunction, cmult, fexp ) data Tetrahedron = Tetrahedron { function_values :: FunctionValues, @@ -71,32 +80,32 @@ instance Show Tetrahedron where -- We just average the four vertices. barycenter :: Tetrahedron -> Point barycenter (Tetrahedron _ v0' v1' v2' v3' _) = - (v0' + v1' + v2' + v3') `scale` (1/4) + (v0' + v1' + v2' + v3') `scale` (1 / 4) {-# INLINE polynomial #-} polynomial :: Tetrahedron -> (RealFunction Point) polynomial t = - V.sum $ V.singleton ((c t 0 0 0 3) `cmult` (beta t 0 0 0 3)) `V.snoc` - ((c t 0 0 1 2) `cmult` (beta t 0 0 1 2)) `V.snoc` - ((c t 0 0 2 1) `cmult` (beta t 0 0 2 1)) `V.snoc` - ((c t 0 0 3 0) `cmult` (beta t 0 0 3 0)) `V.snoc` - ((c t 0 1 0 2) `cmult` (beta t 0 1 0 2)) `V.snoc` - ((c t 0 1 1 1) `cmult` (beta t 0 1 1 1)) `V.snoc` - ((c t 0 1 2 0) `cmult` (beta t 0 1 2 0)) `V.snoc` - ((c t 0 2 0 1) `cmult` (beta t 0 2 0 1)) `V.snoc` - ((c t 0 2 1 0) `cmult` (beta t 0 2 1 0)) `V.snoc` - ((c t 0 3 0 0) `cmult` (beta t 0 3 0 0)) `V.snoc` - ((c t 1 0 0 2) `cmult` (beta t 1 0 0 2)) `V.snoc` - ((c t 1 0 1 1) `cmult` (beta t 1 0 1 1)) `V.snoc` - ((c t 1 0 2 0) `cmult` (beta t 1 0 2 0)) `V.snoc` - ((c t 1 1 0 1) `cmult` (beta t 1 1 0 1)) `V.snoc` - ((c t 1 1 1 0) `cmult` (beta t 1 1 1 0)) `V.snoc` - ((c t 1 2 0 0) `cmult` (beta t 1 2 0 0)) `V.snoc` - ((c t 2 0 0 1) `cmult` (beta t 2 0 0 1)) `V.snoc` - ((c t 2 0 1 0) `cmult` (beta t 2 0 1 0)) `V.snoc` - ((c t 2 1 0 0) `cmult` (beta t 2 1 0 0)) `V.snoc` + V.sum $ singleton ((c t 0 0 0 3) `cmult` (beta t 0 0 0 3)) `snoc` + ((c t 0 0 1 2) `cmult` (beta t 0 0 1 2)) `snoc` + ((c t 0 0 2 1) `cmult` (beta t 0 0 2 1)) `snoc` + ((c t 0 0 3 0) `cmult` (beta t 0 0 3 0)) `snoc` + ((c t 0 1 0 2) `cmult` (beta t 0 1 0 2)) `snoc` + ((c t 0 1 1 1) `cmult` (beta t 0 1 1 1)) `snoc` + ((c t 0 1 2 0) `cmult` (beta t 0 1 2 0)) `snoc` + ((c t 0 2 0 1) `cmult` (beta t 0 2 0 1)) `snoc` + ((c t 0 2 1 0) `cmult` (beta t 0 2 1 0)) `snoc` + ((c t 0 3 0 0) `cmult` (beta t 0 3 0 0)) `snoc` + ((c t 1 0 0 2) `cmult` (beta t 1 0 0 2)) `snoc` + ((c t 1 0 1 1) `cmult` (beta t 1 0 1 1)) `snoc` + ((c t 1 0 2 0) `cmult` (beta t 1 0 2 0)) `snoc` + ((c t 1 1 0 1) `cmult` (beta t 1 1 0 1)) `snoc` + ((c t 1 1 1 0) `cmult` (beta t 1 1 1 0)) `snoc` + ((c t 1 2 0 0) `cmult` (beta t 1 2 0 0)) `snoc` + ((c t 2 0 0 1) `cmult` (beta t 2 0 0 1)) `snoc` + ((c t 2 0 1 0) `cmult` (beta t 2 0 1 0)) `snoc` + ((c t 2 1 0 0) `cmult` (beta t 2 1 0 0)) `snoc` ((c t 3 0 0 0) `cmult` (beta t 3 0 0 0)) @@ -108,7 +117,7 @@ beta t i j k l = coefficient `cmult` (b0_term * b1_term * b2_term * b3_term) where denominator = (factorial i)*(factorial j)*(factorial k)*(factorial l) - coefficient = 6 / (fromIntegral denominator) + coefficient = (6 / (fromIntegral denominator)) :: Double b0_term = (b0 t) `fexp` i b1_term = (b1 t) `fexp` j b2_term = (b2 t) `fexp` k @@ -151,102 +160,102 @@ c !t !i !j !k !l = coefficient :: Int -> Int -> Int -> Int -> Double coefficient 0 0 3 0 = - (1/8) * (i' + f + l' + t' + lt + fl + ft + flt) + (1 / 8) * (i' + f + l' + t' + lt + fl + ft + flt) coefficient 0 0 0 3 = - (1/8) * (i' + f + r + t' + rt + fr + ft + frt) + (1 / 8) * (i' + f + r + t' + rt + fr + ft + frt) coefficient 0 0 2 1 = - (5/24)*(i' + f + t' + ft) + (1/24)*(l' + fl + lt + flt) + (5 / 24)*(i' + f + t' + ft) + (1 / 24)*(l' + fl + lt + flt) coefficient 0 0 1 2 = - (5/24)*(i' + f + t' + ft) + (1/24)*(r + fr + rt + frt) + (5 / 24)*(i' + f + t' + ft) + (1 / 24)*(r + fr + rt + frt) coefficient 0 1 2 0 = - (5/24)*(i' + f) + (1/8)*(l' + t' + fl + ft) - + (1/24)*(lt + flt) + (5 / 24)*(i' + f) + (1 / 8)*(l' + t' + fl + ft) + + (1 / 24)*(lt + flt) coefficient 0 1 0 2 = - (5/24)*(i' + f) + (1/8)*(r + t' + fr + ft) - + (1/24)*(rt + frt) + (5 / 24)*(i' + f) + (1 / 8)*(r + t' + fr + ft) + + (1 / 24)*(rt + frt) coefficient 0 1 1 1 = - (13/48)*(i' + f) + (7/48)*(t' + ft) - + (1/32)*(l' + r + fl + fr) - + (1/96)*(lt + rt + flt + frt) + (13 / 48)*(i' + f) + (7 / 48)*(t' + ft) + + (1 / 32)*(l' + r + fl + fr) + + (1 / 96)*(lt + rt + flt + frt) coefficient 0 2 1 0 = - (13/48)*(i' + f) + (17/192)*(l' + t' + fl + ft) - + (1/96)*(lt + flt) - + (1/64)*(r + d + fr + fd) - + (1/192)*(rt + ld + frt + fld) + (13 / 48)*(i' + f) + (17 / 192)*(l' + t' + fl + ft) + + (1 / 96)*(lt + flt) + + (1 / 64)*(r + d + fr + fd) + + (1 / 192)*(rt + ld + frt + fld) coefficient 0 2 0 1 = - (13/48)*(i' + f) + (17/192)*(r + t' + fr + ft) - + (1/96)*(rt + frt) - + (1/64)*(l' + d + fl + fd) - + (1/192)*(rd + lt + flt + frd) + (13 / 48)*(i' + f) + (17 / 192)*(r + t' + fr + ft) + + (1 / 96)*(rt + frt) + + (1 / 64)*(l' + d + fl + fd) + + (1 / 192)*(rd + lt + flt + frd) coefficient 0 3 0 0 = - (13/48)*(i' + f) + (5/96)*(l' + r + t' + d + fl + fr + ft + fd) - + (1/192)*(rt + rd + lt + ld + frt + frd + flt + fld) + (13 / 48)*(i' + f) + (5 / 96)*(l' + r + t' + d + fl + fr + ft + fd) + + (1 / 192)*(rt + rd + lt + ld + frt + frd + flt + fld) coefficient 1 0 2 0 = - (1/4)*i' + (1/6)*(f + l' + t') - + (1/12)*(lt + fl + ft) + (1 / 4)*i' + (1 / 6)*(f + l' + t') + + (1 / 12)*(lt + fl + ft) coefficient 1 0 0 2 = - (1/4)*i' + (1/6)*(f + r + t') - + (1/12)*(rt + fr + ft) + (1 / 4)*i' + (1 / 6)*(f + r + t') + + (1 / 12)*(rt + fr + ft) coefficient 1 0 1 1 = - (1/3)*i' + (5/24)*(f + t') - + (1/12)*ft - + (1/24)*(l' + r) - + (1/48)*(lt + rt + fl + fr) + (1 / 3)*i' + (5 / 24)*(f + t') + + (1 / 12)*ft + + (1 / 24)*(l' + r) + + (1 / 48)*(lt + rt + fl + fr) coefficient 1 1 1 0 = - (1/3)*i' + (5/24)*f - + (1/8)*(l' + t') - + (5/96)*(fl + ft) - + (1/48)*(d + r + lt) - + (1/96)*(fd + ld + rt + fr) + (1 / 3)*i' + (5 / 24)*f + + (1 / 8)*(l' + t') + + (5 / 96)*(fl + ft) + + (1 / 48)*(d + r + lt) + + (1 / 96)*(fd + ld + rt + fr) coefficient 1 1 0 1 = - (1/3)*i' + (5/24)*f - + (1/8)*(r + t') - + (5/96)*(fr + ft) - + (1/48)*(d + l' + rt) - + (1/96)*(fd + lt + rd + fl) + (1 / 3)*i' + (5 / 24)*f + + (1 / 8)*(r + t') + + (5 / 96)*(fr + ft) + + (1 / 48)*(d + l' + rt) + + (1 / 96)*(fd + lt + rd + fl) coefficient 1 2 0 0 = - (1/3)*i' + (5/24)*f - + (7/96)*(l' + r + t' + d) - + (1/32)*(fl + fr + ft + fd) - + (1/96)*(rt + rd + lt + ld) + (1 / 3)*i' + (5 / 24)*f + + (7 / 96)*(l' + r + t' + d) + + (1 / 32)*(fl + fr + ft + fd) + + (1 / 96)*(rt + rd + lt + ld) coefficient 2 0 1 0 = - (3/8)*i' + (7/48)*(f + t' + l') - + (1/48)*(r + d + b + lt + fl + ft) - + (1/96)*(rt + bt + fr + fd + ld + bl) + (3 / 8)*i' + (7 / 48)*(f + t' + l') + + (1 / 48)*(r + d + b + lt + fl + ft) + + (1 / 96)*(rt + bt + fr + fd + ld + bl) coefficient 2 0 0 1 = - (3/8)*i' + (7/48)*(f + t' + r) - + (1/48)*(l' + d + b + rt + fr + ft) - + (1/96)*(lt + bt + fl + fd + rd + br) + (3 / 8)*i' + (7 / 48)*(f + t' + r) + + (1 / 48)*(l' + d + b + rt + fr + ft) + + (1 / 96)*(lt + bt + fl + fd + rd + br) coefficient 2 1 0 0 = - (3/8)*i' + (1/12)*(t' + r + l' + d) - + (1/64)*(ft + fr + fl + fd) - + (7/48)*f - + (1/48)*b - + (1/96)*(rt + ld + lt + rd) - + (1/192)*(bt + br + bl + bd) + (3 / 8)*i' + (1 / 12)*(t' + r + l' + d) + + (1 / 64)*(ft + fr + fl + fd) + + (7 / 48)*f + + (1 / 48)*b + + (1 / 96)*(rt + ld + lt + rd) + + (1 / 192)*(bt + br + bl + bd) coefficient 3 0 0 0 = - (3/8)*i' + (1/12)*(t' + f + l' + r + d + b) - + (1/96)*(lt + fl + ft + rt + bt + fr) - + (1/96)*(fd + ld + bd + br + rd + bl) + (3 / 8)*i' + (1 / 12)*(t' + f + l' + r + d + b) + + (1 / 96)*(lt + fl + ft + rt + bt + fr) + + (1 / 96)*(fd + ld + bd + br + rd + bl) @@ -286,7 +295,7 @@ det p0 p1 p2 p3 = {-# INLINE volume #-} volume :: Tetrahedron -> Double volume (Tetrahedron _ v0' v1' v2' v3' _) = - (1/6)*(det v0' v1' v2' v3') + (1 / 6)*(det v0' v1' v2' v3') -- | The barycentric coordinates of a point with respect to v0. {-# INLINE b0 #-} @@ -325,7 +334,7 @@ b3 t point = (volume inner_tetrahedron) / (precomputed_volume t) -- | Check the volume of a particular tetrahedron (computed by hand) -- Its vertices are in clockwise order, so the volume should be -- negative. -tetrahedron1_geometry_tests :: Test.Framework.Test +tetrahedron1_geometry_tests :: TestTree tetrahedron1_geometry_tests = testGroup "tetrahedron1 geometry" [ testCase "volume1" volume1 ] @@ -343,7 +352,7 @@ tetrahedron1_geometry_tests = volume1 :: Assertion volume1 = - assertEqual "volume is correct" True (vol ~= (-1/3)) + assertEqual "volume is correct" True (vol ~= (-1 / 3)) where vol = volume t @@ -351,7 +360,7 @@ tetrahedron1_geometry_tests = -- | Check the volume of a particular tetrahedron (computed by hand) -- Its vertices are in counter-clockwise order, so the volume should -- be positive. -tetrahedron2_geometry_tests :: Test.Framework.Test +tetrahedron2_geometry_tests :: TestTree tetrahedron2_geometry_tests = testGroup "tetrahedron2 geometry" [ testCase "volume1" volume1 ] @@ -368,7 +377,7 @@ tetrahedron2_geometry_tests = precomputed_volume = 0 } volume1 :: Assertion - volume1 = assertEqual "volume1 is correct" True (vol ~= (1/3)) + volume1 = assertEqual "volume1 is correct" True (vol ~= (1 / 3)) where vol = volume t @@ -498,17 +507,17 @@ prop_swapping_vertices_doesnt_affect_coefficients4 t = -tetrahedron_tests :: Test.Framework.Test +tetrahedron_tests :: TestTree tetrahedron_tests = - testGroup "Tetrahedron Tests" [ + testGroup "Tetrahedron tests" [ tetrahedron1_geometry_tests, tetrahedron2_geometry_tests ] -p78_24_properties :: Test.Framework.Test +p78_24_properties :: TestTree p78_24_properties = - testGroup "p. 78, Section (2.4) 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] @@ -516,7 +525,7 @@ p78_24_properties = -- | 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 = - weighted_sum `scale` (1/3) + weighted_sum `scale` (1 / 3) where v0' = (v0 t) `scale` (fromIntegral i) v1' = (v1 t) `scale` (fromIntegral j) @@ -542,10 +551,10 @@ p78_24_properties = (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) + 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) + term4 = (3 / 2)*(p t 1 2 0 0) -- | Given in Sorokina and Zeilfelder, p. 78. prop_c1110_identity :: Tetrahedron -> Property @@ -553,16 +562,16 @@ p78_24_properties = (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)) + 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 +tetrahedron_properties :: TestTree tetrahedron_properties = - testGroup "Tetrahedron Properties" [ + testGroup "Tetrahedron properties" [ p78_24_properties, testProperty "b0_v0_always_unity" prop_b0_v0_always_unity, testProperty "b0_v1_always_zero" prop_b0_v1_always_zero,