From: Michael Orlitzky Date: Thu, 6 Feb 2014 06:36:45 +0000 (-0500) Subject: Add the (in progress) FEM.R1 module. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=b9df47247151260ce8941562b53c040d779e8d2f;p=numerical-analysis.git Add the (in progress) FEM.R1 module. --- diff --git a/numerical-analysis.cabal b/numerical-analysis.cabal index 92b716e..b666583 100644 --- a/numerical-analysis.cabal +++ b/numerical-analysis.cabal @@ -18,6 +18,7 @@ data-files: makefile library exposed-modules: + FEM.R1 Integration.Gaussian, Integration.Simpson, Integration.Trapezoid, diff --git a/src/FEM/R1.hs b/src/FEM/R1.hs new file mode 100644 index 0000000..c6a76c8 --- /dev/null +++ b/src/FEM/R1.hs @@ -0,0 +1,169 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE ScopedTypeVariables #-} + +-- | Finite Element Method (FEM) to solve the problem, +-- +-- -(A(x)*u'(x))' + c(x)*u*(x) = f(x) on (a,b) +-- +-- with one of the two types of boundary conditions, +-- +-- * Dirichlet: u(a) = u(b) = 0 +-- +-- * A(a)u'(a) = alpha, A(b)u'(b) = beta +-- +-- if c(x) = 0, then it is assumed that the boundary conditions are +-- Dirichlet. +-- +module FEM.R1 +where + +import qualified Algebra.Algebraic as Algebraic ( C ) +import qualified Algebra.Field as Field ( C ) +import qualified Algebra.RealField as RealField ( C ) +import Data.Vector.Fixed ( Arity, S ) +import NumericPrelude +import qualified Prelude as P + +import Linear.Matrix ( + Col, + Mat(..), + (!!!), + construct, + nrows ) +import Polynomials.Orthogonal ( legendre ) + +-- | Dirichlet boundary conditions. Since u(a)=u(b)=0 are fixed, +-- there's no additional information conveyed by this type. +data Dirichlet = Dirichlet + +-- | Neumann boundary conditions. @alpha@ specifies A(a)u'(b) and +-- @beta@ specifies A(b)u'(b). +data Neumann a = Neumann { alpha :: a, beta :: a } + +-- | Boundary conditions can be either Dirichlet or Neumann. +type BoundaryConditions a = Either Dirichlet (Neumann a) + +type Interval a = (a,a) + +-- | The data for the PDE that we are attempting to solve. +data PDE a = + PDE { + -- | A(x) + big_A :: (a -> a), + -- | c(x) + c :: (a -> a), + -- | f(x) + f :: (a -> a), + -- | The domain in R^1 as an interval + domain :: Interval a, + bdy :: BoundaryConditions a } + + + +-- | Non-PDE parameters for the finite element method. The additional +-- type parameter @n@ should be a type-level representation of the +-- largest element in @max_degrees@. It needs to be known statically +-- for the dimensions of the pointer matrix. +data Params m n a = + Params { + -- | A partition of the domain. + mesh :: Col m (Interval a), + + -- | A list of the highest-degree polynomial that we will use over + -- each interval in our mesh. + max_degrees :: Col m Int } + + +-- | The pointer matrix mapping local basis elements to global +-- ones. The (i,j)th entry of the matrix contains the index of the +-- global basis function corresponding to the jth local basis +-- function over element i. +-- +-- TODO: This works for Dirichlet boundary conditions only. +-- +-- Note that the number of columns in the matrix depends on the +-- +-- Examples: +-- +-- >>> import Data.Vector.Fixed ( N5, N6 ) +-- >>> import Linear.Matrix ( Col5, fromList ) +-- +-- >>> let p = fromList [[3],[3],[5],[4],[5]] :: Col5 Int +-- >>> let mesh = undefined :: Col5 (Int,Int) +-- >>> let params = Params mesh p :: Params N5 N5 Int +-- >>> let row1 = [0,1,5,6,0,0] :: [Int] +-- >>> let row2 = [1,2,7,8,0,0] :: [Int] +-- >>> let row3 = [2,3,9,10,11,12] :: [Int] +-- >>> let row4 = [3,4,13,14,15,0] :: [Int] +-- >>> let row5 = [4,0,16,17,18,19] :: [Int] +-- >>> let expected = fromList [row1,row2,row3,row4,row5] :: Mat N5 N6 Int +-- >>> pointer params == expected +-- True +-- +pointer :: (Arity m, Arity n) => Params m n a -> Mat m (S n) Int +pointer params = + construct lambda + where + -- | The number of polynomials in the local basis for element i. + numpolys :: Int -> Int + numpolys i = ((max_degrees params) !!! (i,0)) + 1 + + lambda i j + + -- The first two (linear) basis functions are easy to do. + | j < 2 = if i + j >= (nrows $ mesh params) then 0 else i + j + + -- No local basis function exists for this j. + | j >= (numpolys i) = 0 + + -- The first row we handle as a special case. + | i == 0 = (nrows $ mesh params) + j - 2 + + -- The first entry for this row not corresponding to one of the + -- linear polynomials (which we did in the first guard). This + -- grabs the biggest number in the previous row and begins counting + -- from there. + | j == 2 = (lambda (i-1) (numpolys (i-1) - 1)) + 1 + + -- If none of the other cases apply, we can compute our value by + -- looking to the left in the matrix. + | otherwise = (lambda i (j-1)) + 1 + + + +-- | An affine transform taking the interval [x1,x2] into [-1,1]. +-- +-- Examples: +-- +-- >>> let phi = affine (-6,9) +-- >>> phi (-6) +-- -1.0 +-- >>> phi (9) +-- 1.0 +-- +affine :: Field.C a => (a,a) -> (a -> a) +affine (x1,x2) x = (fromInteger 2)*(x - x1)/(x2 - x1) - (fromInteger 1) + + +-- | Orthonormal basis functions over [-1,1]. The test case below +-- comes from Sage where the orthogonality of the polynomials' +-- derivatives can easily be tested. +-- +-- Examples: +-- +-- >>> import qualified Algebra.Absolute as Absolute ( abs ) +-- +-- >>> let expected = 6.33910180790284 +-- >>> let actual = big_N 3 1.5 :: Double +-- >>> Absolute.abs (actual - expected) < 1e-12 +-- True +-- +big_N :: forall a. (Algebraic.C a, RealField.C a) => Integer -> a -> a +big_N k x + | k < 0 = error "requested a negative basis function" + | otherwise = + coeff * ( legendre (k+2) x - legendre k x ) + where + four = fromInteger 4 + six = fromInteger 6 + coeff = one / (sqrt (four*(fromInteger k) + six)) :: a