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.

The cross-module rules behind these examples are summarized in Numerical contracts, including the Problem sign convention and the distinction between interpolation and L2 projection.

From Strong to Weak Form

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

\[ \mathcal{L} u = f \quad \text{in } \Omega \]

with boundary conditions on $ \partial\Omega $ . The standard procedure to obtain the weak form is:

  1. Multiply by a test function $ v $
  2. Integrate over the domain $ \Omega $
  3. Integrate by parts to reduce the order of derivatives on $ u $
  4. Apply boundary conditions

This produces a variational problem of the form:

\[ \text{Find } u \in V \text{ s.t. } \forall v \in V, \quad a(u, v) = \ell(v) \]

where $ a(u, v) $ is a bilinear form and $ \ell(v) $ is a linear form.

Example: The Poisson Equation

The Poisson equation $ -\Delta u = f $ with $ u = 0 $ on $ \partial\Omega $ has the weak form:

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

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 $ a(u, v) $ involves both the trial function $ u $ and the test function $ v $ . 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)

Named cell forms

For repeated assembly of common cell operators, Rodin provides MassForm, DiffusionForm, HelmholtzForm, and LinearElasticityForm. They assemble the same discrete operators, including quadrature orders, as the corresponding integrals. The first assembly constructs the sparse pattern; later assemblies reuse its map to matrix values. The named forms currently use local sparse matrices with Sequential or OpenMP assembly.

MassForm mass(c, u, v);                 // Integral(c * u, v)
DiffusionForm diffusion(a, u, v);       // Integral(a * Grad(u), Grad(v))
HelmholtzForm helmholtz(a, c, u, v);    // Diffusion and mass terms
LinearElasticityForm elasticity(lambda, mu, w, z);

Here u and v use scalar spaces for diffusion and Helmholtz; mass also accepts vector spaces. Elasticity uses vector-valued w and z whose dimension matches the mesh. Each form can restrict its cell integral with .over(attribute). The elasticity form is declared by Rodin/Solid.h; the other three are declared by Rodin/Variational.h.

Linear Forms

A linear form $ \ell(v) $ involves only the test function $ v $ :

// 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 $ A\mathbf{x} = \mathbf{b} $ :

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 both A and b

Sign convention explained:

The mathematical variational problem is:

\[ a(u, v) = \ell(v) \]

Rodin moves everything to the left-hand side, giving:

\[ a(u, v) - \ell(v) = 0 \]

This is why bilinear forms (involving both $ u $ and $ v $ ) are added with +, while linear forms (involving only $ v $ ) are subtracted with -. The Problem object recognizes which terms are bilinear (containing both trial and test) and which are linear (containing only the test function), and routes them to the matrix $ A $ or vector $ \mathbf{b} $ accordingly.

Differential Operators

Rodin provides standard differential operators that work on trial and test functions. These are expression templates: they are not evaluated immediately, but recorded as symbolic expressions and evaluated later during assembly.

OperatorMathematical meaningResultCode
Gradient $ \nabla u $ Vector (for scalar $ u $ )Grad(u)
Divergence $ \nabla \cdot \mathbf{u} $ Scalar (for vector $ \mathbf{u} $ )Div(u)
Jacobian $ \nabla \mathbf{u} $ (matrix) $ d \times d $ matrixJacobian(u)
Transpose $ (\nabla \mathbf{u})^T $ $ d \times d $ matrixTranspose(Jacobian(u))
Trace $ \text{tr}(\nabla \mathbf{u}) $ ScalarTrace(Jacobian(u))
Frobenius norm $ \lVert \nabla \mathbf{u} \rVert_F $ ScalarFrobenius(Jacobian(u))
Component $ u_i $ Scalar (from vector $ \mathbf{u} $ )Component(u, i) or u.x()
// Laplacian-like term: \int \nabla u \cdot \nabla v dx
Integral(Grad(u), Grad(v))

// Divergence term (for vector fields): \int (\nabla\cdot u)(\nabla\cdot v) dx
Integral(Div(u), Div(v))

// Elasticity: symmetric gradient ε(u) : ε(v)
Integral(Jacobian(u) + Transpose(Jacobian(u)),
         Jacobian(v) + Transpose(Jacobian(v)))

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();

// Or pass a standalone solver to the problem
Solver::CG solver(problem);
solver.setTolerance(1e-10);
solver.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);

Discontinuous Galerkin Formulations

For problems with discontinuous finite element spaces (or when penalty methods are needed across element interfaces), Rodin provides specialized integral types and operators:

DG Operators

OperatorDescriptionCode
JumpJump across a face: $ [u] = u^+ - u^- $ Jump(u)
AverageAverage across a face: $ \{u\} = \frac{1}{2}(u^+ + u^-) $ Average(u)

DG Integral Types

Integral typeRegionDescription
Integral(...)Cell interiors $ K $ Standard domain integral
BoundaryIntegral(...)Boundary faces $ \partial\Omega $ Boundary-only integration
FaceIntegral(...)All faces $ \mathcal{F}_h $ Interior + boundary faces
InterfaceIntegral(...)Interior faces $ \mathcal{F}_h^i $ Interior faces only

DG Poisson Example

The simplest DG approach uses P0 (piecewise constant) functions. The actual DG Poisson example uses this pattern:

P0 Vh(mesh);           // Discontinuous piecewise-constant space
TrialFunction u(Vh);
TestFunction  v(Vh);

RealFunction f = 1.0;

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

Solver::SparseLU(poisson).solve();

For higher-order DG methods, the symmetric interior penalty (SIP) formulation uses explicit jump and average operators:

\[ \sum_K \int_K \nabla u \cdot \nabla v \, dx - \int_{\mathcal{F}_h^i} \{\nabla u\} \cdot [v] \, ds - \int_{\mathcal{F}_h^i} [u] \cdot \{\nabla v\} \, ds + \frac{\sigma}{h} \int_{\mathcal{F}_h^i} [u] \cdot [v] \, ds = \int_\Omega f v \, dx \] 
Real sigma = 10.0;  // penalty parameter

Problem dg(u, v);
dg = Integral(Grad(u), Grad(v))
   - InterfaceIntegral(Average(Grad(u)), Jump(v))
   - InterfaceIntegral(Jump(u), Average(Grad(v)))
   + InterfaceIntegral(sigma * Jump(u), Jump(v))
   - Integral(f, v)
   + DirichletBC(u, Zero());

See the DG Poisson example for a complete working implementation.

See Also

« Matrix storage and finite-element component ordering | General concepts | Weak Formulations »