PDE Project Course

1. Theory of finite element methods

Johan Jansson

johanjan@math.chalmers.se

Department of Computational Mathematics

Goal of the lecture

• Refresh or give knowledge and notation of theoretical aspects of finite element methods which are relevant for carrying out the projects.

Lecture plan

• Introduction to FEM

• FEM for Poisson’s equation

• Adaptivity for Poisson’s equation

Introduction to FEM

PDE Project Course 04/05 – p.4/25

A method for solving PDEs

The finite element method (FEM), also known as Galerkin’s method, is a general method for solving PDEs (or ODEs) of the form

A(u(x)) = f(x), x ∈ Ω

where A is a differential operator, f is a given force term, u is the solution and Ω is the domain.

Solving PDEs

• Analytic solutions can be obtained only for simple geometries in special cases:

PSfrag replacements −Δu = 0

• Using the computer, we can obtain solutions to general problems with complex geometries:

\begin{aligned}u̇ + u · ∇u − νΔu + ∇p & = & f \\ ∇ · u & = & 0 \end{aligned}

The finite element method

Find an approximate solution U of the form

U(x) = ∑_{j=1}^{N} ξ_{j}φ_{j}(x).

Here U is linear linear combination of (a finite number of) basis functions with local support:

{φ_{j}}_{j=1}^{N}.

PDE Project Course 04/05 – p.7/25

Basis functions

• Typically piecewise linear nodal basis functions. Other bases also possible.

• We define a mesh T_{h} = {K} “triangulating” the domain Ω into elements K (intervals, triangles, tetrahedrons, ...), where the vertices of the elements form nodes N. h is a function defined by h(x) = h_{K}, where h_{K} is the diameter of the ball circumscribing K.

Basis functions

• Nodal basis function defined by:

ϕ_{j}(N_{i}) = \begin{cases}1 & , i = j, \\ 0 & , i ≠ j \end{cases}

Linear between nodes.

Basis functions 1D

• Local basis functions on the reference element (an interval): K_{r} = [0, 1]

• Local basis functions:

\begin{aligned}ϕ_{0}(x) & = & 1 − x \\ ϕ_{1}(x) & = & x \end{aligned}

• The global basis function ϕ_{j}(x) on any element K can be computed via a mapping F : Ω_{0} → Ω. ϕ_{j}(x) on Ω is the union of ϕ_{j}(x) on all K.

PDE Project Course 04/05 – p.10/25

Basis functions 2D

• Local basis functions on the reference element (a triangle): K_{r} = {(0, 0), (1, 0), (0, 1)} • Local basis functions:

\begin{aligned}ϕ_{0}(x) & = & 1 − x − y \\ ϕ_{1}(x) & = & x \\ ϕ_{2}(x) & = & y \end{aligned}

Some notation from functional analysis

• Scalar product for functions v, w: (v, w) = ∫_{Ω} v(x)w(x) dx

• L_{2}(Ω)-norm of a function v:

∥v∥_{L_{2}(Ω)} = (∫_{Ω} v^{2} dx)^{1/2} = \sqrt{(v, v)}

Some notation from functional analysis

• Cauchy’s inequality:

|(v, w)| ≤ ∥v∥ ∥w∥

• v and w are orthogonal iff (v, w) = 0

PSfrag replacements

v
w

PDE Project Course 04/05 – p.13/25

Galerkin’s method

It is often advisable to rewrite the differential equation A(u) = f from operator form to variational form:

a(u, v) = (f, v) ∀v ∈ V,

where a(·, ·) = (A(·), ·) is a bilinear form, and V is a suitable function space. Concretely:

∫_{Ω} A(u)v dx = ∫_{Ω} fv dx ∀v ∈ V

Galerkin’s method

The finite element method is based on Galerkin’s method:

• Let V_{h} denote a finite dimensional trial space.

• Let V̂_{h} denote a finite dimensional test space.

• Find U ∈ V_{h} such that the residual R(U) = A(U) − f is orthogonal to V̂_{h}:

\begin{aligned}(R(U), v) & = & 0 ∀v ∈ V̂_{h} ⇔ \\ a(U, v) − (f, v) & = & 0 ∀v ∈ V̂_{h} \end{aligned}

PDE Project Course 04/05 – p.15/25

Galerkin’s method

Starting from the variational formulation with V_{h} = V̂_{h} = span{φ_{j}}_{j=1}^{N} we have

\begin{aligned}a(U, v) − (f, v) & = & 0, & ∀v ∈ V_{h}, \\ a(∑_{j=1}^{N} ξ_{j}φ_{j}, v) & = & (f, v), & ∀v ∈ V_{h}, \\ ∑_{j=1}^{N} ξ_{j}a(φ_{j}, φ_{i}) & = & (f, φ_{i}), & i = 1, …, N, \end{aligned}

or

A_{h}ξ = b,

where A_{h} = a(φ_{j}, φ_{i}), b = (f, φ_{i}).

PDE Project Course 04/05 – p.16/25

FEM for Poisson’s equation

PDE Project Course 04/05 – p.17/25

Poisson in three different forms

• Equation:

−Δu = f

• Variational formulation: ∫_{Ω} ∇u · ∇v dx = ∫_{Ω} fv dx ∀v ∈ V

• Linear system:

A_{h} = ∫_{Ω} ∇φ_{j} · ∇φ̂_{i} dx, b = ∫_{Ω} fφ̂_{i} dx

PDE Project Course 04/05 – p.18/25

Details

• Green’s formula:

∫_{Ω} ∇u · ∇v dx = ∫_{Γ} ∂_{n}uv ds − ∫_{Ω} Δuv dx • Start from equation:

−Δu = f ⇒
− ∫_{Ω} Δuv dx = ∫_{Ω} fv dx ⇒
∫_{Ω} ∇u · ∇v dx − ∫_{Γ} ∂_{n}uv ds = ∫_{Ω} fv dx ⇒
∫_{Ω} ∇u · ∇v dx = ∫_{Ω} fv dx

PDE Project Course 04/05 – p.19/25

Details

• Continue

∫_{Ω} ∇u · ∇v dx = ∫_{Ω} fv dx
⇒ insert Ansatz
∑_{j=1}^{N} ξ_{j} ∫_{Ω} ∇φ_{j} · ∇φ_{i} dx = ∫_{Ω} fφ_{i} dx ⇒
A_{h}ξ = b

Adaptivity for Poisson

PDE Project Course 04/05 – p.21/25

How large is the error?

We expect the error e = U − u to decrease if we increase the dimension N of V_{h} and V̂_{h}. This can be done in different ways:

• h-adaptivity: decrease the mesh size h

• p-adaptivity: increase the polynomial order p

• hp-adaptivity: a combination of the h and p methods

We will only consider h-adaptivity.

An a posteriori error estimate

Let ∥ · ∥_{E} denote the energy-norm given by ∥v∥_{E} = ∥∇v∥. Then the (piecewise linear) finite element solution U = U(x) satisfies the error estimate

∥e∥_{E} = ∥U − u∥_{E} ≤ C∥h(R_{1}(U) + R_{2}(U))∥,

where R_{1}(U) = |f + ΔU| = |f| and

R_{2}(U) = \frac{1}{2} max_{S⊂∂K} h_{K}^{−1}|[∂_{S}U]|.

Adaptive error control

Find V_{h}, given by a triangulation T_{h}, such that

∥e∥_{E} ≤ TOL,

where TOL is a given tolerance for the error. This is satisfied if

C∥h(R_{1}(U) + R_{2}(U))∥ ≤ TOL.

An adaptive algorithm

1. Choose an initial triangulation T_{h}^{0}.

2. Compute the solution U on the current triangulation.

3. Compute the residuals R_{1}, R_{2}, and the error estimate.

4. If the error estimate is below the tolerance we stop. Otherwise, we refine the elements where R_{1} + R_{2} is large and start again at 2.