src/FEM/R1.hs

   1 {-# LANGUAGE GADTs #-}
   2 {-# LANGUAGE ScopedTypeVariables #-}
   3
   4 -- | Finite Element Method (FEM) to solve the problem,
   5 --
   6 --     -(A(x)*u'(x))' + c(x)*u*(x) = f(x) on (a,b)
   7 --
   8 --   with one of the two types of boundary conditions,
   9 --
  10 --     * Dirichlet: u(a) = u(b) = 0
  11 --
  12 --     * A(a)u'(a) = alpha, A(b)u'(b) = beta
  13 --
  14 --   if c(x) = 0, then it is assumed that the boundary conditions are
  15 --   Dirichlet.
  16 --
  17 module FEM.R1
  18 where
  19
  20 import qualified Algebra.Algebraic as Algebraic ( C )
  21 import qualified Algebra.Field as Field ( C )
  22 import qualified Algebra.RealField as RealField ( C )
  23 import Data.Vector.Fixed ( Arity, S )
  24 import NumericPrelude
  25 import qualified Prelude as P
  26
  27 import Linear.Matrix (
  28   Col,
  29   Mat(..),
  30   (!!!),
  31   construct,
  32   nrows )
  33 import Polynomials.Orthogonal ( legendre )
  34
  35 -- | Dirichlet boundary conditions. Since u(a)=u(b)=0 are fixed,
  36 --   there's no additional information conveyed by this type.
  37 data Dirichlet = Dirichlet
  38
  39 -- | Neumann boundary conditions. @alpha@ specifies A(a)u'(b) and
  40 --   @beta@ specifies A(b)u'(b).
  41 data Neumann a = Neumann { alpha :: a, beta :: a }
  42
  43 -- | Boundary conditions can be either Dirichlet or Neumann.
  44 type BoundaryConditions a = Either Dirichlet (Neumann a)
  45
  46 type Interval a = (a,a)
  47
  48 -- | The data for the PDE that we are attempting to solve.
  49 data PDE a =
  50   PDE {
  51     -- | A(x)
  52     big_A   :: (a -> a),
  53     -- | c(x)
  54     c      :: (a -> a),
  55     -- | f(x)
  56     f      :: (a -> a),
  57     -- | The domain in R^1 as an interval
  58     domain :: Interval a,
  59     bdy    :: BoundaryConditions a }
  60
  61
  62
  63 -- | Non-PDE parameters for the finite element method. The additional
  64 --   type parameter @n@ should be a type-level representation of the
  65 --   largest element in @max_degrees@. It needs to be known statically
  66 --   for the dimensions of the pointer matrix.
  67 data Params m n a =
  68   Params {
  69     -- | A partition of the domain.
  70     mesh :: Col m (Interval a),
  71
  72     -- | A list of the highest-degree polynomial that we will use over
  73     --   each interval in our mesh.
  74     max_degrees :: Col m Int }
  75
  76
  77 -- | The pointer matrix mapping local basis elements to global
  78 --   ones. The (i,j)th entry of the matrix contains the index of the
  79 --   global basis function corresponding to the jth local basis
  80 --   function over element i.
  81 --
  82 --   TODO: This works for Dirichlet boundary conditions only.
  83 --
  84 --   Note that the number of columns in the matrix depends on the
  85 --
  86 --   Examples:
  87 --
  88 --   >>> import Data.Vector.Fixed ( N5, N6 )
  89 --   >>> import Linear.Matrix ( Col5, fromList )
  90 --
  91 --   >>> let p = fromList [[3],[3],[5],[4],[5]] :: Col5 Int
  92 --   >>> let mesh = undefined :: Col5 (Int,Int)
  93 --   >>> let params = Params mesh p :: Params N5 N5 Int
  94 --   >>> let row1 = [0,1,5,6,0,0] :: [Int]
  95 --   >>> let row2 = [1,2,7,8,0,0] :: [Int]
  96 --   >>> let row3 = [2,3,9,10,11,12] :: [Int]
  97 --   >>> let row4 = [3,4,13,14,15,0] :: [Int]
  98 --   >>> let row5 = [4,0,16,17,18,19] :: [Int]
  99 --   >>> let expected = fromList [row1,row2,row3,row4,row5] :: Mat N5 N6 Int
 100 --   >>> pointer params == expected
 101 --   True
 102 --
 103 pointer :: (Arity m, Arity n) => Params m n a -> Mat m (S n) Int
 104 pointer params =
 105   construct lambda
 106   where
 107     -- | The number of polynomials in the local basis for element i.
 108     numpolys :: Int -> Int
 109     numpolys i = ((max_degrees params) !!! (i,0)) + 1
 110
 111     lambda i j
 112
 113       -- The first two (linear) basis functions are easy to do.
 114       | j < 2 = if i + j >= (nrows $ mesh params) then 0 else i + j
 115
 116       -- No local basis function exists for this j.
 117       | j >= (numpolys i) = 0
 118
 119       -- The first row we handle as a special case.
 120       | i == 0 = (nrows $ mesh params) + j - 2
 121
 122       -- The first entry for this row not corresponding to one of the
 123       -- linear polynomials (which we did in the first guard). This
 124       -- grabs the biggest number in the previous row and begins counting
 125       -- from there.
 126       | j == 2 = (lambda (i-1) (numpolys (i-1) - 1)) + 1
 127
 128       -- If none of the other cases apply, we can compute our value by
 129       -- looking to the left in the matrix.
 130       | otherwise = (lambda i (j-1)) + 1
 131
 132
 133
 134 -- | An affine transform taking the interval [x1,x2] into [-1,1].
 135 --
 136 --   Examples:
 137 --
 138 --   >>> let phi = affine (-6,9)
 139 --   >>> phi (-6)
 140 --   -1.0
 141 --   >>>  phi (9)
 142 --   1.0
 143 --
 144 affine :: Field.C a => (a,a) -> (a -> a)
 145 affine (x1,x2) x = (fromInteger 2)*(x - x1)/(x2 - x1) - (fromInteger 1)
 146
 147
 148 -- | Orthonormal basis functions over [-1,1]. The test case below
 149 --   comes from Sage where the orthogonality of the polynomials'
 150 --   derivatives can easily be tested.
 151 --
 152 --   Examples:
 153 --
 154 --   >>> import qualified Algebra.Absolute as Absolute ( abs )
 155 --
 156 --   >>> let expected = 6.33910180790284
 157 --   >>> let actual = big_N 3 1.5 :: Double
 158 --   >>> Absolute.abs (actual - expected) < 1e-12
 159 --   True
 160 --
 161 big_N :: forall a. (Algebraic.C a, RealField.C a) => Integer -> a -> a
 162 big_N k x
 163   | k < 0 = error "requested a negative basis function"
 164   | otherwise =
 165       coeff * ( legendre (k+2) x - legendre k x )
 166       where
 167         four = fromInteger 4
 168         six = fromInteger 6
 169         coeff = one / (sqrt (four*(fromInteger k) + six)) :: a