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 \sigma(\mathbf{u}) &= \mathbf{0} && \text{in } \Omega \\ \mathbf{u} &= \mathbf{0} && \text{on } \Gamma_D \\ \sigma(\mathbf{u}) \, \mathbf{n} &= \mathbf{f} && \text{on } \Gamma_N \end{aligned} \right. \]

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

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 \sigma(\mathbf{u}) : e(\mathbf{v}) \, dx = \int_{\Gamma_N} \mathbf{f} \cdot \mathbf{v} \, ds \]

Implementation

Rodin provides a convenient LinearElasticityIntegral that assembles the bilinear form for the Hooke tensor. The Lamé parameters $ \lambda $ and $ \mu $ are passed as arguments:

#include <Rodin/Solver.h>
#include <Rodin/Geometry.h>
#include <Rodin/Variational.h>
#include <Rodin/LinearElasticity.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 lambda = 0.5769, mu = 0.3846;

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

  VectorFunction load{0, -1};  // Downward force

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

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

  return 0;
}

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

Bilinear forms involve both the trial and test functions. Common examples:

Mathematical ExpressionRodin Code
$ \int_\Omega \nabla u \cdot \nabla v \, dx $ Integral(Grad(u), Grad(v))
$ \int_\Omega u \, v \, dx $ Integral(u, v)
$ \int_\Omega \lambda \nabla u \cdot \nabla v \, dx $ Integral(lambda * Grad(u), Grad(v))
$ \int_\Omega \lambda (\nabla \cdot u)(\nabla \cdot v) \, dx $ Integral(lambda * Div(u), Div(v))

Linear Forms

Linear forms involve only the test function. Common examples:

Mathematical ExpressionRodin Code
$ \int_\Omega f \, v \, dx $ Integral(f, v)
$ \int_{\Gamma_N} g \, v \, ds $ BoundaryIntegral(g, v).over(GammaN)
$ \int_\Omega \mathbf{f} \cdot \mathbf{v} \, dx $ Integral(f, v)

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(mesh, d, 2);  // Velocity: vector H1, order 2
H1 Qh(mesh, 1, 1);  // 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 »