General concepts » Weak Formulations

Understanding and implementing weak formulations in Rodin.

Introduction

The weak formulation (also called variational formulation) is the foundation of the finite element method. Instead of solving a PDE in its classical (strong) form, we reformulate it in a weaker sense that is more amenable to numerical approximation.

The general idea is to multiply the PDE by a test function, integrate over the domain, and apply integration by parts. This process:

  • Reduces the regularity requirements on the solution
  • Naturally incorporates certain boundary conditions
  • Leads to a well-posed problem suitable for discretization

Derivation Procedure

Consider a general second-order PDE on a domain $ \Omega $ :

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

To derive the weak formulation:

  1. Multiply by a test function $ v $ from an appropriate space
  2. Integrate over $ \Omega $ : $ \int_\Omega (\mathcal{L} u) \, v \, dx = \int_\Omega f \, v \, dx $
  3. Integrate by parts to transfer derivatives from $ u $ to $ v $
  4. Apply boundary conditions (essential or natural)

Example: The Poisson Equation

Strong Form

The Poisson equation with homogeneous Dirichlet boundary conditions:

\[ \left\{ \begin{aligned} -\Delta u &= f && \text{in } \Omega \\ u &= 0 && \text{on } \partial\Omega \end{aligned} \right. \]

Derivation

Multiply by $ v \in H^1_0(\Omega) $ and integrate:

\[ -\int_\Omega \Delta u \, v \, dx = \int_\Omega f \, v \, dx \]

Apply Green's formula (integration by parts):

\[ \int_\Omega \nabla u \cdot \nabla v \, dx - \underbrace{\int_{\partial\Omega} \frac{\partial u}{\partial n} v \, ds}_{= 0 \text{ (since } v = 0 \text{ on } \partial\Omega)} = \int_\Omega f \, v \, dx \]

Weak Formulation

Find $ u \in H^1_0(\Omega) $ such that for all $ v \in H^1_0(\Omega) $ :

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

Implementation in Rodin

The implementation in Rodin closely mirrors the mathematical notation:

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

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

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

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

  RealFunction f = 1;

  Problem poisson(u, v);
  poisson = Integral(Grad(u), Grad(v))   // ∫ ∇u · ∇v dx
          - Integral(f, v)               // ∫ f v dx
          + DirichletBC(u, Zero());      // u = 0 on ∂Ω

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

  return 0;
}

Example: Linearized Elasticity

Strong Form

The linearized elasticity system governs the displacement $ \mathbf{u} : \Omega \to \mathbb{R}^d $ of an elastic body:

\[ \left\{ \begin{aligned} -\nabla \cdot \boldsymbol{\sigma}(\mathbf{u}) &= \mathbf{f} && \text{in } \Omega \\ \mathbf{u} &= \mathbf{0} && \text{on } \Gamma_D \\ \boldsymbol{\sigma}(\mathbf{u}) \, \mathbf{n} &= \mathbf{g} && \text{on } \Gamma_N \end{aligned} \right. \]

where $ \boldsymbol{\sigma}(\mathbf{u}) = 2\mu \, \boldsymbol{\varepsilon}(\mathbf{u}) + \lambda \, \mathrm{tr}(\boldsymbol{\varepsilon}(\mathbf{u})) \, \mathbf{I} $ is the Cauchy stress tensor (Hooke's law) and $ \boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T) $ is the linearized (infinitesimal) strain tensor.

Lamé Parameters

The Lamé parameters $ \lambda $ and $ \mu $ characterize the material. They are related to the more common Young's modulus $ E $ and Poisson's ratio $ \nu $ by:

\[ \lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}, \quad \mu = \frac{E}{2(1+\nu)} \]

For example, with $ E = 1 $ and $ \nu = 0.3 $ :

const Real E = 1.0, nu = 0.3;
const Real lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));  // ≈ 0.5769
const Real mu     = E / (2.0 * (1.0 + nu));                     // ≈ 0.3846

Derivation

Multiply by $ \mathbf{v} \in [H^1_0(\Omega)]^d $ , integrate, and use the divergence theorem. The boundary integral on $ \Gamma_N $ produces the traction term:

\[ \int_\Omega \boldsymbol{\sigma}(\mathbf{u}) : \boldsymbol{\varepsilon}(\mathbf{v}) \, dx = \int_\Omega \mathbf{f} \cdot \mathbf{v} \, dx + \int_{\Gamma_N} \mathbf{g} \cdot \mathbf{v} \, ds \]

Expanding Hooke's law, the bilinear form becomes:

\[ a(\mathbf{u}, \mathbf{v}) = \int_\Omega \bigl[ \lambda (\nabla \cdot \mathbf{u})(\nabla \cdot \mathbf{v}) + 2\mu \, \boldsymbol{\varepsilon}(\mathbf{u}) : \boldsymbol{\varepsilon}(\mathbf{v}) \bigr] \, dx \]

Implementation

Rodin provides LinearElasticityIntegral (in <Rodin/Solid.h>) which assembles the bilinear form for the Hooke tensor. The helper class LinearElasticityIntegral(u, v) returns a callable that accepts the Lamé parameters — either as constants or spatially-varying functions:

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

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

int main()
{
  Mesh mesh;
  mesh.load("beam.mesh", IO::FileFormat::MFEM);
  mesh.getConnectivity().compute(1, 2);

  const size_t d = mesh.getSpaceDimension();
  const Real E = 1.0, nu = 0.3;
  const Real lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
  const Real mu     = E / (2.0 * (1.0 + nu));

  P1 Vh(mesh, d);  // Vector-valued P1 space
  TrialFunction u(Vh);
  TestFunction  v(Vh);

  VectorFunction traction{0, -1};  // Traction force (stress vector) on Neumann boundary

  Problem elasticity(u, v);
  elasticity = LinearElasticityIntegral(u, v)(lambda, mu)
             - BoundaryIntegral(traction, v).over(GammaN)
             + DirichletBC(u, VectorFunction{0, 0}).on(GammaD);

  Solver::CG(elasticity).solve();

  return 0;
}

Beyond Small Deformations

For large deformations (finite strain), use the Solid module's hyperelastic framework with MaterialTangent and InternalForce integrators coupled with a NewtonSolver. Available constitutive laws include NeoHookean, SaintVenantKirchhoff, and MooneyRivlin.

The Newton linearization at each iterate $ \mathbf{u}^k $ solves:

\[ K(\mathbf{u}^k) \, \delta\mathbf{u}^k = -F_{\mathrm{int}}(\mathbf{u}^k) + F_{\mathrm{ext}} \]

where $ K $ is the material tangent and $ F_{\mathrm{int}} $ is the internal force vector.

#include <Rodin/Solid.h>
#include <Rodin/Solver/NewtonSolver.h>
#include <Rodin/Solver/SparseLU.h>

Solid::NeoHookean law(lambda, mu);

GridFunction u(Vh);  // accumulated displacement

Solid::MaterialTangent tangent(law, du, v);
tangent.setDisplacement(u);

Solid::InternalForce residual(law, v);
residual.setDisplacement(u);

Problem newton(du, v);
newton = tangent + residual
       - Integral(bodyForce, v)
       + DirichletBC(du, VectorFunction{Zero(), Zero()}).on(fixedBC);

SparseLU linearSolver(newton);
NewtonSolver solver(linearSolver);
solver.solve(u);  // updates u in place

See Solid module for the full API.

General Structure

Most weak formulations in Rodin follow this pattern:

\[ \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 (depends on both $ u $ and $ v $ )
  • $ \ell(v) $ is a linear form (depends only on $ v $ )

In Rodin, the problem is defined by moving everything to one side:

Problem problem(u, v);
problem = [bilinear terms] - [linear terms] + [boundary conditions];

Bilinear Forms

A bilinear form $ a(u, v) $ is a mapping that is linear in both arguments. It involves both the trial function $ u $ and the test function $ v $ , and is assembled into the system matrix $ A $ .

Common examples:

Mathematical expressionRodin codeName
$ \int_\Omega \nabla u \cdot \nabla v \, dx $ Integral(Grad(u), Grad(v))Stiffness
$ \int_\Omega u \, v \, dx $ Integral(u, v)Mass
$ \int_\Omega \lambda \nabla u \cdot \nabla v \, dx $ Integral(lambda * Grad(u), Grad(v))Weighted stiffness
$ \int_\Omega \lambda (\nabla \cdot \mathbf{u})(\nabla \cdot \mathbf{v}) \, dx $ Integral(lambda * Div(u), Div(v))Divergence penalization

Linear Forms

A linear form $ \ell(v) $ is a mapping that is linear in $ v $ only (it does not depend on the unknown $ u $ ). It is assembled into the right-hand side vector** $ \mathbf{b} $ .

Common examples:

Mathematical expressionRodin codeName
$ \int_\Omega f \, v \, dx $ Integral(f, v)Volume source
$ \int_{\Gamma_N} g \, v \, ds $ BoundaryIntegral(g, v).over(GammaN)Neumann load
$ \int_\Omega \mathbf{f} \cdot \mathbf{v} \, dx $ Integral(f, v)Vector body force

Boundary Conditions in Weak Form

Dirichlet (Essential) Conditions

Dirichlet conditions are imposed directly on the function space:

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

// u = g on specific boundary
DirichletBC(u, g).on(boundaryAttribute)

Neumann (Natural) Conditions

Neumann conditions appear as boundary integrals in the weak formulation:

// ∂u/∂n = g on boundary with attribute 2
BoundaryIntegral(g, v).over(2)

Periodic Boundary Conditions

Periodic conditions tie DOFs on opposite boundary segments using an IndexMap. See the periodic Poisson example for a complete walkthrough.

// Build a DOF correspondence map: left edge ↔ right edge
IndexMap<IndexSet> dofs;
for (Index i = 0; i < Vh.getSize(); i += n)
  dofs[i].insert(i + n - 1);

Problem problem(u, v);
problem = Integral(Grad(u), Grad(v))
        - Integral(f, v)
        + PeriodicBC(u, dofs);

Mixed Formulations

Some PDEs involve multiple unknowns from different spaces. The Navier–Stokes example uses Taylor–Hood elements with separate trial and test functions for velocity and pressure:

H1 Vh(std::integral_constant<size_t, 2>{}, mesh, d);  // Velocity: vector H1, order 2
H1 Qh(std::integral_constant<size_t, 1>{}, mesh);     // Pressure: scalar H1, order 1

TrialFunction u(Vh);  TestFunction v(Vh);   // Velocity pair
TrialFunction p(Qh);  TestFunction q(Qh);   // Pressure pair

Problem stokes(u, p, v, q);
stokes = Integral(Grad(u), Grad(v))     // viscous term
       - Integral(p, Div(v))            // pressure–velocity coupling
       + Integral(Div(u), q)            // incompressibility constraint
       - Integral(f, v)                 // body force
       + DirichletBC(u, uD);           // velocity BC

Well-Posedness

The Lax-Milgram theorem guarantees existence and uniqueness of the solution when:

  • The bilinear form $ a(\cdot, \cdot) $ is continuous and coercive** on the function space $ V $
  • The linear form $ \ell(\cdot) $ is continuous on $ V $

For the Poisson equation, the bilinear form $ a(u, v) = \int_\Omega \nabla u \cdot \nabla v \, dx $ is coercive on $ H^1_0(\Omega) $ by the Poincaré inequality.

For mixed formulations (saddle-point problems), the Brezzi conditions (inf-sup stability) replace coercivity. Taylor–Hood elements satisfy the inf-sup condition, making them a standard choice for Stokes-type problems.

See Also

« Variational formulations in Rodin | General concepts | Boundary Conditions »