X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FTetrahedron.hs;h=f3b53198362768b8fdbe83085cc55b8e6306b7cc;hb=1cd0b90dae4b2a0ea35447427e7962b6fe053308;hp=81fe91c11772122d0643b884e4fc3aea52163a43;hpb=6fb9ab6b6068870323e996da931fc04c7710e3e4;p=spline3.git diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs index 81fe91c..f3b5319 100644 --- a/src/Tetrahedron.hs +++ b/src/Tetrahedron.hs @@ -1,26 +1,66 @@ -module Tetrahedron +module Tetrahedron ( + Tetrahedron(..), + b0, -- Cube test + b1, -- Cube test + b2, -- Cube test + b3, -- Cube test + c, + polynomial, + tetrahedron_properties, + tetrahedron_tests, + volume -- Cube test + ) where -import Numeric.LinearAlgebra hiding (i, scale) -import Prelude hiding (LT) +import qualified Data.Vector as V ( + singleton, + snoc, + sum + ) -import Cardinal -import FunctionValues +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.QuickCheck (Arbitrary(..), Gen, Property, (==>)) + +import Comparisons ((~=), nearly_ge) +import FunctionValues (FunctionValues(..), empty_values) import Misc (factorial) import Point import RealFunction import ThreeDimensional -data Tetrahedron = Tetrahedron { fv :: FunctionValues, - v0 :: Point, - v1 :: Point, - v2 :: Point, - v3 :: Point } - deriving (Eq) +data Tetrahedron = + Tetrahedron { function_values :: FunctionValues, + v0 :: Point, + v1 :: Point, + v2 :: Point, + v3 :: Point, + precomputed_volume :: Double + } + deriving (Eq) + + +instance Arbitrary Tetrahedron where + arbitrary = do + rnd_v0 <- arbitrary :: Gen Point + rnd_v1 <- arbitrary :: Gen Point + rnd_v2 <- arbitrary :: Gen Point + rnd_v3 <- arbitrary :: Gen Point + rnd_fv <- arbitrary :: Gen FunctionValues + + -- We can't assign an incorrect precomputed volume, + -- so we have to calculate the correct one here. + let t' = Tetrahedron rnd_fv rnd_v0 rnd_v1 rnd_v2 rnd_v3 0 + let vol = volume t' + return (Tetrahedron rnd_fv rnd_v0 rnd_v1 rnd_v2 rnd_v3 vol) + instance Show Tetrahedron where show t = "Tetrahedron:\n" ++ - " fv: " ++ (show (fv t)) ++ "\n" ++ + " function_values: " ++ (show (function_values t)) ++ "\n" ++ " v0: " ++ (show (v0 t)) ++ "\n" ++ " v1: " ++ (show (v1 t)) ++ "\n" ++ " v2: " ++ (show (v2 t)) ++ "\n" ++ @@ -28,18 +68,56 @@ instance Show Tetrahedron where instance ThreeDimensional Tetrahedron where - center t = ((v0 t) + (v1 t) + (v2 t) + (v3 t)) `scale` (1/4) - contains_point t p = - (b0 t p) >= 0 && (b1 t p) >= 0 && (b2 t p) >= 0 && (b3 t p) >= 0 + center (Tetrahedron _ v0' v1' v2' v3' _) = + (v0' + v1' + v2' + v3') `scale` (1/4) + + 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 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 } polynomial :: Tetrahedron -> (RealFunction Point) polynomial t = - sum [ (c t i j k l) `cmult` (beta t i j k l) | i <- [0..3], - j <- [0..3], - k <- [0..3], - l <- [0..3], - i + j + k + l == 3] + 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` + ((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. @@ -76,176 +154,566 @@ beta t i j k l b3_term = (b3 t) `fexp` l +-- | The coefficient function. c t i j k l returns the coefficient +-- c_ijkl with respect to the tetrahedron t. The definition uses +-- pattern matching to mimic the definitions given in Sorokina and +-- Zeilfelder, pp. 84-86. If incorrect indices are supplied, the +-- function will simply error. c :: Tetrahedron -> Int -> Int -> Int -> Int -> Double -c t 0 0 3 0 = eval (fv t) $ - (1/8) * (I + F + L + T + LT + FL + FT + FLT) - -c t 0 0 0 3 = eval (fv t) $ - (1/8) * (I + F + R + T + RT + FR + FT + FRT) - -c t 0 0 2 1 = eval (fv t) $ - (5/24)*(I + F + T + FT) + - (1/24)*(L + FL + LT + FLT) - -c t 0 0 1 2 = eval (fv t) $ - (5/24)*(I + F + T + FT) + - (1/24)*(R + FR + RT + FRT) - -c t 0 1 2 0 = eval (fv t) $ - (5/24)*(I + F) + - (1/8)*(L + T + FL + FT) + - (1/24)*(LT + FLT) - -c t 0 1 0 2 = eval (fv t) $ - (5/24)*(I + F) + - (1/8)*(R + T + FR + FT) + - (1/24)*(RT + FRT) - -c t 0 1 1 1 = eval (fv t) $ - (13/48)*(I + F) + - (7/48)*(T + FT) + - (1/32)*(L + R + FL + FR) + - (1/96)*(LT + RT + FLT + FRT) - -c t 0 2 1 0 = eval (fv t) $ - (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) - -c t 0 2 0 1 = eval (fv t) $ - (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) - -c t 0 3 0 0 = eval (fv t) $ - (13/48)*(I + F) + - (5/96)*(L + R + T + D + FL + FR + FT + FD) + - (1/192)*(RT + RD + LT + LD + FRT + FRD + FLT + FLD) - -c t 1 0 2 0 = eval (fv t) $ - (1/4)*I + - (1/6)*(F + L + T) + - (1/12)*(LT + FL + FT) - -c t 1 0 0 2 = eval (fv t) $ - (1/4)*I + - (1/6)*(F + R + T) + - (1/12)*(RT + FR + FT) - -c t 1 0 1 1 = eval (fv t) $ - (1/3)*I + - (5/24)*(F + T) + - (1/12)*FT + - (1/24)*(L + R) + - (1/48)*(LT + RT + FL + FR) - -c t 1 1 1 0 = eval (fv t) $ - (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) - -c t 1 1 0 1 = eval (fv t) $ - (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) - -c t 1 2 0 0 = eval (fv t) $ - (1/3)*I + - (5/24)*F + - (7/96)*(L + R + T + D) + - (1/32)*(FL + FR + FT + FD) + - (1/96)*(RT + RD + LT + LD) - -c t 2 0 1 0 = eval (fv t) $ - (3/8)*I + - (7/48)*(F + T + L) + - (1/48)*(R + D + B + LT + FL + FT) + - (1/96)*(RT + BT + FR + FD + LD + BL) - -c t 2 0 0 1 = eval (fv t) $ - (3/8)*I + - (7/48)*(F + T + R) + - (1/48)*(L + D + B + RT + FR + FT) + - (1/96)*(LT + BT + FL + FD + RD + BR) - -c t 2 1 0 0 = eval (fv t) $ - (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) - -c t 3 0 0 0 = eval (fv t) $ - (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) - -c _ _ _ _ _ = error "coefficient index out of bounds" - - - -vol_matrix :: Tetrahedron -> Matrix Double -vol_matrix t = (4><4) $ - [1, 1, 1, 1, - x1, x2, x3, x4, - y1, y2, y3, y4, - z1, z2, z3, z4 ] - where - x1 = x_coord (v0 t) - x2 = x_coord (v1 t) - x3 = x_coord (v2 t) - x4 = x_coord (v3 t) - y1 = y_coord (v0 t) - y2 = y_coord (v1 t) - y3 = y_coord (v2 t) - y4 = y_coord (v3 t) - z1 = z_coord (v0 t) - z2 = z_coord (v1 t) - z3 = z_coord (v2 t) - z4 = z_coord (v3 t) - --- Computed using the formula from Lai & Schumaker, Definition 15.4, --- page 436. +c t i j k l = + coefficient i j k l + where + fvs = function_values t + f = front fvs + b = back fvs + r = right fvs + l' = left fvs + t' = top fvs + d = down fvs + fl = front_left fvs + fr = front_right fvs + fd = front_down fvs + ft = front_top fvs + bl = back_left fvs + br = back_right fvs + bd = back_down fvs + bt = back_top fvs + ld = left_down fvs + lt = left_top fvs + rd = right_down fvs + rt = right_top fvs + fld = front_left_down fvs + flt = front_left_top fvs + frd = front_right_down fvs + frt = front_right_top fvs + i' = interior fvs + + coefficient :: Int -> Int -> Int -> Int -> Double + coefficient 0 0 3 0 = + (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) + + coefficient 0 0 2 1 = + (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) + + coefficient 0 1 2 0 = + (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) + + 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) + + 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) + + 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) + + 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) + + coefficient 1 0 2 0 = + (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) + + 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) + + 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) + + 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) + + 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) + + 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) + + 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) + + 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) + + 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) + + coefficient _ _ _ _ = error "coefficient index out of bounds" + + + +-- | Compute the determinant of the 4x4 matrix, +-- +-- [1] +-- [x] +-- [y] +-- [z] +-- +-- where [1] = [1, 1, 1, 1], +-- [x] = [x1,x2,x3,x4], +-- +-- et cetera. +-- +-- The termX nonsense is an attempt to prevent Double overflow. +-- which has been observed to happen with large coordinates. +-- +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 + 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 + term4 = ((x3 - x4)*y1 - (x1 - x4)*y3 + (x1 - x3)*y4)*z2 + term5 = term1 - term2 + term6 = term3 - term4 + + +-- | Computed using the formula from Lai & Schumaker, Definition 15.4, +-- page 436. volume :: Tetrahedron -> Double volume t - | (v0 t) == (v1 t) = 0 - | (v0 t) == (v2 t) = 0 - | (v0 t) == (v3 t) = 0 - | (v1 t) == (v2 t) = 0 - | (v1 t) == (v3 t) = 0 - | (v2 t) == (v3 t) = 0 - | otherwise = (1/6)*(det (vol_matrix t)) + | v0' == v1' = 0 + | v0' == v2' = 0 + | v0' == v3' = 0 + | v1' == v2' = 0 + | v1' == v3' = 0 + | v2' == v3' = 0 + | otherwise = (1/6)*(det v0' v1' v2' v3') + where + v0' = v0 t + v1' = v1 t + v2' = v2 t + v3' = v3 t +-- | The barycentric coordinates of a point with respect to v0. b0 :: Tetrahedron -> (RealFunction Point) -b0 t point = (volume inner_tetrahedron) / (volume t) +b0 t point = (volume inner_tetrahedron) / (precomputed_volume t) where inner_tetrahedron = t { v0 = point } + +-- | The barycentric coordinates of a point with respect to v1. b1 :: Tetrahedron -> (RealFunction Point) -b1 t point = (volume inner_tetrahedron) / (volume t) +b1 t point = (volume inner_tetrahedron) / (precomputed_volume t) where inner_tetrahedron = t { v1 = point } + +-- | The barycentric coordinates of a point with respect to v2. b2 :: Tetrahedron -> (RealFunction Point) -b2 t point = (volume inner_tetrahedron) / (volume t) +b2 t point = (volume inner_tetrahedron) / (precomputed_volume t) where inner_tetrahedron = t { v2 = point } + +-- | The barycentric coordinates of a point with respect to v3. b3 :: Tetrahedron -> (RealFunction Point) -b3 t point = (volume inner_tetrahedron) / (volume t) +b3 t point = (volume inner_tetrahedron) / (precomputed_volume t) where inner_tetrahedron = t { v3 = point } + + + + +-- Tests + + +-- | Check the volume of a particular tetrahedron (computed by hand) +-- and whether or not it contains a specific point chosen to be +-- outside of it. Its vertices are in clockwise order, so the volume +-- should be negative. +tetrahedron1_geometry_tests :: Test.Framework.Test +tetrahedron1_geometry_tests = + testGroup "tetrahedron1 geometry" + [ 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) + t = Tetrahedron { v0 = p0, + v1 = p1, + v2 = p2, + v3 = p3, + function_values = empty_values, + precomputed_volume = 0 } + + volume1 :: Assertion + volume1 = + assertEqual "volume is correct" True (vol ~= (-1/3)) + where + vol = volume t + + doesnt_contain_point1 :: Assertion + doesnt_contain_point1 = + assertEqual "doesn't contain an exterior point" False contained + where + exterior_point = (5, 2, -9.0212) + contained = contains_point t exterior_point + + +-- | Check the volume of a particular tetrahedron (computed by hand) +-- and whether or not it contains a specific point chosen to be +-- inside of it. Its vertices are in counter-clockwise order, so the +-- volume should be positive. +tetrahedron2_geometry_tests :: Test.Framework.Test +tetrahedron2_geometry_tests = + testGroup "tetrahedron2 geometry" + [ 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) + t = Tetrahedron { v0 = p0, + v1 = p1, + v2 = p2, + v3 = p3, + function_values = empty_values, + precomputed_volume = 0 } + + volume1 :: Assertion + volume1 = assertEqual "volume1 is correct" True (vol ~= (1/3)) + where + vol = volume t + + contains_point1 :: Assertion + contains_point1 = assertEqual "contains an inner point" True contained + where + inner_point = (1, 0, 0.5) + contained = contains_point t inner_point + + +-- | Ensure that tetrahedra do not contain a particular point chosen to +-- be outside of them. +containment_tests :: Test.Framework.Test +containment_tests = + testGroup "containment tests" + [ testCase "doesn't contain point2" doesnt_contain_point2, + testCase "doesn't contain point3" doesnt_contain_point3, + 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) + + 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) + t = Tetrahedron { v0 = p0, + v1 = p1, + v2 = p2, + v3 = p3, + function_values = empty_values, + precomputed_volume = 0 } + contained = contains_point t exterior_point + + + doesnt_contain_point3 :: Assertion + doesnt_contain_point3 = + assertEqual "doesn't contain an exterior point" False contained + where + p0 = (1, 1, 1) + p1 = (1, 0, 1) + t = Tetrahedron { v0 = p0, + v1 = p1, + v2 = p2, + v3 = p3, + function_values = empty_values, + precomputed_volume = 0 } + contained = contains_point t exterior_point + + + doesnt_contain_point4 :: Assertion + doesnt_contain_point4 = + assertEqual "doesn't contain an exterior point" False contained + where + p0 = (1, 0, 1) + p1 = (0, 0, 1) + t = Tetrahedron { v0 = p0, + v1 = p1, + v2 = p2, + v3 = p3, + function_values = empty_values, + precomputed_volume = 0 } + contained = contains_point t exterior_point + + + doesnt_contain_point5 :: Assertion + doesnt_contain_point5 = + assertEqual "doesn't contain an exterior point" False contained + where + p0 = (0, 0, 1) + p1 = (0, 1, 1) + t = Tetrahedron { v0 = p0, + v1 = p1, + v2 = p2, + v3 = p3, + function_values = empty_values, + precomputed_volume = 0 } + contained = contains_point t exterior_point + + +-- | The barycentric coordinate of v0 with respect to itself should +-- be one. +prop_b0_v0_always_unity :: Tetrahedron -> Property +prop_b0_v0_always_unity t = + (volume t) > 0 ==> (b0 t) (v0 t) ~= 1.0 + +-- | The barycentric coordinate of v1 with respect to v0 should +-- be zero. +prop_b0_v1_always_zero :: Tetrahedron -> Property +prop_b0_v1_always_zero t = + (volume t) > 0 ==> (b0 t) (v1 t) ~= 0 + +-- | The barycentric coordinate of v2 with respect to v0 should +-- be zero. +prop_b0_v2_always_zero :: Tetrahedron -> Property +prop_b0_v2_always_zero t = + (volume t) > 0 ==> (b0 t) (v2 t) ~= 0 + +-- | The barycentric coordinate of v3 with respect to v0 should +-- be zero. +prop_b0_v3_always_zero :: Tetrahedron -> Property +prop_b0_v3_always_zero t = + (volume t) > 0 ==> (b0 t) (v3 t) ~= 0 + +-- | The barycentric coordinate of v1 with respect to itself should +-- be one. +prop_b1_v1_always_unity :: Tetrahedron -> Property +prop_b1_v1_always_unity t = + (volume t) > 0 ==> (b1 t) (v1 t) ~= 1.0 + +-- | The barycentric coordinate of v0 with respect to v1 should +-- be zero. +prop_b1_v0_always_zero :: Tetrahedron -> Property +prop_b1_v0_always_zero t = + (volume t) > 0 ==> (b1 t) (v0 t) ~= 0 + +-- | The barycentric coordinate of v2 with respect to v1 should +-- be zero. +prop_b1_v2_always_zero :: Tetrahedron -> Property +prop_b1_v2_always_zero t = + (volume t) > 0 ==> (b1 t) (v2 t) ~= 0 + +-- | The barycentric coordinate of v3 with respect to v1 should +-- be zero. +prop_b1_v3_always_zero :: Tetrahedron -> Property +prop_b1_v3_always_zero t = + (volume t) > 0 ==> (b1 t) (v3 t) ~= 0 + +-- | The barycentric coordinate of v2 with respect to itself should +-- be one. +prop_b2_v2_always_unity :: Tetrahedron -> Property +prop_b2_v2_always_unity t = + (volume t) > 0 ==> (b2 t) (v2 t) ~= 1.0 + +-- | The barycentric coordinate of v0 with respect to v2 should +-- be zero. +prop_b2_v0_always_zero :: Tetrahedron -> Property +prop_b2_v0_always_zero t = + (volume t) > 0 ==> (b2 t) (v0 t) ~= 0 + +-- | The barycentric coordinate of v1 with respect to v2 should +-- be zero. +prop_b2_v1_always_zero :: Tetrahedron -> Property +prop_b2_v1_always_zero t = + (volume t) > 0 ==> (b2 t) (v1 t) ~= 0 + +-- | The barycentric coordinate of v3 with respect to v2 should +-- be zero. +prop_b2_v3_always_zero :: Tetrahedron -> Property +prop_b2_v3_always_zero t = + (volume t) > 0 ==> (b2 t) (v3 t) ~= 0 + +-- | The barycentric coordinate of v3 with respect to itself should +-- be one. +prop_b3_v3_always_unity :: Tetrahedron -> Property +prop_b3_v3_always_unity t = + (volume t) > 0 ==> (b3 t) (v3 t) ~= 1.0 + +-- | The barycentric coordinate of v0 with respect to v3 should +-- be zero. +prop_b3_v0_always_zero :: Tetrahedron -> Property +prop_b3_v0_always_zero t = + (volume t) > 0 ==> (b3 t) (v0 t) ~= 0 + +-- | The barycentric coordinate of v1 with respect to v3 should +-- be zero. +prop_b3_v1_always_zero :: Tetrahedron -> Property +prop_b3_v1_always_zero t = + (volume t) > 0 ==> (b3 t) (v1 t) ~= 0 + +-- | The barycentric coordinate of v2 with respect to v3 should +-- be zero. +prop_b3_v2_always_zero :: Tetrahedron -> Property +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 + where + t' = t { v0 = (v1 t), v1 = (v0 t) } + +prop_swapping_vertices_doesnt_affect_coefficients2 :: Tetrahedron -> Bool +prop_swapping_vertices_doesnt_affect_coefficients2 t = + c t 0 1 1 1 == c t' 0 1 1 1 + where + t' = t { v2 = (v3 t), v3 = (v2 t) } + +prop_swapping_vertices_doesnt_affect_coefficients3 :: Tetrahedron -> Bool +prop_swapping_vertices_doesnt_affect_coefficients3 t = + c t 2 1 0 0 == c t' 2 1 0 0 + where + t' = t { v2 = (v3 t), v3 = (v2 t) } + +prop_swapping_vertices_doesnt_affect_coefficients4 :: Tetrahedron -> Bool +prop_swapping_vertices_doesnt_affect_coefficients4 t = + c t 2 0 0 1 == c t' 2 0 0 1 + where + t' = t { v0 = (v3 t), v3 = (v0 t) } + + + + +tetrahedron_tests :: Test.Framework.Test +tetrahedron_tests = + testGroup "Tetrahedron Tests" [ + tetrahedron1_geometry_tests, + tetrahedron2_geometry_tests, + containment_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] + + +tetrahedron_properties :: Test.Framework.Test +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, + testProperty "b0_v2_always_zero" prop_b0_v2_always_zero, + testProperty "b0_v3_always_zero" prop_b0_v3_always_zero, + testProperty "b1_v1_always_unity" prop_b1_v1_always_unity, + testProperty "b1_v0_always_zero" prop_b1_v0_always_zero, + testProperty "b1_v2_always_zero" prop_b1_v2_always_zero, + testProperty "b1_v3_always_zero" prop_b1_v3_always_zero, + testProperty "b2_v2_always_unity" prop_b2_v2_always_unity, + testProperty "b2_v0_always_zero" prop_b2_v0_always_zero, + testProperty "b2_v1_always_zero" prop_b2_v1_always_zero, + testProperty "b2_v3_always_zero" prop_b2_v3_always_zero, + testProperty "b3_v3_always_unity" prop_b3_v3_always_unity, + testProperty "b3_v0_always_zero" prop_b3_v0_always_zero, + testProperty "b3_v1_always_zero" prop_b3_v1_always_zero, + testProperty "b3_v2_always_zero" prop_b3_v2_always_zero, + testProperty "swapping_vertices_doesnt_affect_coefficients1" $ + prop_swapping_vertices_doesnt_affect_coefficients1, + testProperty "swapping_vertices_doesnt_affect_coefficients2" $ + prop_swapping_vertices_doesnt_affect_coefficients2, + testProperty "swapping_vertices_doesnt_affect_coefficients3" $ + prop_swapping_vertices_doesnt_affect_coefficients3, + testProperty "swapping_vertices_doesnt_affect_coefficients4" $ + prop_swapping_vertices_doesnt_affect_coefficients4 ]