Lecture 02 – Introduction to the Finite Element Method

Notation

For simplicity, we wrote PDEs in vector form, but they can also be expressed in component, e.g. \[ -\nabla \cdot (k\nabla u) = -\partial_x(k \partial_x u) - \partial_y(k \partial_y u) - \partial_z(k \partial_z u). \]

Other helpful differential operators include:

• Gradient: \(\nabla u = \left(\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}, \frac{\partial u}{\partial z}\right)\)
• Divergence: \(\nabla \cdot \mathbf{v} = \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} + \frac{\partial v_z}{\partial z}\) with \(\mathbf{v} = (v_x, v_y, v_z)\)
• Laplacian: \(\Delta u = \nabla^2 u = \nabla \cdot \nabla u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}\)
• Rotation: \(\nabla \times \mathbf{v} = \left(\frac{\partial v_z}{\partial y} - \frac{\partial v_y}{\partial z}, \frac{\partial v_x}{\partial z} - \frac{\partial v_z}{\partial x}, \frac{\partial v_y}{\partial x} - \frac{\partial v_x}{\partial y}\right)\)

Boundary-value-problems (BVPs)

• Domain \(\Omega\) with boundary \(\Gamma\), the domain \(\Omega\) is where the PDE holds. (Mathematically, a domain is an open, non-empty, connected subset of \(\mathbb{R}^n\) (e.g., a 2D region or 3D volume), and the boundary is its edge)
• \(\Gamma = \partial \Omega\) is where we specify boundary conditions, \(n\) is the outward normal vector on \(\Gamma\).

Boundary conditions:

• Dirichlet: \(u = u_D\) on \(\Gamma_D\) (specify value of \(u\))
• Neumann: \(-k\nabla u \cdot n = g_N\) on \(\Gamma_N\) (specify flux across boundary)
• Robin: \(\alpha u + \beta (-k\nabla u \cdot n) = g_R\) on \(\Gamma_R\) (combination of value and flux)

A BVP: \[ \begin{cases} -\nabla \cdot (k\nabla u) = f & \text{in } \Omega, \\ u = u_D & \text{on } \Gamma_D, \\ -k\nabla u \cdot n = g_N & \text{on } \Gamma_N. \end{cases} \]

Exact solution

For some simpler BVPs, we can find an exact solution analytically. \[ -\partial_x(k \partial_x u) = f(x) \quad \text{on } (0,1), \qquad u(0)=u(1)=0. \] with constant \(k \in \mathbb{R}, k \neq 0\) and \(f = \sin(\pi x)\) has exact solution:

\[ u(x) = \tfrac{\sin(\pi x)}{k\pi^2}. \]

Since \(\partial_x(k \partial_x u) = \partial_x\left(k \frac{\pi \cos(\pi x)}{k\pi^2}\right) = -\sin(\pi x)\), and the BCs are satisfied.

import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 1, 100)
k = 1.0
u = np.sin(np.pi * x) / (k * np.pi**2)
plt.plot(x, u, label='Exact solution')
plt.xlabel('x')
plt.ylabel('u(x)')
plt.title('Exact solution of 1D Poisson equation')
# change size
plt.gcf().set_size_inches(6, 3)
plt.grid(); plt.legend(); plt.show()

More complex domains

• For structured domains, the finite difference method (see exercise) can sometimes help.
• However, for more complex domains (e.g., with holes, curved boundaries), we need a more flexible method.

Finite element method (FEM): Can handle complex domains, unstructured meshes, and variable coefficients.

Mesh example

Variational problem

• Define trial space \(V\) and test space \(V_0\)
  • Sobolev space \(H^1(\Omega)\) include functions with square-integrable derivatives, i.e., \(H^1(\Omega) = \{ w\in L^2(\Omega)\,|\,\nabla w \in (L^2(\Omega))^d \}\) meaning \(w\) and its gradient are square-integrable, i.e., \(\int_\Omega |w|^2\,dx < \infty\) and \(\int_\Omega |\nabla w|^2\,dx < \infty\). \[ V = \{ w\in H^1(\Omega)\,|\,w=u_D\text{ on }\Gamma_D\}, \qquad V_0 = \{ w\in H^1(\Omega)\,|\,w=0\text{ on }\Gamma_D\}. \]

Resulting problem: Find \(u\in V\) such that

\[ a(u,v)=L(v)\quad \forall v\in V_0, \]

with a form \(a : V \times V_0 \to \mathbb{R}\) and a linear form \(L: V_0 \to \mathbb{R}\) defined as

\[ a(u,v)=\int_\Omega k\nabla u\cdot\nabla v\,dx,\qquad L(v)=\int_\Omega fv\,dx+\int_{\Gamma_N}g_N v\,ds. \]

The energy viewpoint (homogeneous BCs, \(k=1\))

• Start from an energy functional: \(E(u) = \frac{1}{2}\int_\Omega |\nabla u|^2\,dx - \int_\Omega fu\,dx\).
• Gateaux derivative: \(E'(u;v) = \lim_{\epsilon\to 0} \frac{E(u+\epsilon v) - E(u)}{\epsilon} = \int_\Omega \nabla u \cdot \nabla v\,dx - \int_\Omega fv\,dx\).
• Energy minimization (calculus of variations): \(\min_{u\in H_0^1(\Omega)} E(u)\) leads to \(E'(u;v) = 0\) for all “direction” \(v\in V_0\) as the first-order optimality condition, which is exactly the weak form.

\[ \begin{aligned} E'(u;v) &= \lim_{\epsilon\to 0} \frac{E(u+\epsilon v) - E(u)}{\epsilon} \\ &= \lim_{\epsilon\to 0} \frac{1}{\epsilon} \left( \frac{1}{2} \int_\Omega \nabla (u+\epsilon v) \cdot \nabla (u+\epsilon v)\,dx - \int_\Omega f(u+\epsilon v)\,dx - \left(\int_\Omega |\nabla u|^2\,dx - \int_\Omega fu\,dx\right) \right) \\ &= \int_\Omega \nabla u \cdot \nabla v\,dx - \int_\Omega fv\,dx \overset{!}{=} 0 \quad \forall v\in V_0. \end{aligned} \]

See [Brezzi & Fortin, 1991]

Implementation: Galerkin

• Start from a variational form: find \(u \in V_0\) such that \(a(u,v)=L(v) \forall v\in V_0\).
• Galerkin: choose \(V_h \subset V_0\) and solve \(a(u_h,v_h)=L(v_h) \quad \forall v_h\in V_h\).
• Therefore, \(u_h = \sum_{j=1}^N U_j\phi_j\) with basis functions \(\{\phi_j\}\) of \(V_h\) and unknown coefficients \(U_j\).
• Inseration into the variational form and testing with basis \(\phi_i\) leads to: \[ \sum_{j=1}^N U_j a(\phi_j, \phi_i) = L(\phi_i) \quad \forall i=1,\ldots,N, \]

Elements and forms

• FEniCS supports a wide range of finite elements (Lagrange, Raviart–Thomas, Nédélec, etc.), see BASIX
• UFL allows to express variational forms (linear, nonlinear, time-dependent), for example:

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = k * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx + g_N * v * ulf.ds

Other operators from UFL: ufl.div, ufl.curl, ufl.cross, ufl.inner, ufl.outer, etc.

Einstein summation convention is supported: Dx(v, i)*Dx(w, i) means \(\sum_{i=1}^d \frac{\partial v}{\partial x_i} \frac{\partial w}{\partial x_i}\).

Algorithmic Workflow

1. Generate (or load) mesh \(\mathcal T_h\) and boundary markers.
2. Build finite element space \(V_h\).
3. Define variational forms \(a(\cdot,\cdot), L(\cdot)\).
4. Assemble sparse matrix/vector.
5. Apply boundary conditions.
6. Solve with direct or iterative solver.
7. [Estimate error / refine mesh  postprocess]