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:

\[ \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)

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 $ \| \nabla \mathbf{u} \|_F $ ScalarFrobenius(Jacobian(u))
Component $ u_i $ Scalar (from vector $ \mathbf{u} $ )Component(u, i) or u.x()
// Laplacian-like term: ∫ ∇u · ∇v dx
Integral(Grad(u), Grad(v))

// Divergence term (for vector fields): ∫ (∇·u)(∇·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 examples/DG/Poisson.cpp 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

« Finite element spaces | General concepts | Weak Formulations »