The goal of this chapter is to show how the Poisson equation, the most basic of all PDEs, can be quickly solved with a few lines of FEniCS code. We introduce the most fundamental FEniCS objects such asMesh,Function,FunctionSpace,TrialFunction, andTestFunction, and learn how to write a basic PDE solver, including the specification of the mathematical variational problem, applying boundary conditions, calling the FEniCS solver, and plotting the solution.
Most books on a programming language start with a "Hello, World!" program. That is, one is curious about how a very fundamental task is expressed in the language, and writing a text to the screen can be such a task. In the world of finite element methods for PDEs, the most fundamental task must be to solve the Poisson equation. Our counterpart to the classical "Hello, World!" program therefore solves $$ \begin{align} - \nabla^2 u(\x) &= f(\x),\quad \x\mbox{ in } \Omega, \tag{2.1}\\ u(\x) &= \ub(\x),\quad \x\mbox{ on } \partial \Omega\tp \tag{2.2} \end{align} $$ Here, \( u = u(\x) \) is the unknown function, \( f = f(\x) \) is a prescribed function, \( \nabla^2 \) is the Laplace operator (also often written as \( \Delta \)), \( \Omega \) is the spatial domain, and \( \partial\Omega \) is the boundary of \( \Omega \). A stationary PDE like this, together with a complete set of boundary conditions, constitute a boundary-value problem, which must be precisely stated before it makes sense to start solving it with FEniCS.
In two space dimensions with coordinates \( x \) and \( y \), we can write out the Poisson equation as $$ \begin{equation} - {\partial^2 u\over\partial x^2} - {\partial^2 u\over\partial y^2} = f(x,y)\tp \tag{2.3} \end{equation} $$ The unknown \( u \) is now a function of two variables, \( u = u(x,y) \), defined over a two-dimensional domain \( \Omega \).
The Poisson equation arises in numerous physical contexts, including heat conduction, electrostatics, diffusion of substances, twisting of elastic rods, inviscid fluid flow, and water waves. Moreover, the equation appears in numerical splitting strategies for more complicated systems of PDEs, in particular the Navier–Stokes equations.
Solving a PDE such as the Poisson equation in FEniCS consists of the following steps:
FEniCS is based on the finite element method, which is a general and efficient mathematical machinery for numerical solution of PDEs. The starting point for the finite element methods is a PDE expressed in variational form. Readers who are not familiar with variational problems will get a very brief introduction to the topic in this tutorial, but reading a proper book on the finite element method in addition is encouraged. The section The finite element method contains a list of some suitable books. Experience shows that you can work with FEniCS as a tool to solve your PDEs even without thorough knowledge of the finite element method, as long as you get somebody to help you with formulating the PDE as a variational problem.
The basic recipe for turning a PDE into a variational problem is to multiply the PDE by a function \( v \), integrate the resulting equation over the domain \( \Omega \), and perform integration by parts of terms with second-order derivatives. The function \( v \) which multiplies the PDE is called a test function. The unknown function \( u \) to be approximated is referred to as a trial function. The terms test and trial function are used in FEniCS programs too. Suitable function spaces must be specified for the test and trial functions. For standard PDEs arising in physics and mechanics such spaces are well known.
In the present case, we first multiply the Poisson equation by the test function \( v \) and integrate over \( \Omega \): $$ \begin{equation} \tag{2.4} -\int_\Omega (\nabla^2 u)v \dx = \int_\Omega fv \dx\tp \end{equation} $$ A common rule when we derive variational formulations is that we try to keep the order of the derivatives of \( u \) and \( v \) as low as possible (this will enlarge the collection of finite elements that can be used in the problem). Here, we have a second-order spatial derivative of \( u \), which can be transformed to a first-derivative of \( u \) and \( v \) by applying the technique of integration by parts. A Laplace term will always be subject to integration by parts . The formula reads $$ \begin{equation} \tag{2.5} -\int_\Omega (\nabla^2 u)v \dx = \int_\Omega\nabla u\cdot\nabla v \dx - \int_{\partial\Omega}{\partial u\over \partial n}v \ds , \end{equation} $$ where \( \frac{\partial u}{\partial n} = \nabla u \cdot n \) is the derivative of \( u \) in the outward normal direction \( n \) on the boundary.
3: Integration by parts in more than one space dimension is based on Gauss' divergence theorem. Simply take (2.5) as the formula to be used.
Another feature of variational formulations is that the test function \( v \) is required to vanish on the parts of the boundary where the solution \( u \) is known (the book [10] explains in detail why this requirement is necessary). In the present problem, this means that \( v=0 \) on the whole boundary \( \partial\Omega \). The second term on the right-hand side of (2.5) therefore vanishes. From (2.4) and (2.5) it follows that $$ \begin{equation} \int_\Omega\nabla u\cdot\nabla v \dx = \int_\Omega fv \dx\tp \tag{2.6} \end{equation} $$ If we require that this equation holds for all test functions \( v \) in some suitable space \( \hat V \), the so-called test space, we obtain a well-defined mathematical problem that uniquely determines the solution \( u \) which lies in some (possibly different) function space \( V \), the so-called trial space. We refer to (2.6) as the weak form or variational form of the original boundary-value problem (2.1)--(2.2).
The proper statement of our variational problem now goes as follows: Find \( u \in V \) such that $$ \begin{equation} \tag{2.7} \int_{\Omega} \nabla u \cdot \nabla v \dx = \int_{\Omega} fv \dx \quad \forall v \in \hat{V}\tp \end{equation} $$ The trial and test spaces \( V \) and \( \hat V \) are in the present problem defined as $$ \begin{align*} V &= \{v \in H^1(\Omega) : v = \ub \mbox{ on } \partial\Omega\}, \\ \hat{V} &= \{v \in H^1(\Omega) : v = 0 \mbox{ on } \partial\Omega\}\tp \end{align*} $$
In short, \( H^1(\Omega) \) is the mathematically well-known Sobolev space containing functions \( v \) such that \( v^2 \) and \( |\nabla v|^2 \) have finite integrals over \( \Omega \) (essentially meaning that the functions are continuous). The solution of the underlying PDE must lie in a function space where also the derivatives are continuous, but the Sobolev space \( H^1(\Omega) \) allows functions with discontinuous derivatives. This weaker continuity requirement of \( u \) in the variational statement (2.7), as a result of the integration by parts, has great practical consequences when it comes to constructing finite element function spaces. In particular, it allows the use of piecewise polynomial function spaces; i.e., function spaces constructed by stitching together polynomial functions on simple domains such as intervals, triangles, or tetrahedrons.
The variational problem (2.7) is a continuous problem: it defines the solution \( u \) in the infinite-dimensional function space \( V \). The finite element method for the Poisson equation finds an approximate solution of the variational problem (2.7) by replacing the infinite-dimensional function spaces \( V \) and \( \hat{V} \) by discrete (finite-dimensional) trial and test spaces \( V_h\subset{V} \) and \( \hat{V}_h\subset\hat{V} \). The discrete variational problem reads: Find \( u_h \in V_h \subset V \) such that $$ \begin{equation} \tag{2.8} \int_{\Omega} \nabla u_h \cdot \nabla v \dx = \int_{\Omega} fv \dx \quad \forall v \in \hat{V}_h \subset \hat{V}\tp \end{equation} $$ This variational problem, together with a suitable definition of the function spaces \( V_h \) and \( \hat{V}_h \), uniquely define our approximate numerical solution of Poisson's equation (2.1). The mathematical framework may seem complicated at first glance, but the good news is that the finite element variational problem (2.8) looks the same as the continuous variational problem (2.7), and FEniCS can automatically solve variational problems like (2.8)!
It turns out to be convenient to introduce the following canonical notation for variational problems: $$ \begin{equation} a(u, v) = L(v)\tp \tag{2.9} \end{equation} $$ For the Poisson equation, we have: $$ \begin{align} a(u, v) &= \int_{\Omega} \nabla u \cdot \nabla v \dx, \tag{2.10}\\ L(v) &= \int_{\Omega} fv \dx\tp \tag{2.11} \end{align} $$ From the mathematics literature, \( a(u,v) \) is known as a bilinear form and \( L(v) \) as a linear form. We shall in every linear problem we solve identify the terms with the unknown \( u \) and collect them in \( a(u,v) \), and similarly collect all terms with only known functions in \( L(v) \). The formulas for \( a \) and \( L \) are then coded directly in the program.
FEniCS provides all the necessary mathematical notation needed to express the variational problem \( a(u, v) = L(v) \). To solve a linear PDE in FEniCS, such as the Poisson equation, a user thus needs to perform only two steps: