General concepts » Finite element spaces

Introduction

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

Every finite element computation in Rodin begins by choosing a mesh and a finite element space defined on that mesh:

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

P1 Vh(mesh);  // Creates a P1 finite element space

Available Spaces

Rodin provides several standard finite element spaces:

P0: Piecewise Constant Space

The P0 space consists of functions that are constant on each cell of the mesh. Functions in P0 are generally discontinuous across cell boundaries.

Properties:

• One DOF per cell (the constant value on that cell)
• Number of DOFs = number of cells (retrieved by Wh.getSize())
• Discontinuous across cell boundaries — this is an space, not -conforming
• Used for piecewise constant approximations such as material properties, indicator functions, and cell-based data
• The Grad operator applied to a P0 function is identically zero within each cell

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

P0 Wh(mesh);
std::cout << "P0 DOFs: " << Wh.getSize() << std::endl;
// For an 8×8 grid of triangles: 128 DOFs (= 128 cells)

// Assign a function: each cell stores the function value at its centroid
GridFunction indicator(Wh);
indicator = [](const Geometry::Point& p) {
  return (p.x() < 0.5) ? 1.0 : 0.0;
};

P0g: Global Constant Space

The P0g space is a global constant finite element space — it consists of functions that are constant over the entire mesh (not per-cell like P0).

Properties:

• Scalar P0g: 1 DOF total (the single global constant)
• Vector P0g: vdim DOFs total (one per component)
• Continuous everywhere (trivially, since the function is constant)

P0g Wh(mesh);           // scalar: 1 DOF
P0g Vh(mesh, d);        // vector: d DOFs

P0g is useful for representing globally constant fields (e.g., a uniform property or a Lagrange multiplier for a global constraint).

P1: Piecewise Linear Space

The P1 space is the most commonly used finite element space. It consists of continuous functions that are linear on each cell.

Properties:

• One DOF per vertex (the function value at that vertex)
• Number of DOFs = number of vertices (retrieved by Vh.getSize())
• Continuous across cell boundaries ( -conforming): values at shared vertices are identified, ensuring global continuity
• The gradient Grad(u) is constant within each cell (since a linear function has constant derivatives)
• Approximation error for smooth solutions:
  • norm:
  • seminorm:

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

P1 Vh(mesh);
std::cout << "P1 DOFs: " << Vh.getSize() << std::endl;
// For an 8×8 grid: 81 DOFs (= 9×9 = 81 vertices)

// P1 is the standard choice for most scalar PDEs
TrialFunction u(Vh);
TestFunction  v(Vh);

H1: Higher-Order Space

The H1 space provides -conforming finite elements of arbitrary polynomial order . The polynomial order is specified as a compile-time constant via std::integral_constant:

#include <Rodin/Variational/H1.h>

// Quadratic scalar space (order 2)
H1 Vh(std::integral_constant<size_t, 2>{}, mesh);

// Quadratic vector space (order 2, d components)
H1 Uh(std::integral_constant<size_t, 2>{}, mesh, mesh.getSpaceDimension());

The number of DOFs per element depends on the polynomial degree and the element type:

Convergence rates for smooth solutions:

• norm:
• seminorm:

A common application is the Taylor-Hood pair for incompressible flow, which uses different polynomial orders for velocity and pressure to satisfy the inf-sup stability condition:

• P2 velocity (order 2): H1 uh(std::integral_constant<size_t, 2>{}, mesh, d)
• P1 pressure (order 1): H1 ph(std::integral_constant<size_t, 1>{}, mesh)

Vector-Valued Spaces

P0, P1, and H1 all support vector-valued function spaces by specifying the number of components. This is essential for problems involving vector fields such as displacement (elasticity) or velocity (fluid dynamics).

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

// Scalar space (default, 1 component)
P1 scalarVh(mesh);

// 2D vector space (2 components)
P1 vectorVh(mesh, 2);

// 3D vector space (3 components)
P1 vector3Vh(mesh, 3);

For the elasticity equation, a vector-valued P1 space represents the displacement field :

size_t d = mesh.getSpaceDimension();
P1 Vh(mesh, d);

TrialFunction u(Vh);  // u is a vector function
TestFunction  v(Vh);  // v is a vector function

Complex-Valued Spaces

Rodin supports complex-valued finite element spaces using C++ complex types. This is useful for problems in wave propagation, acoustics, and electromagnetics:

// Complex-valued P1 space
P1<Complex> Vh(mesh);

GridFunction gf(Vh);
gf = [](const Geometry::Point& p) {
  return Complex(p.x(), p.y());  // real + imaginary
};

// Extract real and imaginary parts for output
P1<Real> realVh(mesh);
GridFunction gfRe(realVh), gfIm(realVh);
gfRe = Re(gf);
gfIm = Im(gf);

Grid Functions

A GridFunction is a concrete function living in a finite element space. It stores the coefficients (DOF values) and can be evaluated, saved, and loaded.

P1 Vh(mesh);

// Create a grid function in the P1 space
GridFunction gf(Vh);

// Interpolate a smooth function into the FE space
gf = [](const Geometry::Point& p) {
  return std::sin(M_PI * p.x()) * std::sin(M_PI * p.y());
};

// Assign from a constant
GridFunction constant(Vh);
constant = 1.0;

P0 Wh(mesh);
GridFunction indicator(Wh);
indicator = [](const Point& p) {
  return p.x() * p.x() + p.y() * p.y();
};

P1 scalarVh(mesh);
GridFunction temperature(scalarVh);
temperature = [](const Point& p) {
  return 1 - (p.x() - 0.5) * (p.x() - 0.5)
           - (p.y() - 0.5) * (p.y() - 0.5);
};

// Interpolate the gradient into a vector P1 space
size_t d = mesh.getSpaceDimension();
P1 vectorVh(mesh, d);

GridFunction grad(vectorVh);
grad = Grad(temperature);

GridFunction trace(vectorVh);
trace.project(Region::Faces, Grad(temperature).traceOf(1));

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

Problem poisson(u, v);
// ... define problem ...
Solver::CG(poisson).solve();

// Access the solution grid function
auto& solution = u.getSolution();

P1 Vh(mesh);

auto f = cos(2 * pi() * F::x + pi() / 2) * sin(2 * pi() * F::y);
auto g = sin(2 * pi() * F::x) * sin(2 * pi() * F::y);

GridFunction gf(Vh);
gf = f > g;         // 1 where f > g, 0 otherwise
gf = f < g;         // 1 where f < g, 0 otherwise
gf = (f > g) || (f < g);  // Logical OR
gf = (f > g) && (f < g);  // Logical AND (always 0 here)

// Save
gf.save("temperature.gf", IO::FileFormat::MFEM);

// Load
GridFunction loaded(Vh);
loaded.load("temperature.gf", IO::FileFormat::MFEM);

IO::XDMF xdmf("output");
xdmf.grid().setMesh(mesh).add("temperature", gf);
xdmf.write();

Matrix-valued solution spaces

The local and MPI variants of P0, P0g, P1, and H1 support SpatialMatrix as their range. Specify the row and column counts independently; each must lie between one and three.

P0 cellStress(mesh, 3, 3);
P0g constantStress(mesh, 3, 3);
P1 stressSpace(mesh, 3, 3);
H1 higherOrderStress(std::integral_constant<size_t, 3>{}, mesh, 3, 3);

GridFunction stress(stressSpace);
Math::SpatialMatrix<Real> initial(3, 3);
initial.setZero();
stress = initial;

TrialFunction sigma(stressSpace);
TestFunction tau(stressSpace);
auto mass = Integral(sigma, tau);
auto entryMass = Integral(Component(sigma, 0, 1), Component(tau, 0, 1));

Each scalar basis function is replicated for every matrix entry. For a range with rows and columns, the basis is , the component index is , and the global DOF index is . No symmetry or Voigt scaling is imposed. Matrix contractions use the Frobenius product (conjugating the second operand for complex entries). The scalar solver coefficients remain flat, while basis and grid-function evaluations return matrices of the specified shape. Transpose and matrix-component expressions can be used in forms.

All four families support real and complex matrix entries. Matrix gradients and Jacobians return rank-three SpatialTensor values with the derivative axis last. Matrix divergence contracts the column and derivative indices, and Derivative selects one physical derivative direction. See Matrix solutions and spatial tensors for tensor coefficients and backend details.

Choosing the Right Space

Rule of thumb: Start with P1 for scalar problems. Use vector-valued P1 for elasticity and fluid dynamics. Move to higher-order H1 when you need better accuracy per DOF or when the PDE's solution is smooth enough to benefit from higher polynomial approximation.

See Also

• Ciarlet's definition — Mathematical foundations
• Core Concepts — Practical guide
• MPI distributed computing — Distributed P1 on Mesh<Context::MPI> with local↔global DOF index maps
• PETSc integration — PETSc-backed grid functions with ghost-layer synchronization
• Poisson example — Using P1 in practice
• Elasticity example — Vector-valued spaces