X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FTetrahedron.hs;h=9f68364042e7563a3b10d23f01e72f2845ba74af;hb=4d695b8d0b05a02d562bdee3d2a1b98ce6a5e747;hp=85ac4a372d222cea1eaa56a821f54cf340e08b82;hpb=f07f76b231a3df623aab8b6035ac6000ce2a5eb2;p=spline3.git diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs index 85ac4a3..9f68364 100644 --- a/src/Tetrahedron.hs +++ b/src/Tetrahedron.hs @@ -1,22 +1,48 @@ module Tetrahedron where +import qualified Data.Vector as V ( + singleton, + snoc, + sum + ) import Numeric.LinearAlgebra hiding (i, scale) import Prelude hiding (LT) +import Test.QuickCheck (Arbitrary(..), Gen) 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 + } + 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" ++ @@ -28,18 +54,56 @@ 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) 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,6 +140,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,6 +266,8 @@ 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) [1, 1, 1, 1, @@ -204,21 +275,13 @@ vol_matrix t = (4><4) 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 +293,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 }