What Are FEM

Finite Element Methods are methods for taking a Partial Differential Equation(PDE) from a continuous equation to a system of discrete equations. PDEs are one of the most important types of equations for modeling physical phenomena and FEM are only one scheme to solve PDEs numerically. FEM were developed by structural and mechanical engineers in the 1950's are in essence a variational method sometimes refered to as a method of weighted residuals which dates back to 19th century.

On this page, we overview a bit of the basic mathematics behind FEM. For this tutorial, we will be using the Poisson equation as our model problem:

Find $u$ on $[0,1]$ such that:

(1)

\begin{align} -\frac{d}{dx} u = f \end{align}

where $u(0) = 0, u'(1) = 0$

Variational Methods

The basic idea behind variational methods is exactly what is sounds like; one sets up the equations such there is some quantity that is able to vary according to some rules. Our model problem, Eq 1, is the strong form and does not lend to itself to a variational setting. To put the problem in a variational form, we integrate it against another function, let's say $v$.

(2)

\begin{eqnarray} (f,v) &=& \int_0^1 f(x)\cdot v(x) dx \\ &=& \int_0^1 -u'(x) \cdot v(x) dx \\ &=& -u(1)\cdot v(1) +u(0)\cdot v(0) + \int_0^1 u'(x)\cdot v'(x) dx \qquad \textrm{(Integration by parts) }\\ &=& \int_0^1 u'(x)\cdot v'(x) dx := a(u,v) \end{eqnarray}

By doing so the problem becomes:

For $\forall v \in V$ where $V = \{v | v \in L^2([0,1]) \textrm{ and } v(0)=0\}$ find $u$ such that:

(3)

\begin{equation} a(u,v) = (f,v) \end{equation}

In Eq 3 (often referred to as the weak or variational form), you can see that we now have a quantity (namely $v$) that we can allow to vary to solve the problem. As is the problem is well defined but not exactly computable, how are we going to find the exact $u$ if it must comply with all $v$? Is there some way to limit the number of functions $v$ we must comply with?

This sort of problem is quite common in scientific computing, we cannot really compute the exact answer but instead we focus on ways to solve a problem with nearly the same answer. For Eq 3 we rephrase the problem to only include a finite dimensional subspace $S$, thus we have:

For $\forall v_h \in S$ where $S\subset V$ find $u_h$ such that:

(4)

\begin{equation} a(u_h,v_h) = (f,v_h) \end{equation}

We use a finite dimensional subspace because then we can use a basis of that space and thus we only have to manage the number of equations as the dimension of the space. How do we determine how to pick a $S$ to approximate $V$?

The Finite Element

The finite element is the mathematical tool used to preserve certain features of the large space $V$ in a finite dimensional context. Some examples of features that we may wish to preserve are varying degrees of continuity and smoothness. Also we should try to keep the error within the element as small as possible up to a certain degree of polynomial.

So now that we know what the finite element is suppose to do what is it?

Any numerical simulation will partition the domain with a mesh. FEM are no different, and the mesh is what is used to define the elements. For each element we associate a number of degrees of freedom (dof), each dof will be associated with a discrete equation in our system. Each element will define a function space basis (functions that will span the space $S$), a set of nodes (or dof) to evaluate the basis on, and a geometry which is defined by the mesh.

For example take a line $[0,1]$ partitioned by $n$ points. The geometry of each element, $e$ will be segment $[x_i,x_{i+1}]$. If we are using linears (and thus we will be continuous across element boundaries by only exact up to a linear constant) each element will have two basis functions, with the dofs $x_i, x_{i+1}$

(5)

\begin{eqnarray} \phi^e_0 &=& \frac{x-x_{i+1}}{x_{i}-x_{i+1}} \\ \phi^e_1 &=& \frac{x-x_{i}}{x_{i+1}-x_i} \end{eqnarray}

Another example for our partitioned line would be quadratic functions. The geometry of the element stays the same, but for quadratic functions we need three points to uniquely define a quadratic function, thus the dof will be $x_i, x_{i+.5}, x_{i+1}$ with the basis functions:

(6)

\begin{eqnarray} \phi^e_0 &=& \frac{(x-x_{i+.5})(x-x_{i+1})}{(x_{i}-x_{i+.5})(x_{i}-x_{i+1})}\\ \phi^e_{.5} &=& \frac{(x-x_i)(x-x_{i+1})}{(x_{i+.5}-x_i)(x_{i+.5}-x_{i+1})} \\ \phi^e_1 &=& \frac{(x-x_i)(x-x_{i+.5})}{(x_{i+1}-x_i)(x_{i+1}-x_{i+.5})} \end{eqnarray}

How did I get these basis functions for that geometry and nodal basis?

First since each dof, we are trying to get its value in the global function we want it to have basis

(7)

\begin{align} \phi_i(x_j) = \delta_{ij} = \left\{\begin{array}{cl}1 & \textrm{if } i = j\\ 0 & \textrm{otherwise}\end{array}\right. \end{align}

But this is not a smooth function, so we approximate it by a $n$ order polynomial:

(8)

\begin{align} \phi_i = \frac{(x-x_0)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}{(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n)} \end{align}

The Matrix Equations

Now that we have an idea of how to make a the finite elements for our simulation, we need to put them together to form the matrix equation:

(9)

\begin{equation} K u = F \end{equation}

where

• $K$ is our stiffness matrix, such that $K[i,j] = a(\phi_i,\phi_j)$,
• $F$ is the load vector, such that $F[i] = (f,\phi_i)$, and
• $u$ is the solution vector.

Thus to form the matrix we only need to iterate over each element, compute the appropriate pieces for the stiffness matrix and load vector then add them to the matrix equation and solve the linear equation. Since this might be a quite a bit to understand, we are going to pause and code a bit. Luckily for you we will be using Python which has libraries defined to integrate and solve linear equations.