General concepts » Boundary Conditions

Guide to applying boundary conditions in Rodin.

Introduction

Boundary conditions specify the behavior of the solution on the boundary $ \partial \Omega $ of the computational domain. They are essential for obtaining a well-posed problem — without appropriate boundary conditions, the PDE may have no solution, or infinitely many solutions.

In the finite element method, boundary conditions are classified by how they enter the discrete system:

  • Essential (Dirichlet) conditions: Directly prescribe the value of the solution on the boundary: $ u = g $ on $ \Gamma_D $ . These are imposed by modifying the assembled linear system — the corresponding rows and columns of the system matrix are adjusted to enforce the prescribed values. In Rodin, this is handled by DirichletBC.
  • Natural (Neumann) conditions: Prescribe the flux or normal derivative of the solution: $ \partial u / \partial n = h $ on $ \Gamma_N $ . These appear naturally in the weak formulation as boundary integrals after integration by parts. In Rodin, they are expressed as BoundaryIntegral(h, v).over(GammaN).
  • Periodic conditions: Tie degrees of freedom on opposite boundaries. In Rodin, handled by PeriodicBC with an explicit DOF correspondence map.

Dirichlet Boundary Conditions

Dirichlet conditions specify the value of the solution on the boundary:

\[ u = g \quad \text{on } \Gamma_D \subseteq \partial\Omega \]

Homogeneous Dirichlet (

The most common case sets the solution to zero on the boundary:

Problem problem(u, v);
problem = Integral(Grad(u), Grad(v))
        - Integral(f, v)
        + DirichletBC(u, Zero());       // u = 0 on entire boundary

Non-Homogeneous Dirichlet

To prescribe a non-zero value:

// Using a constant
DirichletBC(u, RealFunction(1.0))

// Using a spatial function
RealFunction g([](const Geometry::Point& p) {
  return p.x() * p.x();
});
DirichletBC(u, g)

Partial Boundary Conditions

To apply Dirichlet conditions only on a specific part of the boundary, use the .on() method with a boundary attribute:

// Apply u = 0 only on boundary with attribute 1
DirichletBC(u, Zero()).on(1)

// Apply u = g only on boundary with attribute 2
DirichletBC(u, g).on(2)

// Apply on multiple boundary attributes at once
DirichletBC(u, Zero()).on({1, 3})

Boundary attributes are integers associated with boundary elements in the mesh. They are typically set during mesh generation or loaded from a mesh file.

Vector Dirichlet Conditions

For vector-valued problems (e.g., elasticity), use vector functions:

// Fix displacement to zero (clamped boundary)
DirichletBC(u, VectorFunction{0, 0}).on(GammaD)

// Prescribe a specific displacement
VectorFunction prescribed{0.1, 0};
DirichletBC(u, prescribed).on(GammaD)

Sparse algebraic elimination

Math::LinearSystem also accepts an arbitrary caller-supplied square CSC operator. Essential elimination does not require symmetric coefficients, stored transpose entries, or an existing constrained diagonal. Duplicate triplets must first be assembled into valid sorted, unique-row sparse storage, as with Eigen's setFromTriplets.

Architecture

Let $C$ denote the prescribed indices, $F$ the remaining indices, and $g_C$ their prescribed coefficients. These are finite-element coefficients, not necessarily nodal samples at high order. The free equations follow by substituting $x_C=g_C$ into the original system:

\[ A_{FF}x_F=b_F-A_{FC}g_C. \]

Define the diagonal constraint projector $Q$ , its complement $P=I-Q$ , and extend $g_C$ by zero on $F$ . Keeping the original system size gives

\[ \widetilde A=PAP+Q,\qquad \widetilde b=P(b-Ag)+g. \]

This leaves $A_{FF}$ unchanged and replaces constrained equations by identities. The formula applies over both real and complex scalar fields; prescribing a coefficient does not conjugate it. Symmetry or Hermitian symmetry is preserved when present but is not required for correctness.

For example, prescribing $x_0=1$ in

\[ A=\begin{pmatrix}2&3\\0&4\end{pmatrix},\qquad b=\begin{pmatrix}5\\6\end{pmatrix} \]

produces $\widetilde A=\operatorname{diag}(1,4)$ and $\widetilde b=(1,6)^T$ , whether or not the zero at $(1,0)$ is explicitly stored. Discovering constrained rows through their transpose entries would miss the coefficient at $(0,1)$ in this example.

Implementation

The operation validates the system dimensions and shifted constraint indices before mutation, then compresses the operator. It traverses every constrained column of the original matrix to accumulate the load lifting, recording its diagonal position after the traversal. Only after all lifting is complete may an absent diagonal be inserted. Missing diagonals are merged together into the existing CSC arrays by moving entries backwards, avoiding repeated shifts or a full matrix copy. Existing entries and explicit zeros are retained.

A temporary byte mask identifies constrained rows. Consecutive constrained columns are cleared in coalesced storage ranges. Free storage between them is examined in blocks: an empty row mask performs no coefficient stores, while a nonempty mask visits only its set bits. The final scalar tail follows the same zeroing rule. Actual zero assignment or bit clearing is required: multiplication by zero would retain NaNs or create them from infinities. Constraint diagonals and prescribed load entries are restored last.

A single contiguous constraint interval can select sorted row-range pruning. The selection uses sixteen samples at equal stored-entry quantiles, requires at least eight stored entries per column on average, and requires at least half the usable samples to reject the interval. These are measured heuristics, not correctness assumptions. Samples only select the clearing kernel; actual column bounds certify every omission. Partial overlaps use binary searches to locate the constrained coefficient range. Fragmented constraints retain the block kernel. No particular FEM boundary numbering is assumed.

For $n$ rows, $m$ stored entries and $k$ constraints, block clearing and batch insertion require $O(n+m+k)$ work and $O(n+k)$ temporary storage. Pruning adds up to two binary searches per overlapping free column. Operator storage grows only for absent constrained diagonals. It remains compressed, with sorted row indices, after the call. The solution vector is untouched, preserving warm starts. Repeating the same constraints is idempotent up to floating-point rounding. An empty constraint map only validates dimensions and compresses storage.

This algorithm is serial and constructs no persistent transpose index or parallel worker pool. Its linear scan replaces repeated transpose searches; it does not make a stronger asymptotic claim than linear stored-entry work. Tests compare the matrix and load against independent projector constructions, including asymmetric storage, missing diagonals, complex data, offsets and mixed finite-element systems. The opt-in RodinSparseEliminationBenchmarks executable measures the public operation with restoration outside timing.

Neumann Boundary Conditions

Neumann conditions prescribe the normal derivative (or flux) on the boundary:

\[ \frac{\partial u}{\partial n} = g \quad \text{on } \Gamma_N \]

These arise naturally from integration by parts in the weak formulation and appear as boundary integrals:

// \int_\GammaN g v ds
BoundaryIntegral(g, v).over(GammaN)

Example: Neumann on One Side

int GammaD = 1;  // Dirichlet boundary
int GammaN = 2;  // Neumann boundary

RealFunction g = 1.0;  // Prescribed flux

Problem problem(u, v);
problem = Integral(Grad(u), Grad(v))
        - Integral(f, v)
        - BoundaryIntegral(g, v).over(GammaN)
        + DirichletBC(u, Zero()).on(GammaD);

Homogeneous Neumann (Natural)

Homogeneous Neumann conditions ( $ \partial u / \partial n = 0 $ ) are the "default" — they are automatically satisfied when no boundary integral is added and no Dirichlet condition is imposed on that part of the boundary.

Vector Neumann (Traction)

In elasticity, a Neumann condition corresponds to an applied traction (force per unit area). The traction vector $ \mathbf{f} $ is prescribed on the loaded boundary:

VectorFunction traction{0, -1};   // Downward pull

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

Periodic Boundary Conditions

Periodic conditions enforce that the solution on one boundary equals the solution on the opposite boundary. They are used for problems on domains that tile (e.g. unit cells in homogenization):

\[ u(0, y) = u(2\pi, y), \quad u(x, 0) = u(x, 2\pi) \]

In Rodin, periodic conditions are imposed by building an IndexMap that pairs the DOF indices on opposite boundaries, then passing it to PeriodicBC:

const size_t n = 200;
Mesh mesh;
mesh = mesh.UniformGrid(Polytope::Type::Triangle, { n, n });
mesh.getConnectivity().compute(1, 2);
mesh.scale(1.0 / (n - 1));
mesh.scale(2 * M_PI);   // Scale to [0, 2π]^2

P1 vh(mesh);

// Build DOF correspondence: left ↔ right, bottom ↔ top
IndexMap<IndexSet> dofs;
for (Index i = 0; i < vh.getSize(); i += n)
  dofs[i].insert(i + n - 1);        // left ↔ right
for (Index i = 0; i < n; i++)
  dofs[i].insert(i + n * (n - 1));  // bottom ↔ top

TrialFunction u(vh);
TestFunction  v(vh);

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

See the periodic Poisson example for a complete walkthrough.

Mixed Boundary Conditions

In practice, different parts of the boundary often have different conditions. For example, in elasticity:

int GammaD = 1;  // Clamped (fixed) boundary
int GammaN = 2;  // Loaded boundary
// Remaining boundary: traction-free (homogeneous Neumann)

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

Complete Example

Here is a complete example with mixed boundary conditions:

#include <Rodin/Solver.h>
#include <Rodin/Geometry.h>
#include <Rodin/Variational.h>
#include <Rodin/IO/XDMF.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 problem(u, v);
  problem = Integral(Grad(u), Grad(v))
          - Integral(f, v)
          + DirichletBC(u, Zero());

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

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

  return 0;
}

Summary Table

TypeMathRodin CodeAppears In
Dirichlet $ u = g $ on $ \Gamma_D $ DirichletBC(u, g).on(GammaD)System modification
Neumann $ \partial u / \partial n = h $ on $ \Gamma_N $ BoundaryIntegral(h, v).over(GammaN)Weak formulation
Homogeneous Neumann $ \partial u / \partial n = 0 $ (nothing needed)Default
Periodic $ u(x_L) = u(x_R) $ PeriodicBC(u, dofMap)System modification

See Also

« Weak Formulations | General concepts | Solvers »