Flux reconstruction and error estimation

Moment reconstruction, conforming potentials and dimension-dependent estimator hypotheses.

The moment construction accepts represented skeletal normal traces and continuous local tests contained in the primal space. Its algebraic conservation identity is distinct from an approximation theorem. Barrenechea et al. (2026), Theorem 4.7 and Remark 5.1 require \(k\ge\ell+d\), \(\ell\le m\le k\) for \(d=2,3\), compatible shape-regular partitions, the stated solution/source Sobolev regularity, and diffusion polynomial on each fine simplex. In particular the minimum local degrees for constant/linear skeletal traces are P2/P3 in 2D and P3/P4 in 3D. Material callbacks do not certify the theorem's regularity. The same publication states the separate reliability and efficiency hypotheses. Fine-cell balance is an independent implementation diagnostic.

All API families

Locally equilibrated H(div) fluxes from a primal Darcy solution.

EquilibratedFlux dataclass

RT0 fluxes, preserving macro traces and balancing each fine source moment.

conservation_residuals() -> tuple[FloatArray, ...]

Return cellwise integrated divergence minus the assembled source.

l2_error(exact: object, order: int = 5) -> float

Integrate the vector flux error with independent quadrature.

equilibrate_flux(solution: DarcySolution, *, solver: str = 'scipy') -> EquilibratedFlux

Minimize weighted RT0 distance to the raw flux subject to exact balances.

On each macrocell, minimize ||q_h + K grad(p_h)||_(K^-1) subject to div(q_h)|_t = average_t(f) and the prescribed skeleton normal flux. This constrained minimization is distinct from the published face-moment reconstruction; no equivalence to that particular operator is claimed. Degree-zero skeleton segments must coincide with fine boundary faces.

Canonical RT moment reconstruction of primal MHM Darcy fluxes.

This is the face/volume construction of Barrenechea et al. (2026), equations (4.9) and (5.1). It is distinct from an energy-minimizing equilibrated flux: its divergence agrees with the source against continuous local test functions, not necessarily against individual discontinuous fine-cell constants.

MomentFluxSolution dataclass

An H(div) RT flux with explicit continuous-test conservation diagnostics.

flux_l2_error(exact: Any, order: int = 6) -> float

Integrate the physical vector-flux error without post-smoothing.

divergence_l2_error(exact: Any, order: int = 6) -> float

Integrate raw divergence error, without projecting onto continuous tests.

continuous_moment_residuals() -> tuple[FloatArray, ...]

Return moments of div(q)-f against C0 P_m on each macrocell.

For m=0 this test space contains one macroconstant, not separate fine constants. The identity requires that the primal local test space contain this continuous space and that its variational equations hold.

fine_conservation_residuals() -> tuple[FloatArray, ...]

Measure fine-cell integral defects; this reconstruction does not impose them.

conservation_residuals() -> FloatArray

Return macrocell integrated flux-minus-source defects.

normal_flux_residuals() -> tuple[FloatArray, ...]

Compare every RT boundary normal moment with the skeletal multiplier.

continuous_divergence_projection() -> tuple[FloatArray, ...]

Compute the macro-local L2 projection of divergence onto continuous P_m.

projected_divergence_l2_error(exact: Any, order: int = 6) -> float

Integrate error of the continuous projection, distinct from raw divergence.

reconstruct_flux_moments(skeleton: SkeletonSpace, trace: FloatArray, local_meshes: Sequence[TriangleMesh], raw_fluxes: Sequence[RawFlux], *, degree: int = 1, source: Any = 0.0, quadrature_order: int = 6, material: Any = None, backend: Literal['serial', 'thread', 'process'] = 'serial', workers: int | None = None) -> MomentFluxSolution

Reconstruct RT0/RT1/RT2 from signed skeleton and one-sided primal fluxes.

Each macrocell supplies raw_flux(points, fine_cell_indices) returning physical vectors with shape (points,2). This contract makes traces at discontinuous material interfaces explicit. Boundary normal moments come from the skeleton; interior normal moments use the arithmetic average of both one-sided fluxes; volume moments match the original primal flux. A CartesianCellField supplied through material partitions interior faces and volumes at its pixel interfaces before integrating those moments.

The source is used only in diagnostics. No fine-cell equilibrium equation or optimization is added. Trace degree must be <= RT degree and every skeleton break must align with fine boundary edges. Process execution requires picklable raw-flux evaluators; native primal evaluators support this contract.

reconstruct_darcy_moments(solution: DarcySolution, *, degree: int = 1, quadrature_order: int | None = None, raw_fluxes: Sequence[RawFlux] | None = None, backend: Literal['serial', 'thread', 'process'] = 'serial', workers: int | None = None) -> MomentFluxSolution

Apply the RT moment reconstruction to a native primal Pk Darcy solution.

Require m<=k so continuous P_m tests belong to the primal space. This identity alone is weaker than the paper's optimal-estimate hypothesis k>=ell+2 in two dimensions, where ell is the skeleton degree. Theorem 4.7 additionally assumes diffusion polynomial on every fine triangle, piecewise H^(k+1) pressure and A grad(p) in piecewise H^(ell+1) intersected with H(div), with compatible shape-regular partitions. These approximation hypotheses are not inferred from a material callback or algebraic conservation.

CartesianCellField interfaces use exact geometric cuts and one-sided pixel values from the incident triangle. For other callbacks, K is evaluated at physical face points. If K has different one-sided values there, provide raw_fluxes with explicit fine-cell indices; the routine never infers a material limit by perturbing coordinates. Discontinuities solely at macrofaces use the skeletal boundary data and do not require such an override.

Three-dimensional canonical RT reconstruction of primal MHM Darcy flux.

Boundary moments are inherited from the physical multiplier, interior face moments from the arithmetic one-sided raw flux average, and volume moments from the raw flux. The resulting equilibrium is against continuous local Pm tests; it does not impose separate fine-cell source balances.

MomentFlux3DSolution dataclass

Canonical tetrahedral RT flux with continuous-test equilibrium diagnostics.

evaluate(cell: int, points: FloatArray) -> tuple[FloatArray, FloatArray]

Return physical flux and divergence at common reference Cartesian points.

flux_l2_error(exact: Any, order: int = 6) -> float

Integrate the physical vector error with independent positive volume quadrature.

normal_flux_residuals() -> tuple[FloatArray, ...]

Compare every boundary RT moment with its original physical skeletal density.

continuous_moment_residuals() -> tuple[FloatArray, ...]

Return moments of div(q)-f against continuous Pm, or one macroconstant at m=0.

fine_conservation_residuals() -> tuple[FloatArray, ...]

Measure fine-cell integral defects, which canonical reconstruction does not impose.

reconstruct_darcy_moments_3d(solution: Darcy3DSolution, *, degree: int | None = None, quadrature_order: int | None = None, backend: Literal['serial', 'thread', 'process'] = 'serial', workers: int | None = None) -> MomentFlux3DSolution

Construct canonical RT_m face/volume moments of a primal tetrahedral MHM field.

The scalar local test space must contain continuous Pm and the multiplier must belong to its RT normal space. Coefficient quadrature must resolve material interfaces; pointwise callbacks do not supply geometric cuts. Both discontinuous and C0 multiplier bases are supported through their declared subtriangle coefficient maps. This algebraic construction alone does not establish the hypotheses of a posteriori estimator theorems. Theorem 4.7 of Barrenechea et al. (2026) requires k>=ell+3 in three dimensions, ell<=m<=k, diffusion polynomial on every fine tetrahedron, piecewise H^(k+1) pressure and A grad(p) in piecewise H^(ell+1) intersected with H(div), and the paper's compatible shape-regular partitions. Callbacks do not certify those approximation hypotheses.

Conforming potential recovery and the unit-diffusion MHM energy estimator.

ConformingPotential dataclass

Oswald nodal average in a global continuous Pk space with zero boundary.

l2_error(exact: Any, order: int = 8) -> float

Integrate the recovered potential error on the global fine mesh.

DarcyEstimator dataclass

Four macro-local estimator contributions and their recovered fields.

flux_defect, nonconformity, divergence_defect and oscillation are eta_1, eta_2, eta_3 and eta_osc in equations (5.3)--(5.7), specialized to unit diffusion. Their quadrature evaluation is not interval-certified.

local_squared: FloatArray property

Return additive squared macro indicators in the reliability bound.

total: float property

Return the square root of the sum of local squared indicators.

energy_error(exact_gradient: Any, order: int = 10) -> float

Integrate the genuine broken energy error for unit diffusion.

recover_potential(solution: DarcySolution, *, homogeneous_dirichlet: bool) -> ConformingPotential

Average each nodal trace over incident fine triangles, then impose zero BC.

This is the Oswald operator in equation (5.2) of Barrenechea et al. (2026), for globally conforming triangular submeshes and uniform local polynomial degree. Repeated values inside one macrocell count once per incident fine triangle, not once per macrocell. No smoothing is applied to the original solution. The caller must declare homogeneous Dirichlet data because DarcySolution does not retain its boundary specification.

estimate_darcy_error(solution: DarcySolution, *, homogeneous_dirichlet: bool, degree: int = 1, quadrature_order: int = 8) -> DarcyEstimator

Evaluate the energy estimator for identity diffusion and zero Dirichlet data.

Use the energy decomposition of Barrenechea et al. (2026).

Require a primal solution, globally conforming fine triangles, local degree k>=ell+2 and ell<=m<=k. Here ell is the largest skeletal degree and m is the nonnegative RT reconstruction degree. Macrocells are convex triangles, hence their Poincare constant is bounded by diameter/pi.

Coefficients must be the literal scalar 1 or identity matrix; arbitrary SPD tensors and callbacks are rejected rather than inheriting an unsupported coefficient-independent bound. Boundary moments and the continuous-test conservation identity are checked using numerical tolerances. The true source is integrated for oscillation; it is not replaced by a projected source. The theorem assumes exact integration; this implementation reports quadrature-based indicators and a separate equilibrium defect, not an interval-certified upper bound. The theorem requires independent polynomials on each skeletal subface; C0 over several segments of one macroface is excluded, while C0 on a single segment has the same polynomial space as DG.

Energy-normalized MHM flux estimators for SPD materials and mixed boundaries.

The decomposition is that of Barrenechea et al. (2026), with dual norms taken in the physical energy inner product. Material-dependent Poincare weights therefore contain the certified lower ellipticity bound. Boundary traces must be representable; no unmeasured boundary residual is silently omitted.

WeightedDarcyEstimator dataclass

Bases: DarcyEstimator

Four material-weighted terms and their explicit reliability assumptions.

The reported total bounds the broken energy error under exact integration, exact represented boundary data, continuous-test flux equilibrium and the supplied ellipticity lower bounds. Evaluated quadrature is not interval arithmetic.

convention: str property

Identify the energy-normalized material convention explicitly.

energy_error(exact_gradient: Any, order: int = 10) -> float

Integrate the broken physical energy error, resolving declared material cuts.

PublishedDarcyIndicator dataclass

Bases: WeightedDarcyEstimator

Literal equations (5.3)--(5.7), without a general-material energy-bound claim.

The flux residual uses the ordinary L2 norm and the divergence/oscillation terms use diameter/pi, without an ellipticity factor. These terms coincide with the energy-normalized estimator for identity diffusion. For general SPD material they are the published numerical indicators; the physical energy_error remains available but no coefficient-independent reliability guarantee is attached to total. ellipticity_lower_bounds is empty because the printed terms use no ellipticity certificate.

convention: str property

Identify the literal printed material convention.

local_squared: FloatArray property

Return equation (5.6)'s additive squared numerical indicators.

recover_dirichlet_potential(solution: DarcySolution, *, dirichlet: Any = 0.0, neumann: dict[int, Any] | None = None) -> ConformingPotential

Oswald-average the broken field and impose only its Dirichlet boundary nodes.

The boundary function must belong to the continuous local Pk trace space. Representability is checked on every fine Dirichlet edge. Nonpolynomial boundary approximation requires an additional lifting estimator and is rejected here. Neumann nodes remain unconstrained except at Dirichlet corners.

estimate_weighted_darcy_error(solution: DarcySolution, *, degree: int = 1, dirichlet: Any = 0.0, neumann: dict[int, Any] | None = None, ellipticity_lower_bound: Any = None, quadrature_order: int = 8, backend: Literal['serial', 'thread', 'process'] = 'serial', workers: int | None = None) -> WeightedDarcyEstimator

Evaluate the energy-normalized decomposition for SPD diffusion.

The decomposition follows Barrenechea et al. (2026).

Macro cells must be convex triangles; fine meshes must be globally conforming. Require k>=ell+2 and ell<=m<=k. The constant/Cartesian coefficient's minimum eigenvalue provides a certified bound automatically. Other callbacks require an explicit positive lower bound, checked at quadrature nodes in addition to the caller's certification between those nodes. Neumann data are outward physical fluxes and must be exactly represented by the skeletal trace.

estimate_darcy_indicator(solution: DarcySolution, *, convention: Literal['published', 'energy'] = 'published', degree: int = 1, dirichlet: Any = 0.0, neumann: dict[int, Any] | None = None, ellipticity_lower_bound: Any = None, quadrature_order: int = 8, backend: Literal['serial', 'thread', 'process'] = 'serial', workers: int | None = None) -> WeightedDarcyEstimator

Select published indicators or the physical-energy normalization.

The published indicators follow Barrenechea et al. (2026).

published evaluates eta1=||K grad(p)+q||, eta2 in the K energy norm, and eta3/oscillation with diameter/pi, exactly as (5.3)--(5.7) are printed. energy uses K^-1 in eta1 and diameter/(pi sqrt(alpha_K)) in eta3 and oscillation. The latter requires certified ellipticity lower bounds for callable coefficients. Both conventions retain represented boundary data, conforming local partitions, space restrictions and continuous equilibrium checks. The published convention is a numerical indicator for general SPD material, not a claim of a coefficient-independent physical-energy bound. backend executes independent reconstruction and indicator calculations; spawn workers require picklable coefficient and source callbacks. Independent polynomial tests on each skeletal subface are required by the cited theorem. Subdivided C0 macrofaces are excluded in either convention; one-segment C0 and DG spaces coincide.

Tetrahedral Oswald recovery and the MHM flux-based energy decomposition.

The four terms use canonical RT moments, continuous local source/divergence projections and diameter/pi Poincare bounds on convex macrotetrahedra. Energy normalization includes certified lower ellipticity bounds; literal published normalization is an indicator for general materials, without that bound claim.

ConformingPotential3D dataclass

Oswald potential and its restrictions to each original local nodal space.

l2_error(exact: Any, order: int = 6) -> float

Integrate the global conforming potential error with positive quadrature.

Darcy3DEstimator dataclass

Four physical terms with explicit equilibrium and material-bound diagnostics.

local_squared: FloatArray property

Return additive local squared indicators with the prescribed cross terms.

total: float property

Return the root of the total squared indicator.

energy_error(exact_gradient: Any, order: int = 8) -> float

Integrate the actual broken energy error using the physical SPD tensor.

recover_potential_3d(solution: Darcy3DSolution, *, dirichlet: Any = 0.0, neumann: dict[int, Any] | None = None) -> ConformingPotential3D

Average incident fine tetrahedra, then impose exactly represented Dirichlet data.

Local fine partitions must form a globally conforming tetrahedral mesh. Neumann faces remain free except at shared Dirichlet nodes. Nonpolynomial boundary data need a separate lifting estimator and are rejected here.

estimate_darcy_error_3d(solution: Darcy3DSolution, *, degree: int = 1, dirichlet: Any = 0.0, neumann: dict[int, Any] | None = None, ellipticity_lower_bound: Any = None, convention: Literal['energy', 'published'] = 'energy', quadrature_order: int = 7) -> Darcy3DEstimator

Evaluate the dimension-independent energy decomposition on tetrahedra.

Use the estimator of Barrenechea et al. (2026).

Require local k>=ell+3 and ell<=m<=k for skeletal degree ell and RT order m. The degree offset is the spatial dimension in Theorem 5.2, not a dimension-independent value of two. The energy bound assumes exact integration, continuous-test equilibrium, represented boundary data, and certified ellipticity. Numerical quadrature supplies measured indicators, not interval-certified upper bounds.

The cited theorem requires independent polynomial tests on each skeletal subface. A C0 macroface with several subtriangles is therefore excluded; C0 with a single subtriangle has the same polynomial space as DG. Canonical RT reconstruction remains available for subdivided C0 traces without transferring this estimator's theorem to that space.

References

• Gabriel R. Barrenechea, Larissa Martins, Weslley Pereira, and Frédéric Valentin (2026). An H(div; Ω)-Conforming Flux Reconstruction for the Multiscale Hybrid-Mixed Method, Multiscale Modeling & Simulation 24(2), 399–428. DOI: 10.1137/24M1673073.