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 cell-average value)
• Number of DOFs = number of cells
• Discontinuous across cell boundaries
• Used for piecewise constant approximations (e.g. material properties, indicator functions)

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

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

GridFunction indicator(Wh);
indicator = [](const Geometry::Point& p) {
  return (p.x() < 0.5) ? 1.0 : 0.0;
};

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 the vertex)
• Number of DOFs = number of vertices
• Continuous across cell boundaries ( -conforming)
• Approximation error: where is the mesh size

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

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

// P1 is the default 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. This is useful for problems requiring higher accuracy or mixed formulations (e.g. Taylor-Hood elements for the Navier-Stokes equations).

The polynomial order is specified as a compile-time 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());

A common application is the Taylor-Hood pair for incompressible flow:

• 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

See Also

• Ciarlet's definition — Mathematical foundations
• Core Concepts — Practical guide
• Poisson example — Using P1 in practice
• Elasticity example — Vector-valued spaces