Solving the Poisson equation

From a boundary-value problem to a finite element solution and visualization.

This example solves the Poisson equation on the unit square with a constant source and zero values on the boundary. It constructs a triangular mesh, expresses the weak form in Rodin, solves the resulting linear system, and exports the solution for visualization.

The mathematical problem

Let . The unknown scalar field satisfies

The source is positive throughout the domain, while the homogeneous Dirichlet condition fixes the field to zero along every edge. The resulting solution is positive in the interior and vanishes on the boundary.

To obtain the weak form, multiply the equation by a test function and integrate over the domain. Integration by parts transfers one derivative from to ; the boundary term vanishes because has zero trace. The weak problem seeks such that

Rodin approximates this problem in a finite-dimensional space of continuous, piecewise linear functions.

1. Construct the mesh

A uniform grid with 16 vertices in each coordinate direction is split into triangles. UniformGrid uses unit spacing, so the coordinates are scaled by to obtain the unit square. The resulting mesh has 256 vertices and 450 triangles.

  constexpr size_t nodes = 16;
  Mesh mesh = Mesh().UniformGrid(Polytope::Type::Triangle, {nodes, nodes});
  mesh.scale(1.0 / (nodes - 1));
  mesh.getConnectivity().compute(1, 2);

The final line computes edge-to-cell connectivity. This allows Rodin to identify the exterior edges on which the boundary condition is imposed.

2. Define the finite element problem

P1 defines the space of functions that are continuous across the mesh and affine on each triangle. The trial function u represents the unknown, and the test function v represents the arbitrary test function in the weak form.

  P1 vh(mesh);
  TrialFunction u(vh);
  TestFunction v(vh);
  RealFunction f = 1;

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

The first Integral is the stiffness term, corresponding to the left-hand side of the weak form. The second is the source term with . DirichletBC(u, Zero()) imposes zero values on the entire exterior boundary. Rodin writes the equation as a residual, so the source term is subtracted from the stiffness term; setting that residual to zero recovers the weak form.

3. Solve the linear system

After the boundary values are imposed, the Poisson system is symmetric positive definite. The conjugate gradient solver is therefore suitable for this problem. solve() assembles and solves the system, and success() checks that the solve completed successfully before the result is exported.

  CG solver(poisson);
  solver.solve();
  if (!solver.success())
    return EXIT_FAILURE;

The computed finite element field is available through u.getSolution().

4. Export and inspect the result

The solution and its mesh are written in XDMF format. The accompanying HDF5 files hold the numerical data and must remain alongside Poisson.xdmf.

  IO::XDMF xdmf("Poisson");
  xdmf.grid().setMesh(mesh).add("u", u.getSolution());
  xdmf.write().flush();

Computed solution on the unit square. Height and color represent the scalar field.

In ParaView:

Open Poisson.xdmf and select Apply.
Color the mesh by the point field u.
Apply Warp By Scalar using u to display the field as a surface.
Select Surface With Edges to make the triangulation visible.

On this mesh, the computed nodal values range from on the boundary to approximately in the interior. The surface shown above uses the same solution; ParaView amplifies its height by a factor of eight for visibility.

Complete program

The complete source combines the steps above and can be built with the other Rodin examples. The output files are written to the working directory.

See Poisson.cpp for the source, the installation guide for build setup, and the I/O guide for additional export options.

Partial differential equations | Solving diffusion with discontinuous Galerkin »