Partial differential equations » Solving PDEs on surface meshes

Reaction-diffusion equation on a 3D surface mesh.

Introduction

Rodin can solve PDEs not only on 2D planar domains or 3D volumes, but also on surface meshes — 2D manifolds embedded in 3D space. This is useful for modeling phenomena on shells, biological membranes, or curved interfaces.

This example solves a reaction-diffusion steady-state problem on a surface mesh:

where is the Laplace–Beltrami operator (the surface Laplacian) and is a curved surface in .

The weak formulation is identical to the planar case — Rodin handles the surface geometry automatically through the mesh metric:

Implementation

Loading the Surface Mesh

The surface mesh is loaded from an MFEM file. The mesh has topological dimension 2 (triangles) but space dimension 3 (vertices live in ):

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

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

Mesh Omega;
Omega.load("surface-meshes-example.mesh", IO::FileFormat::MFEM);

Finite Element Space

We use H1 elements (higher-order continuous Lagrange elements) to better approximate the solution on the curved surface:

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

Right-Hand Side

The source term is defined in terms of the 3D coordinates of points on the surface:

auto f = RealFunction(
  [](const Point& x)
  {
    double l2 = x(0) * x(0) + x(1) * x(1) + x(2) * x(2);
    return 7 * x(0) * x(1) / l2;
  });

Problem Assembly

The problem is assembled exactly as in the planar case. The Integral(Grad(u), Grad(v)) automatically uses the surface gradient:

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

Problem rd(u, v);
rd = Integral(Grad(u), Grad(v))    // Surface stiffness: ∫_Γ ∇_Γ u · ∇_Γ v ds
   + Integral(u, v)                // Mass term: ∫_Γ u v ds
   - Integral(f, v);               // Load: ∫_Γ f v ds
rd.solve(cg);

Saving Results

u.getSolution().save("u.gf", IO::FileFormat::MFEM);
Omega.save("Omega.mesh", IO::FileFormat::MFEM);

Key Observations

No special operators needed: Rodin automatically computes the surface gradient from the mesh geometry. The code is identical to the planar case.
H1 elements on surfaces: Higher-order elements are recommended on curved surfaces for better geometric approximation.
Boundary-free manifolds: Closed surfaces have no boundary, so no boundary conditions are required. The mass term provides coercivity.
Surface integrals: The Integral(...) operator automatically uses the surface measure instead of .

Full Source Code

The complete source code is:

SurfaceMesh.cpp