Boundary Conditions
Guide to applying boundary conditions in Rodin.
Introduction
Boundary conditions specify the behavior of the solution on the boundary of the computational domain. They are essential for obtaining a well-posed problem.
In the finite element method, boundary conditions are classified into two main types:
- Essential (Dirichlet) conditions: Directly prescribe the value of the solution on the boundary. These are imposed by modifying the finite element space or the assembled linear system.
- Natural (Neumann) conditions: Prescribe the flux or derivative of the solution on the boundary. These appear naturally in the weak formulation as boundary integrals.
A third type, periodic conditions, ties degrees of freedom on opposite sides of the domain.
Dirichlet Boundary Conditions
Dirichlet conditions specify the value of the solution on the boundary:
Homogeneous Dirichlet (@f$ g = 0 @f$)
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)
Neumann Boundary Conditions
Neumann conditions prescribe the normal derivative (or flux) on the boundary:
These arise naturally from integration by parts in the weak formulation and appear as boundary integrals:
// ∫_ΓN 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 ( ) 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 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):
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
| Type | Math | Rodin Code | Appears In |
|---|---|---|---|
| Dirichlet | on | DirichletBC(u, g).on(GammaD) | System modification |
| Neumann | on | BoundaryIntegral(h, v).over(GammaN) | Weak formulation |
| Homogeneous Neumann | (nothing needed) | Default | |
| Periodic | PeriodicBC(u, dofMap) | System modification |