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 defines the type of approximation used. Rodin provides several standard spaces:

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

// Linear elements
P1 Vh1(mesh);
std::cout << "P1 space DOFs: " << Vh1.getSize() << std::endl;

// Quadratic elements (if available)
// P2 Vh2(mesh);
// std::cout << "P2 space DOFs: " << Vh2.getSize() << std::endl;

Trial and Test Functions

In the finite element method, we express problems in weak form using trial and test functions.

P1 Vh(mesh);

// Trial function: the unknown we're solving for
TrialFunction u(Vh);

// Test function: used in weak formulation
TestFunction v(Vh);

Mathematical background:**

• We seek (the trial space)
• We test against all (the test space)
• The discrete problem is: Find such that for all

Differential Operators

Rodin provides operators that work on trial and test functions:

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

// Gradient operator
auto grad_u = Grad(u);
auto grad_v = Grad(v);

// Used in integrals
auto bilinear_form = Integral(Grad(u), Grad(v));

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:

// Direct solver (good for small to medium problems)
SparseLU().solve(problem);

// Iterative solver (good for large problems)
CG solver;
solver.setMaxIterations(1000);
solver.setRelativeTolerance(1e-10);
solver.solve(problem);

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