General concepts » Numerical contracts

Interpolation, projection, assembly signs, quadrature, and solver contracts.

Rodin tries to make code read like mathematics, but some operations have precise numerical contracts. This page collects the rules that affect accuracy and solver behavior across modules.

Interpolation vs Projection

Assigning an expression to a GridFunction applies the finite element space's degree-of-freedom functionals to that expression:

GridFunction uh(Vh);
uh = f;        // interpolate f in the DOF sense of Vh

This is interpolation, not an -orthogonal projection. For nodal P1 spaces, the DOF functionals are vertex evaluations, so the assignment stores nodal values. For higher-order H1<K> spaces, Rodin uses higher-order bases and DOF machinery; coefficients should not be read as nodal values. Evaluate the GridFunction at a Point instead.

L2 Projection

When the mathematical operation is the projection, assemble the mass problem:

In Rodin this is written as:

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

Problem l2(u, v);
l2 = Integral(u, v) - Integral(f, v);
Solver::CG(l2).solve();

const auto& uh = u.getSolution();

Use this route for optimal-order transfer of rough data, conservative projection workflows, or whenever the variational definition of the projection matters.

Problem Sign Convention

A Problem assignment writes the residual equation with all terms on the left-hand side. The Poisson weak form

becomes:

Problem poisson(u, v);
poisson = Integral(Grad(u), Grad(v))
        - Integral(f, v)
        + DirichletBC(u, Zero());

Bilinear terms contribute to the operator . Linear terms contribute to the right-hand side with the sign implied by the residual equation. For nonlinear Newton problems, assemble the tangent and residual so the linearized system is

Do not write a second manual negation around a residual term that the Problem assembly already moves to the right-hand side.

LinearSystem Lifetime

A LinearSystem is tied to the finite element spaces and block layout that created it. If the mesh, finite element space, number of fields, or global size changes, create a new problem or linear system.

PETSc makes this rule explicit: once a Mat or Vec has been laid out and assembled, it cannot be resized in place. Reassembly may reuse storage when sizes and sparsity structure match, but changed structure is a new system.

Quadrature Exactness

Quadrature order controls the mathematical form being assembled. For polynomial integrands on affine elements, choose a rule exact for the product degree. Curved transformations, material coefficients, nonlinear constitutive laws, and flow-map evaluations can raise the effective degree or make the integrand non-polynomial.

If a specialized integrator uses a shortcut, verify it against the generic quadrature path or a manufactured solution. Agreement should include the assembled values, convergence rates, and any invariants the form is expected to satisfy.

Residual and Tangent Consistency

Nonlinear residual/tangent pairs must be consistent. For a residual and tangent , the check is:

Rodin provides finite-difference probe helpers under Rodin::Test. Run these checks when adding or changing nonlinear terms, especially in solid mechanics and mesh optimization.

Constraints

Dirichlet and identification boundary conditions are structural constraints. They are represented in the assembly constraint map and eliminated or expanded algebraically. Use penalty terms only when the model is explicitly a penalty method.