Linear Finite Elements

Beginning with this post, I’ll be publishing about the basics of finite element formulations, from personal notes that accumulated over the years. This one is about linear and scalar problems which came to be the “Hello World” for FE. Details regarding spaces and discretization are omitted for the sake of brevity. For those who want to delve into theory, I recommend “The Finite Element Method: Theory, Implementation, and Applications” by Larson and Bengzon.

The weak formulation of a canonical linear problem reads

Find $u\in V$ such that

$$\begin{equation} a(u, v) = b(v) \htmlId{eq:femlinear1}{} \tag{1}\end{equation}$$

for all $v \in V$ where $a(\cdot, \cdot)$ is a bilinear form and $b(\cdot)$ is a linear form.

We define the discretization of $u$ as

$$\begin{equation} u_h := \suml{J=1}{\nnode} u^J N^J ,\quad u_h \in V_h \quad\text{where}\quad V_h\subset V \tag{2}\end{equation}$$

The discretization $u_h$ is a linear combination of basis functions $N^J$ and corresponding scalars $u^J$, $J=1,\dots,\nnode$ so that $V_h$ is a subset of $V$. The discretization of (1) then reads

$$\begin{equation} a(u_h, v_h) = b(v_h) \quad \forall v_h \in V_h . \tag{3}\end{equation}$$

We then have

$$\begin{equation} a\rbr{\suml{J=1}{\nnode} u^J N^J, \suml{I=1}{\nnode} v^I N^I} = b\rbr{\suml{I=1}{\nnode} v^I N^I} \tag{4}\end{equation}$$

Using the linearity properties,

$$\begin{equation} a(\alpha u, \beta v) = \alpha\beta\, a(u,v) \eqand b(\alpha v) = \alpha b(v) \tag{5}\end{equation}$$

we obtain

$$\begin{equation} \suml{I=1}{\nnode} \suml{J=1}{\nnode} u^J v^I a(N^J, N^I) = \suml{I=1}{\nnode} v^I b(N^I) . \htmlId{eq:femlinear2}{} \tag{6}\end{equation}$$

For arbitrary test function values $v^I$, we can express (6) as a system of $\nnode$ equations

$$\begin{equation} \suml{J=1}{\nnode} u^J a(N^J, N^I) = b(N^I) \htmlId{eq:femlinear3}{} \tag{7}\end{equation}$$

for $I = 1,2,\dots,\nnode$. If we expand the summations as

$$\begin{alignat*}{6} & a(N^1, N^1) u^1 &&+ a(N^2, N^1) u^2 &&+ \cdots &&+ a(N^{\nnode}, N^1) u^{\nnode} &&\quad=\quad b(N^1) \\ & a(N^1, N^2) u^1 &&+ a(N^2, N^2) u^2 &&+ \cdots &&+ a(N^{\nnode}, N^2) u^{\nnode} &&\quad=\quad b(N^2) \\ & \qquad\vdots && \qquad\quad\;\vdots && \quad\;\;\vdots && \qquad\qquad\vdots && \qquad\qquad\vdots \\ & a(N^1, N^{\nnode}) u^1 &&+ a(N^2, N^{\nnode}) u^2 &&+ \cdots &&+ a(N^{\nnode}, N^{\nnode}) u^{\nnode} &&\quad=\quad b(N^{\nnode}) \end{alignat*}$$

we can see that the terms with $a$ constitute a matrix $\BA$ and the terms with $b$ constitute a vector $\Bb$, allowing us to write

$$\begin{equation} \BA\Bu = \Bb \htmlId{eq:discrete9}{} \tag{8}\end{equation}$$

where we chose to express the unknown coefficients $u^I$ as a vector $\Bu = [u^1,u^2,\dots,u^{\nnode}]\tra$.

\It can be seen that the components of the $\BA$ and $\Bb$ are defined as

$$\begin{equation} \boxed{ \Aelid{I\!J}{} = a(N^J,N^I) \eqand b^I = b(N^I), } \tag{9}\end{equation}$$

we can express the linear system as

$$\begin{alignat*}{6} & \Aelid{11}{} u^1 &&+ \Aelid{12}{} u^2 &&+ \cdots &&+ \Aelid{1\nnode}{} u^{\nnode} &&\quad=\quad b^1 \\ & \Aelid{21}{} u^1 &&+ \Aelid{22}{} u^2 &&+ \cdots &&+ \Aelid{2\nnode}{} u^{\nnode} &&\quad=\quad b^2 \\ & \quad\vdots && \qquad\;\vdots && \quad\;\;\vdots && \qquad\;\vdots && \qquad\quad\;\;\vdots \\ & \Aelid{\nnode 1}{} u^1 &&+ \Aelid{\nnode 2}{} u^2 &&+ \cdots &&+ \Aelid{\nnode\nnode}{} u^{\nnode} &&\quad=\quad b^{\nnode} \end{alignat*}$$

Note that with the given definitions, (7) becomes

$$\begin{equation} \boxed{ \suml{J=1}{\nnode} \Aelid{I\!J}{} \,u^J = b^I \quad\text{for}\quad I=1,2,\dots\nnode. } \htmlId{eq:discrete10}{} \tag{10}\end{equation}$$

Example: Poisson’s Equation

In the weak form of Poisson’s equation

$$\begin{equation} \begin{alignedat}{4} - \Var u &= f \quad && \text{in} \quad && \Omega \\ u &= 0 \quad && \text{on} \quad && \del\Omega \end{alignedat} \tag{11}\end{equation}$$

The weak formulation reads

Find $u\in V$ such that

$$\begin{equation} - \int_\Omega \Delta(u) v \dv= \int_\Omega f v \dv \tag{12}\end{equation}$$

for all $v\in V$ where $V=H^1_0(\Omega)$.

Applying integration by parts and divergence theorem on the left-hand side

$$\begin{equation} \begin{aligned} \int_\Omega \Delta(u) v \dv &= \int_\Omega \nabla \dtp (\nabla (u) v) \dv - \int_\Omega \nabla u\dtp\nabla v \dv \\ &= \underbrace{\int_{\del\Omega} v (\nabla u\dtp\Bn) \da}_{v = 0 \text{ on } \del\Omega} - \int_\Omega \nabla u\dtp\nabla v \dv \\ \end{aligned} \tag{13}\end{equation}$$

We have the following variational forms:

$$\begin{equation} \begin{aligned} a(u,v) &= \int_{\Omega} \nabla u \dtp \nabla v \dv\\ b(v) &= \int_{\Omega} f \, v \dv\\ \end{aligned} \tag{14}\end{equation}$$

Following (7), we can calculate the stiffness matrix $\BA$ as

$$\begin{equation} \begin{aligned} \Aelid{I\!J}{} = a(N^J, N^I) &= \int_{\Omega} \nabla N^J \dtp \nabla N^I \dv \\ &= \int_{\Omega} \BB^J \dtp \BB^I \dv \end{aligned} \tag{15}\end{equation}$$

where we have defined the gradient of the basis functions as

$$\begin{equation} \BB^I := \nabla N^I\,. \tag{16}\end{equation}$$

Similarly, we integrate the force term into a vector $\Bb$ as

$$\begin{equation} \begin{aligned} b^I &= \int_{\Omega} f N^I \dv \end{aligned} \tag{17}\end{equation}$$