Adaptive scalar and elasticity indicators

Darcy refinement policies, local error controls, metric remeshing and the primal elasticity indicator.

Darcy and transport controllers accept solve_step to assemble each state from user-defined equations. The callable receives the current mesh and skeleton with the problem's keyword options; it returns the physical solution required by the chosen estimator. The default formulation remains available. The reconstruction's algebraic moment conditions, estimator reliability and optimal approximation estimates have separate hypotheses. In particular, the recovery paper's optimal flux estimate assumes polynomial diffusion on each fine simplex and its stated pressure/flux regularity; exact integration of a material cut does not establish that fine-simplex material condition. See Barrenechea et al. (2026) and the material-weighted estimator for the applicable norms, geometry, boundary representation and equilibrium conditions. See the variational guide for composition with mathematical definition factories.

All API families

Estimator-driven conforming macro refinement for primal Darcy problems.

AdaptiveDarcyResult dataclass

Solved states, estimators, marking decisions and actual refinement maps.

totals: FloatArray property

Return the selected estimator or numerical indicator at every solved level.

mark_dorfler(local_squared: Any, theta: float = 0.5) -> np.ndarray

Select a minimal-cardinality deterministic bulk set of squared indicators.

The selected sum is at least theta times the total sum. Ties follow original cell order. Zero indicators produce an empty set. Scaling before summation avoids overflow and preserves the marking for extreme units.

solve_adaptive_darcy(mesh: TriangleMesh, *, iterations: int = 5, theta: float = 0.5, tolerance: float = 0.0, maximum_cells: int = 10000, trace_degree: int = 0, trace_segments: int = 1, trace_space: FaceSpace | None = None, reconstruction_degree: int = 1, estimator_order: int = 8, estimator_convention: Literal['published', 'energy'] = 'energy', estimator_backend: Literal['serial', 'thread', 'process'] = 'serial', estimator_workers: int | None = None, ellipticity_lower_bound: Any = None, skeleton: SkeletonSpace | None = None, local_mesh_factory: Callable[[TriangleMesh, int], TriangleMesh] | None = None, macro_refiner: Callable[[TriangleMesh, Any], TriangleRefinement] = refine_triangles, solve_step: Callable[..., DarcySolution] | None = None, **problem: Any) -> AdaptiveDarcyResult

Solve, estimate, mark and red-green-refine a primal Darcy macro mesh.

This loop uses the weighted energy decomposition of Barrenechea et al. (2026) and Dörfler bulk marking; it is an original adaptive policy, not a reproduction of a published marking schedule. All estimator assumptions remain required: polynomial represented boundary data, conforming fine triangulations, certified ellipticity and local degree >= trace degree + 2. Local refinement does not certify the recovery paper's optimal flux approximation hypotheses: Theorem 4.7 additionally requires polynomial diffusion on each fine triangle, H^(k+1) pressure and A grad(p) in H^(ell+1) intersected with H(div), with compatible shape-regular partitions. These are separate from the energy decomposition's reliability assumptions. Local refinement and degree stay fixed. Neumann face data and macro-local ellipticity bounds follow exact ancestry. iterations counts solves. A proposed refinement exceeding maximum_cells is not solved. No contraction or optimality theorem is inferred from the loop, and indicator changes need not be monotone. local_mesh_factory(macro_mesh, cell_index) can construct material-fitted local partitions afresh at each level. Static local_meshes are rejected because their ancestry changes after refinement. macro_refiner defaults to red-green refinement. Passing refine_longest_edge selects conforming longest-edge propagation; either operation must provide cell/face ancestors for boundary and trace transfer. estimator_convention='published' selects the literal numerical terms (5.3)--(5.7); the default 'energy' uses physical material weights. The printed convention does not imply a general-SPD energy upper bound. estimator_backend and estimator_workers select independent execution of local reconstruction and indicator terms separately from the PDE solver. trace_space declares the scalar face space used initially and on new macrofaces; its literal degree and continuity override trace_degree and trace_segments. Existing faces retain their transferred partitions. solve_step(mesh, skeleton=current_space, **problem_data) may assemble the user's own mathematical equations at each state. It receives the actual transferred boundary data and local partitions; its returned physical field must satisfy the estimator's existing hypotheses. The default retains the built-in solve. Callback exceptions propagate without changing refinement.

Explicit local/macro refinement balance for the reconstructed Darcy estimator.

BalancedDarcyResult dataclass

Solved states and the declared refinement decisions, including budget termination.

solve_balanced_adaptive_darcy(mesh: TriangleMesh, *, iterations: int = 5, local_refinement: int = 1, maximum_local_refinement: int = 8, local_error_ratio: float = 1.0, theta: float = 0.5, tolerance: float = 0.0, maximum_cells: int = 10000, trace_degree: int = 0, trace_segments: int = 1, reconstruction_degree: int = 2, estimator_order: int = 8, estimator_convention: Literal['published', 'energy'] = 'energy', ellipticity_lower_bound: Any = None, skeleton: SkeletonSpace | None = None, local_mesh_factory: Callable[[TriangleMesh, int, int], TriangleMesh] | None = None, local_error_indicator: Callable[[DarcySolution], float] | None = None, on_state: Callable[[BalancedDarcyResult], None] | None = None, macro_refiner: Callable[[TriangleMesh, Any], TriangleRefinement] = refine_triangles, solve_step: Callable[..., DarcySolution] | None = None, **problem: Any) -> BalancedDarcyResult

Alternate globally uniform local refinement and marked macro refinement.

Compare norm(eta_div + eta_osc) with local_error_ratio times sqrt(norm(eta_flux)**2 + norm(eta_nonconf)**2). When the first exceeds the second, double the local refinement everywhere; otherwise apply Dörfler marking and conforming red-green macro refinement. These are refinement proxies, not an exact separation of local and global errors. No estimator-contraction or complexity-optimality theorem is asserted.

Uniform local refinement keeps shared fine boundary partitions compatible. A custom local_mesh_factory(mesh, cell, refinement) must likewise give a globally conforming partition, checked by potential recovery. It can fit Cartesian material interfaces while leaving the macro mesh unchanged. A local-resolution limit terminates the loop explicitly rather than substituting macro refinement for unresolved local indicators. local_error_indicator(solution) may supply a nonnegative alternative, such as an energy difference between nested local solves with unchanged trace data. Such a difference requires a saturation assumption to bound the actual local error; the policy does not assume that bound. on_state receives the complete accumulated result after each validated solve and refinement decision. It may save checkpoints; its exceptions are propagated. No incomplete solution is passed to the callback. macro_refiner selects the conforming macro refinement operation and its exact ancestry. The default is red-green; refine_longest_edge instead applies longest-edge propagation without green descendants. estimator_convention selects the literal published numerical terms or physical-energy material weights, without changing the finite-element PDE. solve_step(mesh, skeleton=current_space, **problem_data) may assemble the user's own mathematical equations at each state. It receives the actual transferred boundary data and local partitions; its returned physical field must satisfy the estimator's existing hypotheses. The default retains the built-in solve. Callback exceptions propagate without changing refinement.

Indicator-ranked macro refinement with an explicit mesh-complexity budget.

DarcyBudgetRefinement dataclass

A conforming refinement, its indicator prefix and the requested cell count.

The actual count can exceed the target because of conformity propagation or the required bulk set. The construction does not claim minimum complexity among all conforming meshes or use an exact/reference solution for marking.

refine_darcy_budget(mesh: TriangleMesh, local_squared: Any, *, target_cells: int, theta: float = 0.5, refiner: Callable[[TriangleMesh, Any], TriangleRefinement] = refine_longest_edge) -> DarcyBudgetRefinement

Enlarge a Dörfler set along decreasing indicators to meet a cell budget.

Start with the deterministic minimal bulk prefix. If its conforming closure already reaches the target, keep it. Otherwise bracket larger prefixes by the full marked set and search between tested prefixes. The returned set contains the bulk prefix and its actual conforming mesh has at least the target count; conformity can cause overshoot. No monotonicity of arbitrary user-supplied refiners or minimal-prefix optimality is assumed. A target beyond the full set's one-call refinement is rejected; use multiple levels. Exact cell/face ancestry is supplied unchanged by the selected refiner.

Estimator-driven conforming tetrahedral MHM adaptation with physical boundary transfer.

AdaptiveDarcy3DResult dataclass

Solved 3D states, indicators, markings and exact refinement ancestry.

totals: FloatArray property

Return the selected indicator at each completed state.

solve_adaptive_darcy_3d(mesh: TetraMesh, *, iterations: int = 5, theta: float = 0.5, tolerance: float = 0.0, maximum_cells: int = 10000, trace_degree: int = 0, trace_subdivisions: int = 1, reconstruction_degree: int = 1, estimator_order: int = 7, estimator_convention: Literal['energy', 'published'] = 'energy', ellipticity_lower_bound: Any = None, on_state: Callable[[int, Darcy3DSolution, Darcy3DEstimator], None] | None = None, solve_step: Callable[..., Darcy3DSolution] | None = None, **problem: Any) -> AdaptiveDarcy3DResult

Solve, estimate, mark and bisect tetrahedral edge stars.

The marking policy is original, using the dimension-independent estimator of Barrenechea et al. (2026) and Dörfler bulk sets. It is not a reproduction of a historical 3D mesh sequence and implies no contraction or optimality guarantee. Physical Neumann data and per-cell ellipticity bounds follow their exact parent entities. The callback persists each solved state before refinement.

Estimation requires k>=ell+3, ell<=m<=k, globally conforming fine tetrahedra, represented boundary data, continuous-test equilibrium and certified ellipticity. Theorem 4.7's optimal flux approximation additionally requires polynomial diffusion on each fine tetrahedron, H^(k+1) pressure and A grad(p) in H^(ell+1) intersected with H(div), with compatible shape-regular partitions. Those extra approximation hypotheses are not inferred from material callbacks.

solve_step(mesh, skeleton=current, **problem) may assemble each state through user-declared public equations. It must return a Darcy3DSolution in the supplied mesh and trace space; all physical data, coefficient bounds and refinement ancestry retain their declared meaning. None preserves the built-in Darcy solver. Selecting a callable changes only state assembly, not the estimator, marking, boundary transfer or stopping policy.

Face-jump refinement indicators for the one-level Darcy MHM analysis.

DarcyJumpEstimator dataclass

Face indicators with the original cell-side multiplicity.

The face indicators follow Araya et al. (2013).

Interior faces contribute once to each of their two incident cells. The published reliability theorem concerns exact local lifts; this quantity alone does not bound the error of finite-dimensional local solves.

total: float property

Return the square root of the sum over all macrocell sides.

estimate_darcy_jumps(solution: DarcySolution, *, dirichlet: Any = 0.0, neumann: dict[int, Any] | None = None, ellipticity_lower_bound: Any = None, calibration: float = 1.0, order: int | None = None) -> DarcyJumpEstimator

Evaluate Araya et al. (2013), equations (5.1)--(5.3).

Use R_F=-jump(p)/2 on interior macrofaces, R_F=g_D-p on Dirichlet faces, and zero on Neumann faces. Then eta_F=c_lc_min||R_F||/sqrt(H_F), where H_F is the original macroface length and c_min is the square root of a certified lower eigenvalue bound for permeability. Integration splits at both fine-edge partitions and all skeletal breakpoints.

calibration is c_l, not a proven constant-one reliability factor. The paper's Section 6 uses 3, 7, 18 and 50 for uniform trace degrees 0--3. Its theorem assumes exact local lifts, regularity and represented boundary data. Finite local errors require an additional estimator, such as the weighted flux-reconstruction estimator; no bound is asserted here.

Nested local energy differences at a fixed physical MHM skeleton trace.

LocalDarcyRefinementEstimate dataclass

Energy differences between nested local solves with the same macro flux.

This is a computed refinement difference, not an error bound without a separately justified saturation estimate. Local means do not enter its norm.

total: float property

Return the square root of the sum of all local energy differences.

estimate_darcy_local_refinement(solution: DarcySolution, *, quadrature_order: int | None = None, local_solver: str = 'scipy', refinement_precision: Literal['double', 'extended'] = 'double', backend: Literal['serial', 'thread', 'process'] = 'serial', workers: int | None = None) -> LocalDarcyRefinementEstimate

Red-refine every local triangle and solve with unchanged physical face data.

The finite-element degree, source, permeability and skeletal coefficients remain fixed. The enriched local pressure has zero mean; the energy norm only compares gradients. Return the difference in the genuine broken K energy norm using material-cut quadrature and exact parent-cell ancestry. Point sources and mixed pressure spaces are outside this contract. No reconstructed-divergence term is used as a surrogate local-error norm. backend distributes independent local solves in cell order. Process workers use spawn and require picklable source/material callbacks; all native factors are constructed and released inside their owning worker.

Residual-indicator mesh sizes and optional FreeFEM/BAMG remeshing.

The nodal metric follows the scalar residual-indicator construction documented by FreeFEM. It is independent of the PDE solver and does not prescribe a bulk marking rule or preserve nested meshes.

ResidualMeshSize dataclass

Current and requested nodal sizes and the lumped indicator projection.

residual_mesh_size(mesh: TriangleMesh, local_squared: Any, *, coefficient: float = 1.0) -> ResidualMeshSize

Compute the documented isotropic residual metric from squared indicators.

The current size at a vertex is the average length of its incident edges, counted once per incident triangle. A continuous P1 indicator is obtained by area-weighted lumped projection of the cellwise square roots. Divide the current size by clip(projected / (coefficient * mean(cellwise)), 1/3, 3). Thus both refinement and coarsening are allowed. A zero indicator leaves sizes unchanged. Unused vertices and invalid data are rejected explicitly.

remesh_freefem(mesh: TriangleMesh, sizes: Any, *, executable: str = 'FreeFem++', max_vertices: int = 10000, timeout: float = 120.0) -> TriangleMesh

Remesh with BAMG using a supplied continuous P1 physical size field.

A separately installed FreeFEM executable is required. Only mesh generation is delegated: the call uses adaptmesh(IsMetric=1, splitpbedge=1) and does not solve a PDE. The requested size is interpolated in the input P1 space; the output need not be nested. Material interfaces and boundary data must be transferred by the caller's physical functions, not by cell-index ancestry. The unchanged domain area is checked after reading the returned mesh.

Face residuals of the one-level primal elasticity estimator.

The estimator follows Harder, Madureira and Valentin (2016).

The reliability theorem assumes exact local Neumann lifts. With computed local finite elements this is a face indicator, not a certified estimator of the full two-level error. In particular, an unresolved continuous local field can have a zero jump indicator and nonzero interior error.

PrimalElasticityIndicator dataclass

Per-segment displacement jumps scaled by the declared ellipticity constant.

face_squared contains c_min²/H_F times the squared displacement residual. Internal faces have the factor one half in the displacement jump and count twice in eta. The denominator is the original macroface length. The object records no assertion that finite local lifts satisfy the one-level reliability theorem; that additional hypothesis is not numerically inferred.

eta: float property

Return the face indicator with the macrocell-boundary multiplicity.

The indicator follows Harder, Madureira and Valentin (2016).

estimate_primal_elasticity_error(solution: PrimalElasticitySolution, *, dirichlet: Any = (0.0, 0.0), c_min: float, full_dirichlet: bool, quadrature_order: int = 8) -> PrimalElasticityIndicator

Evaluate the face estimator under an explicit full-Dirichlet contract.

Use equations (5.4)--(5.6) of Harder, Madureira and Valentin (2016).

Supply c_min such that A(x) epsilon:epsilon >= c_min² |epsilon|² everywhere. The supplied bound is checked at face integration points, but only the user can establish its global validity for a generic callable. The theorem's constants depend on the material; no incompressibility-uniform or unit reliability constant is claimed. For finite local approximations, refine the local space separately to assess the unmeasured local consistency error.

References

• Rodolfo Araya, Christopher Harder, Diego Paredes, and Frédéric Valentin (2013). Multiscale Hybrid-Mixed Method, SIAM Journal on Numerical Analysis 51(6), 3505–3531. DOI: 10.1137/120888223.

• 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.

• Christopher Harder, Alexandre L. Madureira, and Frédéric Valentin (2016). A hybrid-mixed method for elasticity, ESAIM: M2AN 50, 311–336. DOI: 10.1051/m2an/2015046.

Source modules

Conforming local closure for adaptively segmented flow skeletons.

refine_flow_local_meshes(skeleton: SkeletonSpace, local_meshes: tuple[TriangleMesh, ...], marked_local: np.ndarray, maximum_cells: int, *, fine_marked: tuple[np.ndarray, ...] | None = None) -> tuple[TriangleMesh, ...] | None

Resolve every new trace vertex while refining only the required local boundary cells.

Local-error dominance red-refines all triangles of that local mesh once, unless fine_marked supplies selected cells for longest-edge refinement. Subsequent longest-edge propagation contains every skeletal breakpoint, preserving conformity and avoiding uniform closure of unrelated interiors. Return None before a local cell cap is exceeded. The input meshes are never mutated, so the caller can retain its last solved configuration.