General concepts » Your First Problem: Solving the Poisson Equation

A complete walkthrough of solving a PDE with Rodin.

In this section, we'll solve a complete partial differential equation using Rodin: the Poisson equation. This is often the "Hello World" of finite element methods.

The Mathematical Problem

We want to solve the following boundary value problem:

\[ \left\{ \begin{aligned} -\Delta u &= f && \text{in } \Omega \\ u &= 0 && \text{on } \partial\Omega \end{aligned} \right. \]

where:

  • $ \Omega $ is our domain (a square [0,1] × [0,1])
  • $ u $ is the unknown function we're solving for
  • $ f = 1 $ is a known source term
  • $ \partial\Omega $ is the boundary of the domain

Weak Formulation

The weak (or variational) formulation of this problem is:

Find $ u \in H^1_0(\Omega) $ such that for all $ v \in H^1_0(\Omega) $ :

\[ \int_{\Omega} \nabla u \cdot \nabla v \, dx = \int_{\Omega} f v \, dx \]

This is the form that finite element methods work with, and it's exactly what we'll implement in Rodin.

The Complete Code

Here's the complete program to solve the Poisson equation:

#include <Rodin/Geometry.h>
#include <Rodin/Variational.h>
#include <Rodin/Solver.h>

using namespace Rodin;
using namespace Rodin::Geometry;
using namespace Rodin::Variational;
using namespace Rodin::Solver;

int main() {
  // 1. Create a mesh
  Mesh mesh;
  mesh = mesh.UniformGrid(Polytope::Type::Triangle, {16, 16});

  // 2. Compute connectivity for boundary operations
  mesh.getConnectivity().compute(1, 2);

  // 3. Create P1 finite element space
  P1 Vh(mesh);

  // 4. Define trial and test functions
  TrialFunction u(Vh);
  TestFunction v(Vh);

  // 5. Define the source term f = 1
  RealFunction f = 1.0;

  // 6. Define the problem
  Problem poisson(u, v);
  poisson = Integral(Grad(u), Grad(v))
          - Integral(f, v)
          + DirichletBC(u, Zero());

  // 7. Solve the problem
  CG(poisson).solve();

  // 8. Save the solution
  u.getSolution().save("solution.gf");
  mesh.save("mesh.mesh");

  return 0;
}

Code Breakdown

Let's understand each part of the code:

Step 1-2: Create the Mesh

Mesh mesh;
mesh = mesh.UniformGrid(Polytope::Type::Triangle, {16, 16});
mesh.getConnectivity().compute(1, 2);
  • Creates a 2D triangular mesh with a 16×16 grid
  • Computes face-to-cell connectivity (dimension 1 to dimension 2)
  • This connectivity is needed to identify boundary elements

Step 3-4: Define Function Space

P1 Vh(mesh);
TrialFunction u(Vh);
TestFunction v(Vh);
  • P1 creates a space of piecewise linear functions
  • u is the unknown function (trial function)
  • v is the test function for the weak formulation

Step 5: Define Source Term

RealFunction f = 1.0;
  • Creates a constant function $ f = 1 $
  • Could also be a more complex function: RealFunction f([](const Point& p) { return p.x() * p.y(); })

Step 6: Formulate the Problem

Problem poisson(u, v);
poisson = Integral(Grad(u), Grad(v))
        - Integral(f, v)
        + DirichletBC(u, Zero());

This directly translates the weak formulation:

  • Integral(Grad(u), Grad(v)) represents $ \int_{\Omega} \nabla u \cdot \nabla v \, dx $
  • Integral(f, v) represents $ \int_{\Omega} f v \, dx $
  • DirichletBC(u, Zero()) enforces $ u = 0 $ on the boundary

Note how the code closely resembles the mathematical notation!

Step 7: Solve the System

CG(poisson).solve();
  • Uses the Conjugate Gradient (CG) iterative solver
  • For symmetric positive definite systems (like the Poisson equation)
  • Alternative solvers: SparseLU, BiCGSTAB, etc.

Step 8: Save Results

u.getSolution().save("solution.gf");
mesh.save("mesh.mesh");
  • Saves the solution in MFEM GridFunction format
  • Saves the mesh in MEDIT format
  • Can be visualized with tools like GLVis

Building and Running

Save the code as poisson.cpp and create a CMakeLists.txt:

cmake_minimum_required(VERSION 3.16)
project(PoissonExample CXX)

set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

find_package(Rodin REQUIRED)

add_executable(poisson poisson.cpp)
target_link_libraries(poisson PRIVATE
  Rodin::Geometry
  Rodin::Variational
  Rodin::Solver
)

Build and run:

cmake .
make
./poisson

Visualizing the Solution

If you have GLVis installed, visualize the solution:

glvis -m mesh.mesh -g solution.gf

You should see a smooth solution that represents the deflection of a membrane under uniform load.

Common Variations

Non-Homogeneous Boundary Conditions

To set $ u = g $ on the boundary instead of $ u = 0 $ :

RealFunction g([](const Point& p) {
  return p.x() * p.x();
});

poisson = Integral(Grad(u), Grad(v))
        - Integral(f, v)
        + DirichletBC(u, g);

Partial Boundary Conditions

To apply boundary conditions only on specific boundaries:

// Mark boundary with attribute 1
int boundaryAttribute = 1;

poisson = Integral(Grad(u), Grad(v))
        - Integral(f, v)
        + DirichletBC(u, Zero()).on(boundaryAttribute);

Variable Source Term

For a spatially varying source:

RealFunction f([](const Point& p) {
  return std::sin(M_PI * p.x()) * std::sin(M_PI * p.y());
});

poisson = Integral(Grad(u), Grad(v))
        - Integral(f, v)
        + DirichletBC(u, Zero());

Using Different Solvers

// Sparse LU direct solver (for small to medium problems)
SparseLU().solve(poisson);

// Conjugate Gradient with options
CG solver;
solver.setMaxIterations(1000);
solver.setRelativeTolerance(1e-10);
solver.printIterations(true);
solver.solve(poisson);

What's Next?

Congratulations! You've solved your first PDE with Rodin. To deepen your understanding, proceed to Understanding Core Concepts to learn more about the underlying concepts.

You can also explore:

« First Steps with Rodin | General concepts | Understanding Core Concepts »