Partial differential equations » Solving the Poisson equation

Resolution of the Poisson equation using Rodin.

Image

Introduction

To introduce the Rodin::Variational module it is important to start from a type of "Hello World!" example. In this context, the Poisson equation fits exactly our purpose. Hence, in this example we will seek to solve:

\[ \left\{ \begin{aligned} - \Delta u &= f && \mathrm{in} \ \Omega \\ u &= g && \mathrm{on} \ \Gamma \end{aligned} \right. \]

where $ \Omega $ is a domain with boundary $ \Gamma := \partial \Omega $ . To make matters easier we will set $ f = 1 $ and $ g = 0 $ .

In general, when solving an equation like this we usually utilize its variational formulation. Multiplying by a test function $ v $ and integrating by parts we may write the variational formulation

\[ \text{Find } u \in H^1_0(\Omega) \text{ s.t. } \forall v \in H^1_0(\Omega), \quad \int_{\Omega} \lambda \nabla u \cdot \nabla v \ dx = \int_{\Omega} f v \ dx \]

Obtaining the variational (or weak) formulation of a partial differential equation is a very standard procedure in the Finite Element Method literature. Hence we won't go into specifics. Rather, we deem it best to show how to implement said weak formulation and solve it.

Includes and Namespaces

In order to start using Rodin we first need to include the relevant parts which will aid us in the resolution of our equation. If you are familiar with finite elements then you will know that in particular we will need three main components:

  • A mesh approximating our domain $ \Omega $ .
  • A solver object that allows us to solve the problem
  • A way to specify the problem

Hence we import the required modules in Rodin:

#include "Rodin/Mesh.h"
#include "Rodin/Solver.h"
#include "Rodin/Variational.h"

using namespace Rodin;
using namespace Rodin::Variational;

The first part includes the following modules:

While the second part will simply introduce all the objects inside the Rodin and Rodin::Variational namespaces to the global scope, so we don't have to write their namespace again.

Implementing the solution

Next we would like to write the code which actually solves the equation. For our purposes we will write it inside the main function like so:

int main()
{
  // Code solving the equation
}

One of the first things you usually do is to keep track of the segment of the boundary where we impose the boundary conditions. In this case, we have homogenous Dirichlet boundary conditions on the boundary elements with attribute equal to 1. Additionally, we load the mesh from some file given its filename:

int Gamma = 1;
Mesh Omega = Mesh::load(filename);

With the mesh loaded, we can build our finite element space $ V_h $ and create some functions belonging to it. For our purposes, we need that our space contain functions belonging to the $ H^1 (\Omega) $ collection of functions. This is achieved by the following line of code:

H1 Vh(Omega);

After creating this space, we construct trial and test functions $ u $ and $ v $ which will be used to define our variational problem:

TrialFunction u(Vh);
TestFunction  v(Vh);

Hence we can start implementing the weak formulation associated to the Poisson equation. In general, Rodin expects that you specify the weak formulation by moving all terms to the left hand side of the equation. In symbols, this would look like:

\[ \text{Find } u \in H^1_0(\Omega) \text{ s.t. } \forall v \in H^1_0(\Omega), \quad a(u, v) - b(v) = 0 \ , \]

where $ a(u, v) $ and $ b(v) $ are our respective bilinear and linear forms. The left hand side is then used to define the variational problem which we want to solve. In the code this is achieved by the following declarations:

auto f = ScalarFunction(1.0);
auto g = ScalarFunction(0.0);

Problem poisson(u, v);
poisson = Integral(Grad(u), Grad(v))
        - Integral(f, v)
        + DirichletBC(u, g).on(Gamma);

In this case the Integral(Grad(u), Grad(v)) object represents the bilinear term

\[ \int_\Omega \nabla u \cdot \nabla v \ dx \ , \]

while the Integral(f, v) object represents the linear term

\[ \int_\Omega f v \ dx . \]

Lastly, the line DirichletBC(u, g).on(Gamma) specifies that the solution will satisfy the Dirichlet condition:

\[ u = g \quad \text{on} \quad \Gamma . \]

With the problem defined, we now need to solve the underlying linear system given by the associated stiffness matrix. There are many different types of solvers in Rodin but for this example we will utilize the Conjugate Gradient (CG) solver, as shown in the following code snippet:

Solver::CG().setMaxIterations(200)
            .setRelativeTolerance(1e-12)
            .printIterations(true)
            .solve(poisson);

This piece of code creates a CG object in place, while setting some parameters, and finally calling the solve method on the problem.

Ultimately, with the problem solved, we will want to visualize our solution. We can do this by saving both the solution and mesh objects to file, and then using GLVis to view the results.

u.getGridFunction().save("u.gf");
Omega.save("Omega.mesh");

To visualize the results, we can call GLVis from the command line:

glvis -m Omega.mesh -g u.gf

The full source code may be found below:

Partial differential equations | Solving the elasticity equation »