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)

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

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

// Project 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);
};

// Project 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();

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