2 {-# LANGUAGE ScopedTypeVariables #-}
4 -- | Finite Element Method (FEM) to solve the problem,
6 -- -(A(x)*u'(x))' + c(x)*u*(x) = f(x) on (a,b)
8 -- with one of the two types of boundary conditions,
10 -- * Dirichlet: u(a) = u(b) = 0
12 -- * A(a)u'(a) = alpha, A(b)u'(b) = beta
14 -- if c(x) = 0, then it is assumed that the boundary conditions are
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 )
25 import qualified Prelude as P
27 import Linear.Matrix (
33 import Polynomials.Orthogonal ( legendre )
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
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 }
43 -- | Boundary conditions can be either Dirichlet or Neumann.
44 type BoundaryConditions a = Either Dirichlet (Neumann a)
46 type Interval a = (a,a)
48 -- | The data for the PDE that we are attempting to solve.
57 -- | The domain in R^1 as an interval
59 bdy :: BoundaryConditions a }
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.
69 -- | A partition of the domain.
70 mesh :: Col m (Interval a),
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 }
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.
82 -- TODO: This works for Dirichlet boundary conditions only.
84 -- Note that the number of columns in the matrix depends on the
88 -- >>> import Data.Vector.Fixed ( N5, N6 )
89 -- >>> import Linear.Matrix ( Col5, fromList )
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
103 pointer :: (Arity m, Arity n) => Params m n a -> Mat m (S n) Int
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
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
116 -- No local basis function exists for this j.
117 | j >= (numpolys i) = 0
119 -- The first row we handle as a special case.
120 | i == 0 = (nrows $ mesh params) + j - 2
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
126 | j == 2 = (lambda (i-1) (numpolys (i-1) - 1)) + 1
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
134 -- | An affine transform taking the interval [x1,x2] into [-1,1].
138 -- >>> let phi = affine (-6,9)
144 affine :: Field.C a => (a,a) -> (a -> a)
145 affine (x1,x2) x = (fromInteger 2)*(x - x1)/(x2 - x1) - (fromInteger 1)
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.
154 -- >>> import qualified Algebra.Absolute as Absolute ( abs )
156 -- >>> let expected = 6.33910180790284
157 -- >>> let actual = big_N 3 1.5 :: Double
158 -- >>> Absolute.abs (actual - expected) < 1e-12
161 big_N :: forall a. (Algebraic.C a, RealField.C a) => Integer -> a -> a
163 | k < 0 = error "requested a negative basis function"
165 coeff * ( legendre (k+2) x - legendre k x )
169 coeff = one / (sqrt (four*(fromInteger k) + six)) :: a