The finite element method

This is an edited section from my thesis “Solving Poisson’s equation with Dataflow computing”. It might be a bit more formal than most of my other posts.


The finite element method is a method to solve boundary value problems. Boundary value problems pose the problem of finding a function that satisfies a given differential equation on and given boundary conditions on the boundary of the domain. We will use Poisson’s problem on a two-dimensional domain with Dirichlet boundary conditions:

For a given , , and , find a function such that

The represents a differential operator called the Laplace operator. For twice-differentiable functions on a domain , the Laplace operator is defined as

The finite element method is a numerical method. This means that the solution is approximated by a function . The function is taken to be in some finite-dimensional approximation space . The space is spanned by a basis . When the basis functions , …, are fixed, can be expressed by a vector :

The elements are known as degrees of freedom.

Conceptually, the finite element method can be separated into several stages:

  1. A method to convert the boundary value problem to a discrete system of symbolic equations. The most popular method is the Bubnov-Galerkin method, which is often simply called the Galerkin method.
  2. The discretization strategy, which describes how to obtain an actual system of equations for a given symbolic system of equations. This involves picking the basis functions, and evaluating or approximating the integrals in the symbolic system. Also, one or more methods of refinement are usually provided, which can be used to increase the number of basis functions and thus the dimension of , so that can approximate better.
  3. The solver, which solves the system of equations. If the boundary value problem is linear, the system of equations will be linear too, and linear solvers like the conjugate gradient method, the generalized minimum-residual method, or the stabilized biconjugate gradient method are popular choices. Nonlinear systems are harder to solve, and require more advanced methods like the Newton-Rhapson method or the Picard method. Typical nonlinear solvers iteratively linearize the system, and solve the linearized system with a linear solver.

Galerkin method

The Galerkin method derives a system of symbolic equations for a given boundary value problem. The Galerkin method is based on the weak formulation of the boundary value problem. First, let us define the Sobolev space of functions for which the first partial derivatives are square-integrable. The weak formulation of a differential equation can be obtained by multiplying both sides by a test function , integrating over , and demanding equality regardless of the choice of . Applying this to gives

One can then apply integration by parts and use to make the left-hand side symmetric:

With this equation, the complete weak formulation becomes:

For a given , , and , find a function such that:

The Galerkin method replaces by a finite-dimensional approximation , spanned by basis functions , and seeks an approximate solution . Ignoring the boundary conditions for now, we get the system

So can be expressed as . Using this, we get

Bringing the summation outside the integral yields

Which can be written as a linear system of size , where the finite element matrix is defined as

and the right-hand side vector is defined as

Furthermore, we define the stiffness matrix as

and the mass matrix as

We can now express the finite element matrix in terms of the stiffness matrix and mass matrix as

The Dirichlet boundary conditions can be implemented by a Dirichlet lift. Each degree of freedom , , …, which corresponds to a basis function which is nonzero on the boundary is prescribed in such a way that on . As such, the boundary condition is approximately satisfied. To pick in such a way, one can use projection on the boundary, or make sure that interpolates in a set of points , , … on the boundary.


Suppose that the mass matrix is to be assembled. One could separately evaluate the integrals . However, the basis functions are usually defined on sections of the domain called elements. The domain can be partitioned into elements so that we have

The usual approach is to loop over the elements and evaluate for the basis functions which are nonzero on that element. Typically, the basis functions will have local support, and there will be a small number of nonzero basis functions on each element. For an element on which there are nonzero basis functions , an element matrix with entries is assembled. After the element matrix is assembled, each entry is added to .

Pseudocode for the assembly is as follows. Note that there are many details omitted (for example, how to evaluate the integral in the assembly of the element matrix).

let ψ[0], ψ[1], ..., ψ[N − 1] be the basis functions
let e[0], e[1], ..., e[M − 1] be the elements
let glob_mat be an N × N matrix with each element equal to zero
for k = 0 .. M − 1
	let n be the number of basis functions with support on e[k]
	let s[k, 0], s[k, 1], ..., s[k, n − 1] be the indices of the basis
		functions with support on e[k]
	let elem_mat be an n × n matrix with each element equal to zero
	// assemble element matrix elem_mat
	for i = 0 .. n − 1
		for j = 0 .. n − 1
			elem_mat[i, j] += integral of ψ[s[k, i]] ψ[s[k, j]] over e[k]
		end for
	end for
	// scatter element matrix elem_mat to global matrix glob_mat
	for i = 0 .. n − 1
		for j = 0 .. n − 1
			glob_mat[s[k, i], s[k, j]] += elem_mat[i, j]
		end for
	end for
end for