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 qualified Algebra.ToRational as ToRational ( C )
24 import Data.Vector.Fixed ( Arity, S )
26 import qualified Prelude as P
28 import Integration.Gaussian ( gaussian )
29 import Linear.Matrix (
36 import Polynomials.Orthogonal ( legendre )
38 -- | Dirichlet boundary conditions. Since u(a)=u(b)=0 are fixed,
39 -- there's no additional information conveyed by this type.
40 data Dirichlet a = Dirichlet { domain_dirichlet :: Interval a }
42 -- | Neumann boundary conditions. @alpha@ specifies A(a)u'(b) and
43 -- @beta@ specifies A(b)u'(b).
45 Neumann { domain_neumann :: Interval a,
49 -- | Boundary conditions can be either Dirichlet or Neumann.
50 type BoundaryConditions a = Either (Dirichlet a) (Neumann a)
52 type Interval a = (a,a)
54 -- | The data for the PDE that we are attempting to solve.
64 -- | The boundary conditions. The domain also specifies the
66 bdy :: BoundaryConditions a }
70 -- | Non-PDE parameters for the finite element method. The additional
71 -- type parameter @n@ should be a type-level representation of the
72 -- largest element in @max_degrees@. It needs to be known statically
73 -- for the dimensions of the pointer matrix. The parameter @l@ is
74 -- the number of global basis functions. It's equal to the number of
75 -- /internal/ mesh nodes (i.e. m-1), plus the sum of (p_i - 1) for
76 -- each p_i in max_degrees.
79 -- | A partition of the domain.
80 mesh :: Col m (Interval a),
82 -- | A list of the highest-degree polynomial that we will use over
83 -- each interval in our mesh.
84 max_degrees :: Col m Int }
87 -- | The pointer matrix mapping local basis elements to global
88 -- ones. The (i,j)th entry of the matrix contains the index of the
89 -- global basis function corresponding to the jth local basis
90 -- function over element i.
92 -- TODO: This works for Dirichlet boundary conditions only.
94 -- Note that the number of columns in the matrix depends on the
98 -- >>> import Data.Vector.Fixed ( N5, N6 )
99 -- >>> import Linear.Matrix ( Col5, fromList )
100 -- >>> import Naturals ( N19 )
102 -- >>> let p = fromList [[3],[3],[5],[4],[5]] :: Col5 Int
103 -- >>> let mesh = undefined :: Col5 (Int,Int)
104 -- >>> let params = Params mesh p :: Params N5 N5 N19 Int
105 -- >>> let row1 = [0,1,5,6,0,0] :: [Int]
106 -- >>> let row2 = [1,2,7,8,0,0] :: [Int]
107 -- >>> let row3 = [2,3,9,10,11,12] :: [Int]
108 -- >>> let row4 = [3,4,13,14,15,0] :: [Int]
109 -- >>> let row5 = [4,0,16,17,18,19] :: [Int]
110 -- >>> let expected = fromList [row1,row2,row3,row4,row5] :: Mat N5 N6 Int
111 -- >>> pointer params == expected
114 pointer :: (Arity m, Arity n, Arity l) => Params m n l a -> Mat m (S n) Int
118 -- | The number of polynomials in the local basis for element i.
119 numpolys :: Int -> Int
120 numpolys i = ((max_degrees params) !!! (i,0)) + 1
124 -- The first two (linear) basis functions are easy to do.
125 | j < 2 = if i + j >= (nrows $ mesh params) then 0 else i + j
127 -- No local basis function exists for this j.
128 | j >= (numpolys i) = 0
130 -- The first row we handle as a special case.
131 | i == 0 = (nrows $ mesh params) + j - 2
133 -- The first entry for this row not corresponding to one of the
134 -- linear polynomials (which we did in the first guard). This
135 -- grabs the biggest number in the previous row and begins counting
137 | j == 2 = (lambda (i-1) (numpolys (i-1) - 1)) + 1
139 -- If none of the other cases apply, we can compute our value by
140 -- looking to the left in the matrix.
141 | otherwise = (lambda i (j-1)) + 1
145 -- | An affine transform taking the interval [x1,x2] into [-1,1].
149 -- >>> let phi = affine (-6,7)
155 affine :: Field.C a => (a,a) -> (a -> a)
156 affine (x1,x2) x = (fromInteger 2)*(x - x1)/(x2 - x1) - (fromInteger 1)
158 -- | The inverse of 'affine'. It should send [-1,1] into [x1,x2].
162 -- >>> let phi = affine_inv (-6,7)
168 affine_inv :: Field.C a => (a,a) -> (a -> a)
169 affine_inv (x1,x2) x =
170 x*(x2 - x1)/two + (x1 + x2)/two
175 -- | Orthonormal basis functions over [-1,1]. The test case below
176 -- comes from Sage where the orthogonality of the polynomials'
177 -- derivatives can easily be tested.
181 -- >>> import qualified Algebra.Absolute as Absolute ( abs )
183 -- >>> let expected = 2.99624907925257
184 -- >>> let actual = big_N 4 1.5 :: Double
185 -- >>> Absolute.abs (actual - expected) < 1e-12
188 big_N :: forall a. (Algebraic.C a, RealField.C a) => Integer -> a -> a
190 | k < 0 = error "requested a negative basis function"
191 | k == 0 = (one - x) / (fromInteger 2)
192 | k == 1 = (one + x) / (fromInteger 2)
194 coeff * ( legendre k x - legendre (k-2) x )
198 coeff = one / (sqrt (four*(fromInteger k) - two)) :: a
201 -- | A matrix containing 'big_N' functions indexed by their
202 -- element/number. Each row in the matrix represents a finite element
203 -- (i.e. an interval in the mesh). Within row @i@, column @j@ contains
204 -- the @j@th 'big_N' basis function.
206 -- Any given 'big_N' will probably wind up in this matrix multiple
207 -- times; the basis functions don't change depending on the
208 -- interval. Only the /number/ of basis functions does. Computing
209 -- them this way allows us to easily construct a lookup \"table\" of
210 -- the proper dimensions.
212 -- The second example below relies on the fact that @big_N 3@ and
213 -- @big_N 6@ expand to Legendre polynomials (2,4) and (5,7)
214 -- respectively and so should be orthogonal over [-1,1].
218 -- >>> import Data.Vector.Fixed ( N5 )
219 -- >>> import Integration.Gaussian ( gaussian )
220 -- >>> import Linear.Matrix ( Col5, fromList )
221 -- >>> import Naturals ( N19 )
223 -- >>> let p = fromList [[3],[3],[5],[4],[5]] :: Col5 Int
224 -- >>> let mesh = undefined :: Col5 (Double,Double)
225 -- >>> let params = Params mesh p :: Params N5 N5 N19 Double
226 -- >>> let big_ns = big_N_matrix params
227 -- >>> let n1 = big_ns !!! (1,0)
228 -- >>> let n4 = big_ns !!! (4,0)
229 -- >>> n1 1.5 == n4 1.5
231 -- >>> let n1 = big_ns !!! (1,3)
232 -- >>> let n2 = big_ns !!! (2,4)
233 -- >>> gaussian (\x -> (n1 x) * (n2 x)) < 1e-12
236 big_N_matrix :: (Arity m, Arity n, Arity l, Algebraic.C a, RealField.C a)
238 -> Mat m (S n) (a -> a)
239 big_N_matrix params =
242 maxdeg i = (max_degrees params) !!! (i,0)
243 lambda i j = if j > maxdeg i
244 then (const $ fromInteger 0)
245 else big_N (toInteger j)
249 -- | Compute the global load vector F.
253 -- >>> import Data.Vector.Fixed ( N3, N4 )
254 -- >>> import Linear.Matrix ( Col4, frobenius_norm, fromList )
255 -- >>> import Naturals ( N7 )
257 -- >>> let big_A = const (1::Double)
258 -- >>> let c x = sin x
259 -- >>> let f x = x*(sin x)
260 -- >>> let bdy = Left (Dirichlet (0,1::Double))
261 -- >>> let pde = PDE big_A c f bdy
263 -- >>> let i1 = (0.0,1/3)
264 -- >>> let i2 = (1/3,2/3)
265 -- >>> let i3 = (2/3,4/5)
266 -- >>> let i4 = (4/5,1.0)
267 -- >>> let mesh = fromList [[i1], [i2], [i3], [i4]] :: Col4 (Double,Double)
268 -- >>> let pvec = fromList [[2],[3],[2],[1]] :: Col4 Int
269 -- >>> let params = Params mesh pvec :: Params N4 N3 N7 Double
271 -- >>> let f1 = [0.0418]
272 -- >>> let f2 = [0.0805]
273 -- >>> let f3 = [0.1007]
274 -- >>> let f4 = [-0.0045]
275 -- >>> let f5 = [-0.0332]
276 -- >>> let f6 = [-0.0054]
277 -- >>> let f7 = [-0.0267]
278 -- >>> let expected = fromList [f1,f2,f3,f4,f5,f6,f7] :: Col N7 Double
279 -- >>> let actual = big_F pde params
280 -- >>> frobenius_norm (actual - expected) < 1e-4
283 big_F :: forall m n l a.
284 (Arity l, Arity m, Arity n,
285 Algebraic.C a, RealField.C a, ToRational.C a)
290 ifoldl2 accum zero (big_N_matrix params)
292 accum :: Int -> Int -> Col l a -> (a -> a) -> Col l a
293 accum i j prev_F this_N =
297 (x1,x2) = (mesh params) !!! (i,0)
298 q = affine_inv (x1,x2)
299 integrand x = ((f pde) (q x)) * (this_N x)
301 -- The pointer matrix numbers from 1 so subtract one here to
302 -- get the right index.
303 global_idx = ((pointer params) !!! (i,j)) - 1
304 this_F = construct lambda
305 lambda k _ = if k == global_idx
306 then (gaussian integrand)*(x2 - x1) / two