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:
- 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.
- 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.
- 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.
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 ψ, ψ, ..., ψ[N − 1] be the basis functions let e, e, ..., 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