General concepts » Understanding Core Concepts

The Mesh Object

A Mesh is the fundamental geometric object in Rodin. It represents a discretization of your domain into simple shapes (polytopes).

Key concepts:**

• Vertices (dimension 0): Points in space
• Edges/Faces (dimension 1 in 2D, 1-2 in 3D): Line segments or faces
• Cells (dimension 2 in 2D, 3 in 3D): Triangles, quads, tetrahedra, etc.

Mesh mesh;
mesh = mesh.UniformGrid(Polytope::Type::Triangle, {16, 16});

std::cout << "Dimension: " << mesh.getDimension() << std::endl;
std::cout << "Vertices: " << mesh.getVertexCount() << std::endl;
std::cout << "Cells: " << mesh.getCellCount() << std::endl;

Polytopes

A Polytope is a geometric shape. Rodin supports various polytope types:

// Create different mesh types
Mesh tri_mesh = Mesh().UniformGrid(Polytope::Type::Triangle, {8, 8});
Mesh tet_mesh = Mesh().UniformGrid(Polytope::Type::Tetrahedron, {4, 4, 4});

Connectivity

Connectivity describes the relationships between polytopes of different dimensions. You need to compute connectivity before certain operations.

Mesh mesh;
mesh = mesh.UniformGrid(Polytope::Type::Triangle, {16, 16});

// Face-to-cell connectivity (needed for boundary operations)
mesh.getConnectivity().compute(1, 2);

// Vertex-to-cell connectivity
mesh.getConnectivity().compute(0, 2);

// Cell-to-cell connectivity (for finding neighbors)
mesh.getConnectivity().compute(2, 2);

Common connectivity patterns:**

• (1, 2): Edges to cells - needed for boundary conditions in 2D
• (2, 3): Faces to cells - needed for boundary conditions in 3D
• (0, d): Vertices to cells - for vertex-based operations
• (d, d): Cell-to-cell adjacency

Mesh Iteration

You can iterate over polytopes of a specific dimension:

Mesh mesh;
mesh = mesh.UniformGrid(Polytope::Type::Triangle, {8, 8});

// Iterate over all cells
for (auto it = mesh.getCell(); it; ++it) {
  // Access cell information
  std::cout << "Cell " << it->getIndex() << std::endl;
}

// Iterate over all vertices
for (auto it = mesh.getVertex(); it; ++it) {
  Point p = it->getCoordinates();
  std::cout << "Vertex at (" << p.x() << ", " << p.y() << ")" << std::endl;
}

// Iterate over faces (edges in 2D)
for (auto it = mesh.getFace(); it; ++it) {
  // Access face information
}

Finite Element Spaces

A finite element space is a finite-dimensional subspace of a function space (such as or ), built by assembling local polynomial spaces on each cell of the mesh. In Rodin, finite element spaces are the bridge between the mesh geometry and the variational formulation of a PDE.

Rodin provides the following standard spaces:

Mesh mesh;
mesh = mesh.UniformGrid(Polytope::Type::Triangle, {16, 16});

// Piecewise constant space: 1 DOF per cell
P0 Wh(mesh);
std::cout << "P0 DOFs: " << Wh.getSize() << std::endl;  // = number of cells

// Piecewise linear space: 1 DOF per vertex
P1 Vh(mesh);
std::cout << "P1 DOFs: " << Vh.getSize() << std::endl;  // = number of vertices

// Quadratic H1 space (order 2)
H1 Qh(std::integral_constant<size_t, 2>{}, mesh);
std::cout << "H1<2> DOFs: " << Qh.getSize() << std::endl;

Trial and Test Functions

In the finite element method, we express problems in weak form using trial and test functions. These are abstract symbolic objects that represent functions in the finite element space — they do not store numerical values themselves.

P1 Vh(mesh);

// Trial function: represents the unknown solution u ∈ Vh
TrialFunction u(Vh);

// Test function: represents the weighting function v ∈ Vh
TestFunction v(Vh);

Mathematical background:**

The finite element method seeks an approximate solution in the finite-dimensional space by requiring the residual to be orthogonal to all test functions:

• The trial function (represented by TrialFunction) is the unknown being solved for. After solving, its coefficients are stored in the associated GridFunction, accessible via u.getSolution().
• The test function (represented by TestFunction) is used to generate the equations of the discrete system. Each basis function of produces one equation , yielding the linear system .

Differential Operators

Rodin provides differential operators that work on trial and test functions. These operators are expression templates: they are not evaluated immediately but are stored as symbolic expressions and evaluated during assembly.

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

// The gradient of a scalar function is a vector
auto grad_u = Grad(u);
auto grad_v = Grad(v);

// Used in integrals: the dot product is taken implicitly
auto stiffness = Integral(Grad(u), Grad(v));  // ∫ ∇u · ∇v dx

Integrals

The Integral function represents integration over the domain or boundary:

// Domain integral: ∫_Ω ∇u · ∇v dx
auto stiffness = Integral(Grad(u), Grad(v));

// Domain integral with function: ∫_Ω f v dx
RealFunction f = 1.0;
auto load = Integral(f, v);

// The dot product is implicit for vector quantities
auto inner_product = Integral(Grad(u), Grad(v));

Boundary Conditions

Rodin supports different types of boundary conditions:

Dirichlet (essential) boundary conditions:**

// Homogeneous: u = 0 on boundary
DirichletBC(u, Zero())

// Non-homogeneous: u = g on boundary
RealFunction g([](const Point& p) { return p.x(); });
DirichletBC(u, g)

// On specific boundary: u = 0 on boundary with attribute 1
DirichletBC(u, Zero()).on(1)

Neumann (natural) boundary conditions:**

// Neumann conditions appear as boundary integrals
RealFunction g = 1.0;
auto neumann = Integral(g, v).over(BoundaryAttribute(2));

Problem Assembly

The Problem class represents a complete variational problem:

Problem problem(u, v);
problem = Integral(Grad(u), Grad(v))  // Bilinear form
        - Integral(f, v)               // Linear form
        + DirichletBC(u, Zero());      // Boundary conditions

What happens internally:**

1. Rodin assembles the stiffness matrix from Integral(Grad(u), Grad(v))
2. Assembles the load vector from Integral(f, v)
3. Applies boundary conditions by modifying the system
4. Prepares the linear system

Solving

After defining a problem, solve it with an appropriate solver from the Solver module. The most common pattern constructs the solver with the problem and calls solve() inline:

// Inline: construct solver with problem, solve immediately
Solver::CG(problem).solve();

// Or: direct solver
Solver::SparseLU(problem).solve();

For finer control, construct the solver, configure it, then solve:

Solver::CG solver(problem);
solver.setMaxIterations(1000);
solver.setTolerance(1e-10);
solver.solve();

if (solver.success())
  std::cout << "Converged." << std::endl;

You can also pass a standalone solver to problem.solve():

Solver::SparseLU lu;
problem.solve(lu);

After solving, access the solution as a GridFunction:

auto& solution = u.getSolution();

The solution is stored in the GridFunction associated with the trial function. You can evaluate it at any point, compute norms, or export it for visualization.

Creating Grid Functions

P1 Vh(mesh);

// Create a grid function in the space (coefficients initialized to zero)
GridFunction u(Vh);

// Project a lambda expression onto the space
u = [](const Geometry::Point& p) { return p.x() * p.y(); };

// Project a RealFunction
RealFunction f = 1;
GridFunction v(Vh);
v.project(f);

Accessing Values

Individual degrees of freedom can be read or written using operator[]:

GridFunction u(Vh);
u[0] = 1.0;                  // Set DOF 0
Real val = u[42];            // Read DOF 42

The full coefficient vector is available via getData():

auto& coeffs = u.getData();  // Returns Math::Vector<Real>&

Evaluating at Points

A grid function can be evaluated at any Point on the mesh. The getValue() method performs finite element interpolation using the basis functions of the element containing the point:

GridFunction u(Vh);
u = [](const Geometry::Point& p) { return std::sin(p.x()); };

// Evaluate at a quadrature point during assembly (inside integrals)
// or explicitly via getValue():
// Real val = u.getValue(point);

For vector-valued grid functions, component extraction is available:

P1 Vh(mesh, 2);              // 2D vector P1 space
GridFunction u(Vh);
auto ux = u.x();             // Component 0 as scalar function
auto uy = u.y();             // Component 1 as scalar function

Saving and Loading

Grid functions can be saved and loaded in MFEM format:

u.save("solution.gf", IO::FileFormat::MFEM);
u.load("solution.gf", IO::FileFormat::MFEM);

Naming Grid Functions

Assigning a name is required when adding grid functions to XDMF output:

GridFunction u(Vh);
u.setName("velocity");

IO::XDMF xdmf("Output");
auto grid = xdmf.grid();
grid.setMesh(mesh);
grid.add(u);                 // Uses the name "velocity"
xdmf.write();

RealFunction

RealFunction represents a scalar-valued function . It can be constructed from a constant, an integer, or a callable:

// Constant function
RealFunction f = 1;
RealFunction g = 3.14;

// Lambda function of the spatial coordinates
RealFunction h([](const Geometry::Point& p) {
  return std::sin(M_PI * p.x()) * std::cos(M_PI * p.y());
});

VectorFunction

VectorFunction represents a vector-valued function . Components can be constants, RealFunctions, or mixed:

// Constant vector
VectorFunction gravity{0, -9.81};

// Mixed: combine constants with symbolic expressions
VectorFunction load{0, RealFunction([](const Geometry::Point& p) {
  return -std::sin(p.x());
})};

Zero

Zero is a class (not a free function) that produces the zero function. It supports both scalar and vector forms:

// Scalar zero (CTAD: Zero() deduces Zero<Real>)
DirichletBC(u, Zero())

// Vector zero: Zero(dimension)
DirichletBC(u, Zero(3))     // Zero vector in 3D

Mathematical Functions

Rodin provides pointwise mathematical function operators that compose with any scalar function expression:

| Absolute value | Abs(f) | | | Maximum | Max(f, g) | | | Minimum | Min(f, g) | |

RealFunction u_exact([](const Geometry::Point& p) {
  return std::sin(M_PI * p.x());
});
// Or equivalently, composing Rodin operators:
// auto u_exact = Sin(RealFunction([](const Geometry::Point& p) { ... }));

Coordinate Functions

The F namespace provides built-in coordinate projection functions that extract spatial coordinates at evaluation points. These are particularly convenient for defining spatially varying source terms and exact solutions:

#include <Rodin/Variational/F.h>

// Define a source term using coordinate functions
auto f = cos(2 * pi() * F::x + pi() / 2) * sin(2 * pi() * F::y);

// Boolean operations using coordinate functions
GridFunction gf(Vh);
gf = F::x > 0.5;  // 1 where x > 0.5, 0 elsewhere

Complex Functions

For complex-valued problems (e.g. Helmholtz, electromagnetics):

P1<Complex> Vh(mesh);
GridFunction gf(Vh);
gf = [](const Geometry::Point& p) { return Complex(p.x(), p.y()); };

// Extract real and imaginary parts
P1<Real> realVh(mesh);
GridFunction re(realVh), im(realVh);
re = Re(gf);
im = Im(gf);

Jump and Average

The Jump and Average operators act on inter-element faces. For a function with values and on the two sides of an internal face:

Face and Interface Integrals

In addition to Integral (cell interior) and BoundaryIntegral (boundary faces), Rodin provides:

• FaceIntegral — integrates over all** faces (interior and boundary)
• InterfaceIntegral — integrates over interior faces only (excludes boundary)

These are essential for DG penalty and consistency terms:

// DG Poisson example (interior penalty method)
P1 Vh(mesh);
TrialFunction u(Vh);
TestFunction  v(Vh);

Real sigma = 10.0;  // penalty parameter

Problem poisson(u, v);
poisson = 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());

Flow Operator

The Flow operator computes the composition , where is the flow map generated by a velocity field :

Flow traces characteristics backward in time and evaluates the function at the departure point. It is the core mechanism behind the semi-Lagrangian advection in the Advection module. The tracing uses robust cell-by-cell intersection with adaptive substepping and configurable boundary policies.

VectorFunction velocity{1.0, 0.0};
auto departure = Flow(-dt, phi, velocity);  // backward trace

Potential Operator

The Potential operator computes integral convolutions of the form:

where is a kernel function. This is fundamental for boundary element methods (BEM) and integral equation formulations.

// Newtonian potential (3D Laplace kernel)
auto K = [](const Geometry::Point& x, const Geometry::Point& y) -> Real {
  return 1.0 / (4.0 * M_PI * (x.getPhysicalCoordinates() - y.getPhysicalCoordinates()).norm());
};

Problem eq(u, v);
eq = Integral(0.0001 * Grad(u), Grad(v))
   + Integral(Potential(K, u), v)
   - Integral(v);

TraceOperator

The TraceOperator restricts a function to a specified mesh boundary: . For , the trace belongs to .

auto u_boundary = TraceOperator(u, boundaryAttribute);