X-Git-Url: https://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FTetrahedron.hs;h=87332b94440afdb999b173c13f681095d48ac796;hb=1f45b6c0880c08809163d18f5f7b63d773ddcbe7;hp=6c13f4bbb429f2a185cdd1eaae3a541f1b147cb8;hpb=b9c6d5f4f9cf3e0e8c32499959c879e6300717f6;p=spline3.git diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs index 6c13f4b..87332b9 100644 --- a/src/Tetrahedron.hs +++ b/src/Tetrahedron.hs @@ -1,9 +1,14 @@ 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 Test.QuickCheck (Arbitrary(..), Gen, choose) import Cardinal import Comparisons (nearly_ge) @@ -13,12 +18,20 @@ 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 @@ -28,11 +41,18 @@ instance Arbitrary Tetrahedron where rnd_v2 <- arbitrary :: Gen Point rnd_v3 <- arbitrary :: Gen Point rnd_fv <- arbitrary :: Gen FunctionValues - return (Tetrahedron rnd_fv rnd_v0 rnd_v1 rnd_v2 rnd_v3) + 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" ++ @@ -41,21 +61,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) `nearly_ge` 0 && - (b1 t p) `nearly_ge` 0 && - (b2 t p) `nearly_ge` 0 && - (b3 t p) `nearly_ge` 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. @@ -247,27 +302,27 @@ volume 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 }