X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FTetrahedron.hs;h=5b459fea2870de5c2fed88ee9f9e5b57c7e5cfc9;hb=5a21d504c18c0fa28d17286c0c5c9e7929d7f59e;hp=81fe91c11772122d0643b884e4fc3aea52163a43;hpb=6fb9ab6b6068870323e996da931fc04c7710e3e4;p=spline3.git diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs index 81fe91c..5b459fe 100644 --- a/src/Tetrahedron.hs +++ b/src/Tetrahedron.hs @@ -3,23 +3,51 @@ where import Numeric.LinearAlgebra hiding (i, scale) import Prelude hiding (LT) +import Test.QuickCheck (Arbitrary(..), Gen, choose) import Cardinal +import Comparisons (nearly_ge) import FunctionValues 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 { fv :: FunctionValues, + v0 :: Point, + v1 :: Point, + v2 :: Point, + v3 :: Point, + precomputed_volume :: Double, + + -- | Between 0 and 23; used to quickly determine which + -- tetrahedron I am in the parent 'Cube' without + -- having to compare them all. + number :: Int + } + 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 + rnd_no <- choose (0,23) + + -- 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 rnd_no + let vol = volume t' + return (Tetrahedron rnd_fv rnd_v0 rnd_v1 rnd_v2 rnd_v3 vol rnd_no) + instance Show Tetrahedron where show t = "Tetrahedron:\n" ++ + " no: " ++ (show (number t)) ++ "\n" ++ " fv: " ++ (show (fv t)) ++ "\n" ++ " v0: " ++ (show (v0 t)) ++ "\n" ++ " v1: " ++ (show (v1 t)) ++ "\n" ++ @@ -28,9 +56,32 @@ instance Show Tetrahedron where instance ThreeDimensional Tetrahedron where - center t = ((v0 t) + (v1 t) + (v2 t) + (v3 t)) `scale` (1/4) + center (Tetrahedron _ v0' v1' v2' v3' _ _) = + (v0' + v1' + v2' + v3') `scale` (1/4) + contains_point t p = - (b0 t p) >= 0 && (b1 t p) >= 0 && (b2 t p) >= 0 && (b3 t p) >= 0 + 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 = p } + + b1_unscaled :: Double + b1_unscaled = volume inner_tetrahedron + where inner_tetrahedron = t { v1 = p } + + b2_unscaled :: Double + b2_unscaled = volume inner_tetrahedron + where inner_tetrahedron = t { v2 = p } + + b3_unscaled :: Double + b3_unscaled = volume inner_tetrahedron + where inner_tetrahedron = t { v3 = p } polynomial :: Tetrahedron -> (RealFunction Point) @@ -76,6 +127,11 @@ 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) @@ -197,28 +253,22 @@ c _ _ _ _ _ = error "coefficient index out of bounds" +-- | The matrix used in the tetrahedron volume calculation as given in +-- Lai & Schumaker, Definition 15.4, page 436. vol_matrix :: Tetrahedron -> Matrix Double -vol_matrix t = (4><4) $ +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. + (x1, y1, z1) = v0 t + (x2, y2, z2) = v1 t + (x3, y3, z3) = v2 t + (x4, y4, z4) = v3 t + +-- | Computed using the formula from Lai & Schumaker, Definition 15.4, +-- page 436. volume :: Tetrahedron -> Double volume t | (v0 t) == (v1 t) = 0 @@ -230,22 +280,29 @@ volume t | otherwise = (1/6)*(det (vol_matrix 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 }