Finite Element Method
Consider an ODE of the form:
with the following Boundary Conditions:
and
1.) Substitute for U an unspecified trial function into governing equation, i.e.
where
n = Number of nodes in system (or, locally, in element).
Uj = Unknown coefficients (to be determined), which are not a function of space, in general.
Nj = Known (user-specified) Basis function, which are a function of space, but not a function of time, in general.
2.) Force the Residual to zero in the global sense, i.e. satisfy the governing equation in the "weak form".
a.) Multiply by a Weighting function:
where Wi = Known (user-specified) Weight function.
b.) Integrate (take inner product) over global domain and force to zero.
Definition of an Inner Product:
Implies that W is orthogonal to R
< ( ) > =
(dimensionally dependent)
If you can solve this equation then you have solved the governing equations, at least in the weak form.
3.) Discretize global domain into subdomains (Elements) such that S (Elements) = contiguous, non-overlapping representation of the global domain.
4.) Select a basis function for - a few helpful hints in this selection are
a.) select an orthogonal series basis function
b.) choose a computationally efficient series
Ex. Lagrange Polynomials applied locally on each element
Assume Linear Element (N = 2)
Assume Quadratic Element (N = 3)
5.) Integration by parts:
Required for order PDEs using linear basis functions; Provides a direct means to apply Boundary Conditions
more generally,
Adding in the other components of the ODE yields:
6.) Select Weighting functions
Ex: Choose "Galerkin"
(However, for clarity retain the W symbol. It will help when assembling matrices later.)
Note that for each weighting function, Wi, one must assemble all contributions from each basis function Nj.
7.) Assemble Matrices – consider only 2 elements, and assemble all components associated with Row I (for Wi):
One row of equations for each Wi
Note: the local definitions of Trial and Weighting functions
for all
= 0 for all
for row i only contributions from
Left Elem Right Elem
= + 0 =
= + =
= 0 + =
= + 0 =
= + =
= 0 + =
= + =
Assembly of row I
If uniform spacing then h1 = h2 = h
+ Simpson’s Rule = g
Note the similarities to F.D.
In general, we are forming an expression
In a FE code: Loop over each element and assemble the local (elemental) contributions for the Lhs matrix [A] and the Rhs vector.
After this loop is completed, then (and only then) apply the boundary conditions to the set of equations.
FE codes contain a lot of bookkeeping. One common mapping array is the element connectivity (or Incidence) list. It has the form:
IN(K,L) where
L = Element Number, L = 1, NE ( and NE = Number of Elements)
K = Local Node Number, K = 1, 2, … to the number of nodes in an element.
IN(K,L) = global node number ( maps global node number to the local node number within element L)
Consider the following 1-D example:
The Element number can have significance if using a frontal matrix solver.
The Node numbering can have significance if using a banded matrix solver.
Node and Element numberings have less significance if using a sparse, iterative matrix solver.
The mapping for this example is:
IN(1,1) = 3 IN(2,1) = 1
IN(1,2) = 4 IN(2,2) = 6
IN(1,3) = 1 IN(2,3) = 4
IN(1,4) = 7 IN(2,4) = 2
IN(1,5) = 5 IN(2,5) = 7
IN(1,6) = 6 IN(2,6) = 5
The FE approximation to the governing equation is accumulated (summed) in 3 nested loops (L, I, and J)
Do L = 1, NE (Loop over all elements)
h = x(in(2,L)) – x(in(1,L)) (Calculate element length)
Do I = 1, 2 (Loop for all Weighting Functions)
dWdx = (1/h)
if (I = 1) dWdx = - dWdx
Do J = 1, 2 (Loop for all Trial Functions)
dNdx = (1/h)
if (J = 1) dNdx = - dNdx
term1 = - dNdx * dWdx * h
term2 = f * h/6
if (I = J) term2 = 2 * term2
AL(J, I) = term1 + term2 (Local 2x2 construction)
Irow = IN(I,L) (Global Mapping)
Jcol = IN(J,L)
AG(Jcol, Irow) += AL(J, I) (Global Matrix Construction)
End Do (End of J loop)
Rhs(Irow) += g * h/2 (Rhs Construction)
End Do (End of I loop)
End Do (Finished with all Elements)
Everything done at the element level!
Easy to automate
a) One Element (problem-specific)
b) Assembly (problem-independent)
8) Add in the Boundary Conditions
a.) Type 1 BC : Satisfy Exactly U(0) = 1
(Note this is at global node number 3 for our example.)
One less unknown in algebraic system.
Remove row_3 and then move column_3 of matrix to Rhs since it is known (not usually done, but possible)
This process is somewhat cumbersome. It can change the bandwidth or structure of your system. It can cause other solution-solving anomalies if one is not careful.
Cleaner process:
Remove row 3 completely, i.e. place zeros in all entries.
Then place a 1 on diagonal of [A]
Put the solution (i.e. U(0)=1) into the Rhs
b.) Type II B.C. Satisfy approximately in {BCs} vector.
Recall: and for node 2 (in our example) the flux at the boundary is 5. Also note:
At the boundary node Wi = 1 for I = 2 (global) and Wi = 0 for all other nodes.
Therefore to apply Type II Boundary Conditions in our example:
Rhs(2) = Rhs(2) – 5
The entire formulation is complete. Call a matrix solver and obtain the solutions for Uj.
ME 515 Finite Element Procedure (Skeleton) (Sullivan) 13