MeanField_SERiF/src/include/serif/backend/mfem/README.md

MFEM backend

MFEM implementations live here and in the matching src/lib/serif/backend/mfem tree. Physics descriptions, equations, EOS relations, normalization rules, Newton acceptance, and continuation scheduling stay in their ordinary SERiF modules. The directory migration preserves physics-oriented namespaces; include paths identify the implementation backend. Helpers whose responsibility changed now follow their owner: model::make_discretization, geometry::surface::AutomaticPressureSurface, generic preconditioner diagonal regularization, and discretization element coordinates.

Directory	Responsibility
discretization	Bind physical terms to MFEM spaces; own coefficient ordering, element orientation, basis evaluation, quadrature, traces, and MFEM mesh extraction.
mapping	Evaluate physical mappings and their variations from supplied reference geometry and current coefficients.
geometry	Adapt surface and deformation descriptions to MFEM resources, including pressure traces and boundary metrics.
operators	Integrate equation contributions and expose their linear actions, transposes, diagonals, and assembly. No EOS-specific pressure construction belongs here.
model	Connect the preceding adapters into the model's computational graph and publish owning evaluated systems. The deformation-to-equation bridge belongs here.
solver	Supply MFEM linear solvers, inverses, geometry admissibility, and stellar adapters to the common Newton and continuation algorithms.
projection	Project reusable radial profiles into MFEM finite-element spaces.
shared	Small backend facilities with multiple independent consumers: live mesh/space snapshots, boundary traversal, vector spans, and matrix/vector checks.

Ownership and reuse

The computational graph owns immutable prepared data, including fixed geometry at equation quadrature points. A mapper does not own a second mutable reference-geometry cache. Standalone point and face mapping evaluate their own reference geometry directly. Accepted systems retain their preparation owners and remain independent of subsequent candidate preparation.

MFEM resources are externally mutable. Snapshot policies distinguish topology, material/boundary labels, finite-element layout, and physical/logical coordinate revisions. Input factories capture current revisions; reusing stale input metadata is rejected. Notify MFEM of direct coordinate changes (for example with NodesUpdated) and rebuild the input context after changing its resources. Revision tracking cannot detect unreported writes to external coordinate storage.

A DofMap owns the meaning of reduced coefficients. Public coefficient access validates that meaning. Internal element traversals obtain a validated Elements view once, then retain element bounds, support, and vector-size checks without rescanning all mesh attributes for every element. Do not mutate the mesh or space during a traversal. Element orientation is implemented once in discretization/element_coordinates.

Surface normalization and surface preconditioning share validated boundary traversal, but keep their different integrands: vector radial/harmonic metric diagonals versus scalar H1 mass and stiffness matrices. Likewise, model block selection preserves native child assembly and repeated ranges; restriction of an arbitrary preconditioner operator retains its separate unique-index contract.

Core boundary

serif_core_lib and serif_core_dep have no MFEM dependency. The core header check and the standard-vector Newton tests compile without MFEM include paths. The public serif library includes the core and adds this backend. Convenience entry headers may select MFEM adapters, while physical configuration and numerical algorithms remain available separately.

Continuation families, physical parameter laws, coordinate stepping, and the driver share serif::solver::continuation in solver/continuation. Its entry header is backend independent; only stellar correction wiring and finite-element state transfer live in backend/mfem/solver/continuation.

The common Newton loop asks its vector adapter only for owned vectors and coefficient spans. The MFEM wrapper supplies linear-solve wiring; it does not duplicate acceptance, backtracking, or checkpoint ownership. A future backend should implement its discretization and execution adapters and reuse these common rules, rather than copy the nonlinear or continuation algorithms.

Current term bindings require serial, conforming three-dimensional MFEM spaces. Supported value, gradient, divergence, trace, and tensor/dense basis paths retain their explicit compile-time and runtime checks. This layout does not imply that every MFEM space or another finite-element library is already supported.

Optional surface density reconstruction

Pass the numerical choice directly to the stellar factory, after any deformation:

namespace reconstruction = serif::discretization::reconstruction;
using Density = serif::discretization::blocks::density::mass;
auto finite_elements = serif::model::make_stellar_discretization(
    star, mesh_options, moving_boundary,
    reconstruction::SurfacePower<Density>{
        .density_exponent = n,
        .region = reconstruction::SurfaceAdjacent<
            serif::discretization::domain::StellarSurfaceBoundary>{}});

Omitting the argument preserves ordinary interpolation. Exponent 1 takes that same path without a reconstruction plan. Stored coefficients remain physical nodal densities. Only material elements touching the selected boundary change: with inward reference coordinate t=0 at the surface and t=1 at their inner face, their zero-trace radial space becomes t^n P_(p-1). Tangential interpolation and deeper elements remain ordinary. The full basis retains surface degrees of freedom and reproduces constants; the vacuum solve enforces their zero values explicitly.

The plan belongs to the field binding and is immutable throughout a solve. Its exponent is an explicit numerical choice; changing an EOS during continuation does not silently change the representation. A different profile, mesh or exponent requires a new binding and deliberate field transfer. Construction validates matching radial coordinates and the actual hanging-node constraints. Tangential/equatorial refinement is supported when those constraints preserve the reconstructed traces. Incompatible radial hanging layers fail with a SERiF diagnostic.

discretization::backend::mfem::FieldValues<Term> reads the same physical representation as the equation operators. Use it for arbitrary-point values, scalar H1 gradients or batches of samples. Its binding, coefficients and mesh resources must outlive the reader. Extra deformation/compactification maps must be supplied for physical gradients and H(div) values. Ordinary MFEM GridFunction interpolation of reconstructed nodal values is not the represented field.

Reference preparation replaces only the radial tensor factors; forward and transpose contractions retain the existing kernels. Automatic volume rules include a fractional-power allowance before quadrature-policy resolution, so an explicit fixed order remains authoritative. Gaussian quadrature is still approximate and should be checked by increasing its order for a new regime. Nonnegativity certifies the polynomial carrier with its positive weight and requires the selected zero-density trace. This implementation does not provide a certificate for finite-pressure surface lifting or a general EOS surface profile.

To add a numerical form, provide an ADL-visible make_reconstruction(specification, binding) overload, a specification with using Term, and an immutable RadialProfile that evaluates complete nodal functions and inward derivatives, including endpoints. The profile must preserve nodal values, constants and compatible face traces. Add its admissibility certificate and archive codec explicitly; missing capabilities report errors. No MFEM element subclass or physics-equation rewrite is required. The current surface-power form supports scalar H1 mass density only.

Weighted archive bindings use schema version 2; older readers reject them rather than reinterpret their coefficients. The live Rayfem adapter sends a higher-order sampled display field labeled density (sampled). This display approximation is never used by the solver or TOP exporter.