General concepts » Variational formulations in Rodin

How to define, assemble, and solve variational problems.

Introduction

The finite element method works by converting a partial differential equation (PDE) from its strong form into a weak (variational) form. Rodin's Variational module provides a form language that lets you express these variational formulations in C++ code that closely mirrors the mathematical notation.

From Strong to Weak Form

Consider a generic second-order PDE in its strong form:

with boundary conditions on . The standard procedure to obtain the weak form is:

Multiply by a test function
Integrate over the domain
Integrate by parts to reduce the order of derivatives on
Apply boundary conditions

This produces a variational problem of the form:

where is a bilinear form and is a linear form**.

Example: The Poisson Equation

The Poisson equation with on has the weak form:

In Rodin, this translates directly to:

Problem poisson(u, v);
poisson = Integral(Grad(u), Grad(v))   // a(u,v) = integral of grad(u) . grad(v)
        - Integral(f, v)                // l(v)   = integral of f * v
        + DirichletBC(u, Zero());       // u = 0 on the boundary

Forms in Rodin

Bilinear Forms

A bilinear form involves both the trial function and the test function . In Rodin, bilinear forms are created using Integral with both trial and test function arguments:

// Stiffness: integral of grad(u) . grad(v)
Integral(Grad(u), Grad(v))

// Mass: integral of u * v
Integral(u, v)

// Scaled: integral of k * u * v
Integral(k * u, v)

Linear Forms

A linear form involves only the test function :

// Load: integral of f * v
Integral(f, v)

// Boundary load: integral of g * v over a boundary
BoundaryIntegral(g, v).over(GammaN)

Problem Assembly

The Problem class assembles bilinear and linear forms into a linear system :

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

Problem problem(u, v);
problem = Integral(Grad(u), Grad(v))   // bilinear form (goes into A)
        - Integral(f, v)                // linear form (goes into b)
        + DirichletBC(u, Zero());       // modifies A and b

Note the sign convention: bilinear forms are added with +, while linear forms are subtracted with -. This corresponds to moving everything to the left-hand side:

Differential Operators

Rodin provides standard differential operators that work on trial and test functions:

Operator	Mathematical Meaning	Code
Gradient		`Grad(u)`
Divergence		`Div(u)`
Jacobian	(matrix)	`Jacobian(u)`
Transpose		`Jacobian(u).T()`

// Laplacian-like term
Integral(Grad(u), Grad(v))

// Divergence term (for vector fields)
Integral(Div(u), Div(v))

// Elasticity strain energy
Integral(Jacobian(u) + Jacobian(u).T(), Jacobian(v) + Jacobian(v).T())

Boundary Conditions

Dirichlet Conditions

Dirichlet** (essential) boundary conditions fix the value of the solution on part of the boundary:

// Homogeneous: u = 0 on entire boundary
DirichletBC(u, Zero())

// Non-homogeneous: u = g on entire boundary
DirichletBC(u, g)

// On specific boundary segment (by attribute)
DirichletBC(u, Zero()).on(1)

// Vector-valued: u = (0, 0) on a segment
DirichletBC(u, VectorFunction{0, 0}).on(GammaD)

Neumann Conditions

Neumann** (natural) boundary conditions specify the flux or traction on part of the boundary. They appear as boundary integrals in the variational formulation:

// Scalar Neumann: integral of g * v over boundary segment
BoundaryIntegral(g, v).over(GammaN)

// Vector Neumann (traction): integral of f . v over boundary segment
BoundaryIntegral(f, v).over(GammaN)

Solving the Problem

After defining a problem, use a solver from the Solver module to solve the resulting linear system:

// Conjugate Gradient (for SPD systems)
Solver::CG(problem).solve();

// Direct sparse LU factorization
Solver::SparseLU(problem).solve();

The solution is then accessible via the trial function:

auto& solution = u.getSolution();

Saving and Visualizing Results

After solving, export the results using XDMF for visualization in ParaView:

#include <Rodin/IO/XDMF.h>

// Export to XDMF
IO::XDMF xdmf("Results");
xdmf.grid().setMesh(mesh).add("u", u.getSolution());
xdmf.write();

Complete Example

Here is a complete example that defines and solves the Poisson equation, then exports the result using XDMF:

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

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

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

  // 2. Finite element space
  P1 Vh(mesh);
  TrialFunction u(Vh);
  TestFunction  v(Vh);

  // 3. Source term
  RealFunction f = 1;

  // 4. Variational problem
  Problem poisson(u, v);
  poisson = Integral(Grad(u), Grad(v))
          - Integral(f, v)
          + DirichletBC(u, Zero());

  // 5. Solve
  CG(poisson).solve();

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

  return 0;
}

Mixed Formulations

Some problems require solving for multiple unknowns simultaneously, such as velocity and pressure in the Navier-Stokes equations. Rodin's Problem class supports mixed formulations with multiple trial/test function pairs:

// Taylor-Hood finite element spaces
H1 uh(std::integral_constant<size_t, 2>{}, mesh, d);  // P2 velocity
H1 ph(std::integral_constant<size_t, 1>{}, mesh);      // P1 pressure

TrialFunction u(uh);  // Velocity
TrialFunction p(ph);  // Pressure
TestFunction  v(uh);
TestFunction  q(ph);

// Mixed problem with multiple trial/test pairs
Problem stokes(u, p, v, q);
stokes = nu * Integral(Jacobian(u), Jacobian(v))
       - Integral(p, Div(v))
       + Integral(Div(u), q)
       - Integral(f, v)
       + DirichletBC(u, g);