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∈Vu\in V such that

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

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

We define the discretization of uu as

uh:=∑J=1nnuJNJ,uh∈VhwhereVh⊂V\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 uhu_h is a linear combination of basis functions NJN^J and corresponding scalars uJu^J, J=1,…,nnJ=1,\dots,\nnode so that VhV_h is a subset of VV. The discretization of (1) then reads

a(uh,vh)=b(vh)∀vh∈Vh.\begin{equation} a(u_h, v_h) = b(v_h) \quad \forall v_h \in V_h . \tag{3}\end{equation}

We then have

a(∑J=1nnuJNJ,∑I=1nnvINI)=b(∑I=1nnvINI)\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,

a(αu,βv)=αβ a(u,v)andb(αv)=αb(v)\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

∑I=1nn∑J=1nnuJvIa(NJ,NI)=∑I=1nnvIb(NI).\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 vIv^I, we can express (6) as a system of nn\nnode equations

∑J=1nnuJa(NJ,NI)=b(NI)\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,…,nnI = 1,2,\dots,\nnode. If we expand the summations as

a(N1,N1)u1+a(N2,N1)u2+⋯+a(Nnn,N1)unn=b(N1)a(N1,N2)u1+a(N2,N2)u2+⋯+a(Nnn,N2)unn=b(N2)⋮  ⋮    ⋮⋮⋮a(N1,Nnn)u1+a(N2,Nnn)u2+⋯+a(Nnn,Nnn)unn=b(Nnn)\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 aa constitute a matrix A\BA and the terms with bb constitute a vector b\Bb, allowing us to write

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

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

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

AI ⁣J=a(NJ,NI)andbI=b(NI),\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

A11u1+A12u2+⋯+A1nnunn=b1A21u1+A22u2+⋯+A2nnunn=b2⋮  ⋮    ⋮  ⋮    ⋮Ann1u1+Ann2u2+⋯+Annnnunn=bnn\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

∑J=1nnAI ⁣J uJ=bIforI=1,2,…nn.\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

−Δu=finΩu=0on∂Ω\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∈Vu\in V such that

−∫ΩΔ(u)v dv=∫Ωfv dv\begin{equation} - \int_\Omega \Delta(u) v \dv= \int_\Omega f v \dv \tag{12}\end{equation}

for all v∈Vv\in V where V=H01(Ω)V=H^1_0(\Omega).

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

∫ΩΔ(u)v dv=∫Ω∇⋅(∇(u)v) dv−∫Ω∇u⋅∇v dv=∫∂Ωv(∇u⋅n) da⏟v=0 on ∂Ω−∫Ω∇u⋅∇v dv\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:

a(u,v)=∫Ω∇u⋅∇v dvb(v)=∫Ωf v dv\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 A\BA as

AI ⁣J=a(NJ,NI)=∫Ω∇NJ⋅∇NI dv=∫ΩBJ⋅BI dv\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

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

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

bI=∫ΩfNI dv\begin{equation} \begin{aligned} b^I &= \int_{\Omega} f N^I \dv \end{aligned} \tag{17}\end{equation}