formoniq: higher-order FEEC, the trimmed family P⁻_rΛᵏ

[luiswirth]

Motivation

The engine implements exactly one finite element space, the lowest-order trimmed
family P⁻₁Λᵏ, the Whitney forms. That is the space the de Rham map is dual to
and the one the whole commuting-diagram story is cleanest on, but it fixes the
convergence rate at first order in the energy norm regardless of how smooth the
solution is. Refinement is currently the only way to buy accuracy, and it is the
expensive way.

The trimmed family P⁻_rΛᵏ is the generalization that keeps everything the
design rests on: it is a subcomplex of the de Rham complex for every r ≥ 1, it
has a geometric decomposition into face-supported degrees of freedom, and it
carries a bounded commuting projection. Nothing about the exterior-calculus
formulation has to be given up to get higher order, which is precisely why this
is worth doing here rather than bolting on a classical high-order element.

The mathematics

What already fits

The element mass matrix generalizes without a new idea.
WhitneyCoeffs::pullback already factors the integrand into a blade half H
and a barycentric half Q, because the blades of a Whitney form are constant
and the coordinates carry the whole x-dependence. A higher-order basis
function is λ^α W_σ, which changes only the second factor: unit_bary_gramian
(the pairwise ∫λ_v λ_w) becomes the exact monomial integral

∫_K̂ λ^α = n! α! / (|α| + n)!,

indexed by a Composition rather than by a vertex pair. Composition is
already the right object and already indexes the Freudenthal lattice in atlas.

SimplexQuadRule::degree is arbitrary-degree already, so the moment degrees of
freedom and the L² error have their integration. The local DOF set is
unit_subsimps × Composition::all, in colex, with no new combinatorics.

Where the real derivation is

DifElmat currently takes d to be unit_boundary_operator(dim, grade)
transposed, a purely combinatorial matrix, because dW_σ is Whitney again by
the deletion formula. At higher order d: P⁻_rΛᵏ -> P⁻_rΛᵏ⁺¹ still holds — that
is what makes the trimmed family a complex — but the matrix is no longer the
coboundary. It is a constant rectangular matrix computed once on the reference
cell, so it is tractable, but it is a derivation rather than a relabelling, and
it is where a sign or a basis-ordering slip would hide.

The conformity question, and why colex answers it

For r > 1 the face degrees of freedom are moments against a polynomial basis
on the face, and such a basis needs a vertex ordering of that face. Two cells
sharing a face must agree on it or the space is not conforming.

Colex supplies this for free: a face's colex vertex order is the same read from
either cell, since colex is preserved under passing to a subset.
Cochain::trace already relies on exactly that. So assembly does not acquire a
dependence on CellOrdering, and the separation between a mesh's ordering and
its assembly-relevant content survives higher order intact.

The one thing to get right is that the moment basis must be a function of the
face's own colex-ordered vertices alone, never of the containing cell. A basis
defined through the cell would be the point where conformity silently breaks.

The engine

The DOF layer the cochain currently is

A degree of freedom is currently identical to a simplex, and not by
convention. Cochain is one coefficient per k-simplex indexed by
SimplexRef. assemble_galmat obtains the global index as
cell.faces(grade).kidx(), so the local-to-global map is the face
incidence. matfree::ElementOperator is built from two FaceIncidence on the
same identification. bc, prolongate, derham_map, WhitneyInterpolant, the
.vtu export and studio's Scene all consume Cochain.

For r > 1 a degree of freedom is a pair: a subsimplex σ of dimension at
least k, together with a moment index, an element of a polynomial space on
σ. So the identification has to break, and a DOF numbering layer has to sit
between the Complex and the linear algebra: per grade-and-dimension blocks
with offsets, and a gather returning (simplex, multi-index) -> usize.

This layer is the bulk of the work. The mathematics above it is comparatively
tame.

Static condensation, and where it stops being automatic

The geometric decomposition attaches degrees of freedom to every face of the
cell, the cell itself included. Those interior ones are shared with no
neighbour, and their count grows one power of r faster than the interface
ones, so the fraction of the local space that is element-private grows with the
order. Putting them into the global system and then solving for them there is
the wasteful way round: eliminated element-locally, they never enter it. This
is not an optimization to defer, it is the reason high-order methods are
affordable at all.

The elimination is exact. Partition the local indices into interior I and
interface Γ:

[A_II  A_IΓ] [u_I]   [f_I]
[A_ΓI  A_ΓΓ] [u_Γ] = [f_Γ],

and the condensed element matrix is the Schur complement
S = A_ΓΓ − A_ΓI A_II⁻¹ A_IΓ with load f_Γ − A_ΓI A_II⁻¹ f_I. It is block
Gaussian elimination on a matrix that is already local, so it costs one dense
interior solve per cell, introduces no approximation, and the interior values
are recovered per cell after the global solve.

It fits the existing seam without widening it. ElMatProvider::eval returns an
element matrix from a metric and a chart, so condensation is a combinator over
that trait and its ElVecProvider counterpart rather than a change to
assembly. What it needs from the DOF layer is the interior/interface partition,
which is the geometric decomposition, so the two are one piece of work and
not two.

Two things are not automatic, and they are why this is a caveat rather than a
step:

A_II has to be invertible. For a mass matrix it is a principal submatrix
of an SPD matrix and it is. For the mixed Hodge–Laplace system it is not
automatic, and that system is not one bilinear form: HodgeBlocks assembles the
masses M_{k−1}, M_k, M_{k+1} and the coboundaries D^{k−1}, D^k separately
and mixed_hodge_laplacian builds the saddle point from them, so there is no
single element matrix to condense. Condensing the blocks and then forming the
saddle point is a different operation from condensing the assembled system, and
whether the two agree is to be checked rather than assumed.

The elimination must respect the complex. The trimmed family being a
subcomplex is the whole reason to use it, and d carries the local space at
grade k into the one at grade k+1. Eliminating interior degrees of freedom
grade by grade is harmless only if it commutes with d; otherwise the condensed
system is assembled from something that is no longer a complex, and the
cohomology, the harmonic space and the commuting projection all sit downstream
of that. Settle this before writing any of it. Where it does not commute,
condensation is confined to the primal elliptic forms and the mixed formulation
keeps its interior degrees of freedom.

Scope, in order

• The DOF numbering layer, at r = 1 first, so it reproduces the current
  results bit for bit before any new mathematics enters. The existing Whitney
  path is the regression test.
• The generalized element mass matrix over Composition-indexed monomials.
• The local d matrix for P⁻_r on the reference cell.
• The canonical projection with moment degrees of freedom, generalizing
  derham_map.
• The interior/interface partition and the condensing combinator over
  ElMatProvider/ElVecProvider, with per-cell recovery after the solve.
• Downstream: bc, prolongate, realize, studio, the multigrid transfers.
  Broad and mechanical, but it is where the elapsed time goes.

Laws / validation

• d ∘ d = 0 on the local bases, at every r, grade and dimension.
• The commuting property d ∘ Π = Π ∘ d for the canonical projection, which is
  what makes the family a subcomplex rather than merely a bigger space.
• Unisolvence: the degrees of freedom on the reference cell are a basis, so the
  local generalized Vandermonde matrix is invertible.
• r = 1 reproduces the Whitney results exactly.
• The condensed system reproduces the interface part of the uncondensed
  solution, and recovery reproduces the interior part, to round-off at every
  r. Block elimination is exact, so anything else is a bug in the partition.
• Convergence at rate r in the energy norm on a smooth solution, swept over
  grades — the end-to-end statement, run by hand as an example.
• The Hodge–Laplace spectrum and the Betti numbers are unchanged by raising r,
  since the cohomology is.

Listed under "higher-order (trimmed polynomial) elements" in the roadmap; this
is what that would actually take.