User-defined variational problems

The principal workflow binds macro/local meshes, declared spaces and user-written local/global forms with MeshHierarchy, bind_interface and bind_problem. The local provider receives LocalContext: it supplies the local mesh, native space binding, supported trace pairings and coordinate maps. The global callback receives GlobalContext for represented boundary data and additional forms. assemble and solve reuse the shared compilation, constrained elimination, ordered assembly and reconstruction owners.

The mathematical operator, physical trace meaning, boundary convention, retained modes and gauge remain explicit. Numbering, incidence signs and supported native coordinate conversions belong to their mesh/space adapters. Scalar, vector and mixed fields use the same contracts.

For the abstract MHM decomposition and local/global interpretation, see Harder and Valentin (2016). The four-block API below also accepts independently declared trial/test forms.

The current contract is a linear system with real coefficients. Nonlinear iterations and time integration can compose these systems; their update rules and scientific qualification belong to the application. Complex formulations use the general realify_operator, realify_vector and complexify_vector operations, with interleaved real and imaginary coordinates for each coefficient. This representation applies to every operator block and load.

The provider notebooks give executable scalar and mixed examples. The API gives the complete signatures. The case gallery also documents predefined physical formulations imported from their implementation owners. Their numerical qualification applies to the stated data and spaces, rather than arbitrary variational forms.

Bind spaces before writing the equations

hierarchy = MeshHierarchy(macro_mesh, local_meshes)
interface = bind_interface(skeleton, convention="normal")


def local_equations(local):
    space = local.native_space(local_element)
    u, v = ufl.TrialFunction(space.space), ufl.TestFunction(space.space)
    dx = ufl.Measure("dx", domain=space.mesh)
    a = (ufl.inner(K * ufl.grad(u), ufl.grad(v)) + reaction * u * v) * dx
    L = source * v * dx
    b = local.trace_pairings(lambda phi, ds: phi * v * ds)
    c = local.trace_pairings(lambda phi, ds: -phi * u * ds, axis="rows")
    local.field("potential", space)
    return local.equations(a=a, L=L, b=b, c=c)


problem = bind_problem(
    hierarchy, interface, local_equations,
    global_equation=global_equation, retained=retained_modes,
)
system = assemble(problem)
solution = solve(system)
potential_fields = solution.field("potential")

This snippet assumes the displayed coefficient symbols are defined on the native domain and the declared operator has no unhandled kernel. A complete runnable example, including a singular local operator and exterior pressure moments, appears in the overview.

retained is the declared mode count per macroelement, either one count for all selected cells or a tuple following the hierarchy order. The binding derives coordinate sizes and offsets without solving local problems on the coordinator. Native trace pairing capabilities have stated geometric/space limits; generic UFL does not automatically couple unrelated meshes.

For an additional supported global interface form, pass a builder to global_problem.interface_equation. Its arguments are the native interface trial/test functions and measure; the adapter transports the assembled form to the declared face coordinates and retains the complete global layout. Boundary loads and gauges are separate explicit data.

The built-in native trace adapters support polynomial edges and triangular faces in 2D and 3D. ComponentTraceSpace(scalar_space, components=3) declares Cartesian vector values without changing the scalar space's continuity. TangentialTraceSpace instead uses the declared canonical tangent frame. Facewise C0 Bernstein coefficients remain shared within their macroface; PressureTraceSpace and PressureTraceSpace3D also share macrovertices and macroedges. Native boundary facets must resolve the selected trace partition. This geometric condition does not replace a method's degree or stability conditions. On the triangular global integration mesh, subfaces are separate geometric cells; coefficient embeddings enforce the declared continuity. Cross-face facet terms are supplied explicitly as additional global blocks.

For already assembled formulations, the same basis is available directly:

from pymhm import project_trace, trace_bilinear_form, trace_linear_form

boundary_load = trace_linear_form(skeleton, prescribed_value, faces=boundary_faces)
boundary_operator = trace_bilinear_form(skeleton, coefficient, faces=boundary_faces)
fixed_coefficients = project_trace(skeleton, prescribed_trace, faces=essential_faces)

These operations integrate the physical scalar/vector value represented by the trace. They choose no pressure, normal-flux or traction meaning and insert no incidence sign. The coefficient in trace_bilinear_form can be a scalar or a matrix in physical value components. Shared nodes accumulate every selected face contribution; project_trace solves their joint Gram system. The same operations serve reaction, Robin, impedance and other interface terms.

The custom-space guide shows how to own all numbering, basis and orientation conventions. Fully manual coefficient records remain available below and use the same numerical owners.

Explicit coefficient equations

A provider returns the following four blocks for each local item:

\[ \begin{aligned} a_K(u_K,v_K)+b_K(\lambda_K,v_K)&=L_K(v_K),\\ c_K(u_K,\mu_K)+d_K(\lambda_K,\mu_K)&=g_K(\mu_K). \end{aligned} \]

Their coefficient representation is

\[ \begin{aligned} A_Ku_K+B_K\lambda_K&=f_K,\\ C_Ku_K+D_K\lambda_K&=g_K. \end{aligned} \]

Rows are test coordinates and columns are trial coordinates. b and c are independent: the package does not infer an adjoint, sign or physical boundary variable. dofs maps the columns of \(B_K\) to global trace coordinates; test_dofs maps the rows of \(C_K\). The maps can differ in order and length. d has shape (len(test_dofs), len(dofs)) and g has length len(test_dofs). Omitting test_dofs uses dofs.

from pymhm import Equation, LocalEquations, MultiscaleProblem, solve

def local(cell):
    return LocalEquations(
        a=[[2.0]], L=[1.0], b=[[1.0]], c=[[-1.0]],
        d=[[1.0]], g=[0.0], dofs=[0],
    )

problem = MultiscaleProblem(
    global_equation=Equation(0, 0),
    local_provider=local,
    items=[0],
    trace_size=1,
    coarse_sizes=(0,),
)
solution = solve(problem)
assert abs(solution.trace[0] - 1 / 3) < 1e-14
assert abs(solution.fields[0][0] - 1 / 3) < 1e-14

This is an algebraic illustration, not a discretization study. For invertible \(A_K\), local elimination contributes

\[ \begin{aligned} S_K&=D_K-C_KA_K^{-1}B_K,\\ r_K&=g_K-C_KA_K^{-1}f_K. \end{aligned} \]

The implementation factors the local operator once for the source and coupling columns. It uses the same checked condensation owner as the established formulations. A global Equation(a_G, L_G) adds an independent global operator and load after all local contributions. It can contain cross-cell couplings, boundary terms or an already assembled facet operator.

If \(P_K\) gathers a cell's trial coordinates and \(Q_K\) gathers its test coordinates, the reduced global equation is

\[ a_G(x,y)+\sum_K(Q_Ky)^TS_K(P_Kx) =L_G(y)+\sum_K(Q_Ky)^Tr_K. \]

The additional global matrix and vector use the complete reduced coordinate order: all trace coefficients first, then each item's retained coefficients in items order. They must have that full shape. Equation(0, 0) denotes zero additions with dimensions supplied by the problem.

The assembled system is a data object. solve_multiscale_system(system) solves its declared boundary and gauge problem; reconstruct_multiscale(system, coordinates) recovers fields from executed coefficients without another global solve. Its solve and reconstruct methods delegate to those free functions. with_global_load(system, rhs) reuses the local responses for a new global balance load while preserving retained source-compatibility rows. It does not change a local volume source.

For repeated volume sources, preserve the operator and its factors explicitly:

from pymhm import OfflineMultiscaleSystem

with OfflineMultiscaleSystem(system) as offline:
    updated = offline.with_loads(local_loads, global_load=additional_global_load)
    solution = solve(updated)

local_loads contains one new \(L_K\) in each assembled test-coordinate order. Optional balance_loads replaces each \(g_K\) in its compiled trace-coordinate union. global_load replaces the additional global functional, rather than the complete condensed right-hand side. Local source responses and retained compatibility rows are recomputed while harmonic lifts and basis matrices are reused literally. A source-dependent physical gauge is computed from updated before solving. This operation supplies no time integrator or source projection. Leaf systems use this cache directly; recursive sources are updated in their child hierarchy before the parent is assembled. The ordinary factorization solver argument also permits explicit global factor reuse.

Named fields can combine eliminated, trace and prescribed coordinates. A FieldDefinition stores the explicit reconstruction

\[ w_K=R_Ku_K+T_K\lambda[\mathrm{trace\_dofs}_K]+w_K^0. \]

Declare those maps with reconstruction, trace_reconstruction, trace_dofs and offset in local.field(...), or use the public pymhm.postprocessing.nodal.nodal_field helper for portable nodal fields. The trace map uses the stated global coordinates; include any required orientation in that matrix. This supports hybrid moment reconstructions and essential-boundary liftings without a method-specific solution object. solution.field(name) returns a tuple of recovered DiscreteField views, each using its archived basis. Evaluate a view at physical points with an explicit cell selection when one-sided interface values are needed:

views = solution.field("pressure")
values = views[0].evaluate(points, cells=local_cells)

Field names do not select an operator, gauge or reconstruction formula. For a hierarchy, solution.field(name, recursive=True) descends into cells without a definition and returns their recovered leaf fields in deterministic order. Each view uses the trace vector at its own level. A cell declaring the field takes precedence over its children, so the result does not duplicate it.

FieldDefinition.basis_matrix exposes an immutable copy when the declared nodal, modal or native basis has a single polynomial matrix. A composite Piola basis additionally needs moment maps and orientation matrices; archive the complete FieldDefinition and its basis_digest for replay. The public nodal_field, modal_field, piola_field and hdiv_field constructors retain their executed coordinate contracts and support one-sided evaluation. Explicit cell owners evaluate only the requested point/cell pairs, with storage proportional to the number of points. For repeated queries on a planar triangle mesh, TrianglePointLocator provides a reusable search and validates each chosen cell's barycentric coordinates.

Formulations composed from the same API

The importable providers under examples/formulations/ state the local spaces, operators, boundary data and recovery maps. Their main paths compose public numerical operations with LocalEquations, Equation, assemble and solve. They require no predefined physical solver or private numerical helper.

Formulation family	Provider and explicit mathematical ingredients
Primal Darcy, tensor diffusion and Cartesian/tetrahedral variants	`darcy`, `tensor_darcy`, `cartesian_darcy`, `tetrahedral_darcy`: stiffness, signed normal pairings, constant kernel and pressure moment.
Mixed Darcy on triangles, tetrahedra, prisms and hexahedra	`mixed_darcy`, `mixed_darcy_3d`, `mapped_darcy`: inverse-material mass, divergence, normal-flux constraints and pressure gauge; trilinear hexahedra use the mapped provider.
Analytical metric-quadratic Darcy	`analytic_darcy`: exact scalar energy and signed face moments, polynomial source response and separate mean-zero nodal source potential.
Robin MH	`robin`: volume and weighted boundary bilinear forms, independently signed trace equations and physical pressure moment.
Three-field MH²M	`three_field`: independent pressure/conormal spaces, moment complement and their L2 pairing.
MsHHO	`moments`, `moments_3d`: cell/face moments and energy reconstruction, including its trace-dependent field map.
Residual Petrov–Galerkin MHM and stationary transport	`penalty`, `scalar`, `transport`, `residual_transport`, `scalar_3d`, `transport_3d`: residual operators, test/trial spaces, stabilization and explicit additional face terms, with declared strong or weak boundary data.
Displacement, pressure and weak-stress elasticity	`elasticity`, `tetrahedral_vector`, `primal_elasticity`, `weak_stress`, `vector`: strain/compliance operators, rotation or pressure constraints, rigid modes and their physical moments.
Stokes, Brinkman and Oseen	`flow`, `vector`: velocity/pressure pairings, reaction/advection, stated stabilization and pressure gauge.
Helmholtz and PML	`examples/tutorial_helmholtz_equations.py`: complex physical blocks represented in general real coordinates, boundary impedance and point loads.
Maxwell and elastodynamics	Wave tutorial providers: mass/curl or mass/strain operators and explicit time-step equations with reusable factors.
Hierarchical and adaptive variants	`NestedEquations`, `adaptive`: declared child restrictions, recovered physical fields, estimators and marking.
Classical conforming references	`conforming`: element or exact separated-factor energy, strong nodal data, outward-flux load and physical integral gauge.

A fully global formulation declares MultiscaleProblem.from_global(Equation(A, L), size, fixed=..., constraints=...). It uses the same assembly and solver contracts without a dummy local provider. The numerical owners also permit direct composition with solve_linear when a user wants to manage the coefficient elimination explicitly.

Adaptive controllers accept an optional solve_step callable. It receives the current mesh, transferred skeleton and declared problem options, and returns the physical solution used by the estimator. The callable can assemble and solve user-written LocalEquations through the same API; the shared controller owns estimation, marking and mesh transfer. Omitting it retains the existing predefined formulation. examples/formulations/adaptive.py shows this composition with an editable mathematical definition factory.

The notebooks under the corresponding physical-problem directories are the executable entry points. The notebook catalogue gives locked environments and reproduction commands. An analytical control qualifies its stated data and spaces; literature reproductions and historical external comparisons retain their own data and provenance requirements.

Use UFL or assembled blocks

a, L, b, c, d, g and the global Equation accept assembled arrays or sparse operators. compile_form also accepts real linear and bilinear UFL forms through the optional DOLFINx adapter. A scalar zero obtains its dimensions from the surrounding block; nonzero scalars are not broadcast.

For example, on a declared scalar space V, a reaction–diffusion volume form can be written directly:

import ufl
from pymhm import LocalEquations, columns, rows

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
dx = ufl.Measure("dx", domain=local_mesh, metadata={"quadrature_degree": 6})
ds = ufl.Measure("ds", domain=local_mesh)
a = (ufl.inner(K * ufl.grad(u), ufl.grad(v)) + reaction * u * v) * dx
L = source * v * dx

local = LocalEquations(
    a=a,
    L=L,
    b=columns(*(phi * v * ds for phi in signed_trace_functions)),
    c=rows(*(-phi * u * ds for phi in signed_trace_functions)),
    dofs=trace_dofs,
)

V, coefficients, measures, integration regions and the signed trace functions are supplied by the provider. Each phi represents one explicitly chosen trace basis function on the local integration mesh, including its support and orientation. columns assembles test-side linear forms as \(B_K\); rows assembles trial-side linear forms as \(C_K\). This direct pairing path supplies its own support and coordinate convention. The bound LocalContext.trace_pairings path supplies those data for supported mesh-associated spaces; arbitrary cross-mesh identification remains explicit.

The scalar example has no declared kernel. If reaction=0 and the local volume operator is a Neumann diffusion operator, provide its constant coefficients and physical mean moments as described below.

Vector and mixed UFL forms follow the same contract. For a mixed velocity and pressure space W, one possible symmetric Brinkman volume form is

u, p = ufl.TrialFunctions(W)
v, q = ufl.TestFunctions(W)
a = (
    2 * viscosity * ufl.inner(ufl.sym(ufl.grad(u)), ufl.sym(ufl.grad(v)))
    + ufl.inner(resistance * u, v)
    - p * ufl.div(v) - q * ufl.div(u)
) * dx
L = ufl.inner(body_force, v) * dx

This declares the volume operator. The trace may represent physical traction or a method-specific multiplier according to the separately declared boundary forms. Stable velocity/pressure spaces, local kernels, compatible data and the physical pressure gauge still have to be supplied. The API does not choose Taylor–Hood, stabilization or a boundary prescription from the field names.

Mixed \(H(\mathrm{div})\times L^2\) volume forms can similarly use inner(K_inverse * flux, test_flux), divergence pairings and a pressure source. Flux-prescribing MHM also needs its boundary variables and normal-trace constraints. A mixed volume form alone does not define that construction.

Native compilation limits

pymhm.backends.forms.assemble_form accepts real rank-one and rank-two UFL forms on single-rank DOLFINx domains, normally MPI.COMM_SELF. Bilinear forms may have distinct test and trial spaces and rectangular shapes; the chosen local pivot in LocalEquations must nevertheless be square. A declared 0-by-0 pivot has no eliminated field coordinates: its d and g contribute directly without a local factorization. Linear forms with UFL argument number zero or one are supported. The complete reduced global operator is also square.

The direct adapter assembles binary64 coefficients. Essential boundary elimination, gauge, quadrature and transpose relations remain explicit. The higher-level context applies its declared interface orientation maps once. Cross-mesh integration requires explicit DOLFINx entity_maps; use assemble_form or a custom compiler to pass these maps and compilation options. DOLFINx 0.9 uses a mapping from native meshes to entity-index arrays; DOLFINx 0.10 uses a sequence of native EntityMap objects. The adapter copies the supplied container to a dictionary or list respectively. Mesh correspondence, trace orientation and the underlying map data remain explicitly supplied inputs. Zero forms whose arguments were simplified away require an explicit shape. Complex-valued blocks, nonzero rank-zero functionals and higher-rank forms are outside this compiler contract. Declare complex formulations in explicit real and imaginary coordinates, as in the Helmholtz notebook; the compiler assembles the resulting real blocks.

Retain local modes and physical constraints

For a singular local operator, declare kernel=Z and its physical column pairings moments=C. For distinct left modes, supply left_kernel, test_moments and the explicit test-side pairing. These arrays use the actual local coefficient basis. They are checked without rediscovering or rotating the basis. A nonsingular moment pairing is required for the chosen modes.

moment_complement(C, pivots=...) constructs coordinates satisfying C.T @ coefficients == 0. Columns of C are the declared linear functionals; free coordinates use identity rows in ascending order, and the pivot rows enforce the moments. For example, quadrature integrals of basis functions define a zero-integral complement. The operation is equally applicable to several conservation or orthogonality constraints. Specify admissible pivots for a fixed coordinate convention and archive the executed basis when saving its coefficients.

coarse_basis can retain modes that need not be exact null vectors, with an independent test_basis when needed. This uses the same complementary elimination as the shared hybrid contract. kernel and coarse_basis are mutually exclusive. coarse_sizes declares the retained width of every item, including zero.

fixed prescribes explicit trace coefficients. constraints contains global physical rows and values. With assemble(problem), the returned system also provides mean_constraint for physical integral weights in reconstructed fields; then pass that row to system.solve(constraints=[...]). Integrals must include the correct physical fields and exclude auxiliary multipliers. No mean pressure, rigid-motion constraint or boundary datum is inferred.

These checks establish consistency with the represented algebraic contract. They do not establish mesh-independent inf-sup stability, uniqueness of a continuous PDE or a convergence theorem. Keep physical field errors and original equation residuals separated by their field blocks in scientific verification.

Assemble serially or in bounded batches

from pymhm import ExecutionConfig, assemble

system = assemble(
    problem,
    execution=ExecutionConfig(
        backend="process", workers=10, batch_size=10, native_threads=1,
    ),
)
solution = system.solve()

Serial execution builds, eliminates and contributes one item before taking the next. Thread or spawn-process execution uses bounded batches. Only the coordinator reduces shared global coordinates, in item order, so two cells contributing to the same face do not concurrently mutate a global matrix. Items, providers, compilers and callbacks must be picklable for spawn. Executable process entry points use the usual __main__ guard.

Create native meshes and forms inside the provider invocation. Compile and release them there before returning CompiledLocalEquations with owned arrays and portable metadata. Live native meshes, matrices and communicators must not cross a worker boundary. SolverConfig.local_solver can use an external numerical solver or response model, subject to the unchanged original-equation checks. A learned field without the required operator, coordinates and response columns does not satisfy that solver contract.

Reuse the interface at finer scales

Set LocalEquations.a to another MultiscaleProblem when its condensed global operator is the local operator at the parent scale. Its reduced coordinates become the parent's local coefficient basis; L adds a parent-scale source to the child's assembled load. Parent trace pairings and moment arrays use those same coordinates.

Child assembly runs serially inside its owning worker, avoiding nested process pools. The child leaves prescribed trace coefficients and physical gauges to its parent: child fixed and constraints must be empty. Express required local essential constraints as explicit local blocks. Parent coupling and additional loads must leave the child's retained compatibility rows unchanged.

NestedEquations provides the equality-constrained boundary construction for a child hybrid problem. With child equation \(A x=f\), boundary selector \(E\) and an explicit oriented restriction \(M\), its local system is

\[ \begin{bmatrix}-A&E\\E^T&0\end{bmatrix} \begin{bmatrix}x\\\mu\end{bmatrix} +\begin{bmatrix}0\\-M\end{bmatrix}\lambda =\begin{bmatrix}-f\\0\end{bmatrix}. \]

boundary_dofs selects distinct child trace coordinates, excluding retained cell coordinates. trace_map must be injective on the parent trace space and contain every orientation and basis restriction. Declare kernels and moments in the child's reduced coefficient basis; the shared nesting operation lifts their reaction coordinates and extends moments by zero reaction weights. The physical nesting moment is \(-M^T\mu\). The shared hybrid saddle negates that pairing in its reduced trace equations, so the four-block convention on this page has \(C[x;\mu]=+M^T\mu\). Global loads use that same row convention. Reconstruction checks the child's actual load \(f+E\mu\). This construction also requires empty child fixed and constraints.

MultiscaleSolution.children[cell] contains a recovered child solution, or None at a leaf. The hierarchy reconstructs against the actual parent-induced load without solving the child's global system again. This is recursive algebraic elimination; an error estimate for a particular multilevel method requires its own space, regularity and compatibility hypotheses.

leaf_moment(system, weights) represents a physical integral in finest-scale fields as row @ coordinates + offset. The callback supplies integral weights in each leaf's executed coefficient basis. The operation includes source offsets and pads boundary-reaction coordinates by zero. A prescribed physical moment value therefore uses the global constraint (row, value - offset). It does not replace that integral by a mean of coarse coefficients or include auxiliary reactions in a physical pressure gauge.

The three-level coefficient notebook checks operator recursion against an independently assembled full matrix. The spatial recursive MHM notebook declares Q2 leaf forms, P1 trace restrictions and NestedEquations, with a physical leaf-integral gauge for Neumann data. Its live analytical checks are separate from the archived convergence acquisitions and their provenance.

Represent distinct method constructions

The interface does not select a method by name. These recipes identify the blocks and data needed to express the published constructions; the existing method implementations provide their separately verified cases.

Construction	Required form and space choices
Neumann-local MHM — Harder, Paredes and Valentin (2013)	Local operator on the complement of its declared kernel; normal-flux trace, source lift, physical moments and retained compatibility equations.
Robin-local MH — Barrenechea, Gomes and Paredes (2024)	Robin terms inside the local operator and their matching global balance; the multiplier must be distinguished from physical flux.
Three-field MH²M — de Barros, Madureira and Valentin (2026, v3)	Local solution and duplicated conormal coefficients as mixed local fields; continuous pressure trace globally; complete source lifting and a boundary-mean local normalization.
MsHHO — Chaumont-Frelet, Ern, Lemaire and Valentin (2022)	Constrained energy reconstruction with declared cell and face moments; a local moment saddle and subsequent cell/face condensation, with the selected source projection.
Residual Petrov–Galerkin MHM — Fernando, Martins, Pereira and Valentin (2023)	Independently declared trial/test pairings, enriched local responses and cross-cell macroface terms, with the published stabilization and reconstruction.

For MsHHO, the reconstruction uses the mixed block

\[ \begin{bmatrix}A_K&C_K\\C_K^T&0\end{bmatrix} \begin{bmatrix}u_K\\\eta_K\end{bmatrix} =\begin{bmatrix}0\\m_K\end{bmatrix}. \]

Here \(C_K\) denotes the declared moment columns and \(m_K\) the prescribed cell and face moments, independently of the balance block notation above. The energy matrix and subsequent elimination are part of this method's construction. They are not obtained merely by changing the name of a flux trace. The free function energy_reconstruction(A, C) supplies the moment reconstruction \(R\) and its reduction \(R^T A R\) through the shared original moment saddle. Its energy-minimization interpretation requires symmetry and coercivity on the moment-zero space; nonsymmetric input preserves the full bilinear reduction. with_global_equation(system, equation) adds forms after the executed local response maps are available, so an explicitly assembled jump of reconstructed fields can supply a Petrov facet term. It preserves those response maps, bases and children in a new system. Automatic coupling of reconstructed fields on different meshes, jump signs and penalties are not inferred by Equation.

The literature imposes additional conditions. The MHM–MsHHO equivalence in Chaumont-Frelet, Ern, Lemaire and Valentin (2022) needs exact local solves and the prescribed polynomial source space; the source-projected variant has a different equivalence statement. The original equivalence also excludes the face-only \(m=-1\) variant.

The three-field analysis of de Barros, Madureira and Valentin (2026, v3) requires its local injectivity and uniform norm bounds, not just an observed full-rank matrix.

The residual Petrov–Galerkin error estimates of Fernando, Martins, Pereira and Valentin (2023) include dimension-dependent enrichment and stabilization bounds: the stated theorem uses trace degree \(\ell\geq1\), local degree \(k\geq\ell+d\) in dimension \(d\), and \(0<\alpha\leq\alpha_0\). The admissible upper constant is part of the analysis, not a numerical value inferred from a successful solve.

Mixed \(H(\mathrm{div})\) methods require compatible divergence, pressure and normal-trace spaces; see Durán, Devloo, Gomes and Valentin (2019) for the mixed-local multiscale construction. The primary literature and scope give the specific hypotheses and the case gallery gives executed evidence. None of these theorems follows from the ability to assemble four blocks.

References

Christopher Harder and Frédéric Valentin (2016). Foundations of the MHM Method, in Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations, Lecture Notes in Computational Science and Engineering 114, Springer. DOI: 10.1007/978-3-319-41640-3_13.
Christopher Harder, Diego Paredes, and Frédéric Valentin (2013). A family of Multiscale Hybrid-Mixed finite element methods for the Darcy equation with rough coefficients, Journal of Computational Physics 245, 107–130. DOI: 10.1016/j.jcp.2013.03.019.
Gabriel R. Barrenechea, Antonio Tadeu A. Gomes, and Diego Paredes (2024). A Multiscale Hybrid Method. SIAM Journal on Scientific Computing 46(3), A1628–A1657. DOI: 10.1137/22M1542556.
Franklin de Barros, Alexandre L. Madureira, and Frédéric Valentin (2026). A three-field Multiscale Method. arXiv preprint, version 3, 5 August 2026; first submitted 25 April 2024. arXiv: 2404.16978v3.
Théophile Chaumont-Frelet, Alexandre Ern, Simon Lemaire, and Frédéric Valentin (2022). Bridging the Multiscale Hybrid-Mixed and Multiscale Hybrid High-Order Methods, ESAIM: M2AN 56, 261–285. DOI: 10.1051/m2an/2021082.
Honório Fernando, Larissa Martins, Weslley Pereira, and Frédéric Valentin (2023). A Petrov–Galerkin multiscale hybrid-mixed method for the Darcy equation on polytopes. Computational and Applied Mathematics 42, article 173. DOI: 10.1007/s40314-023-02304-y.
Omar Durán, Philippe R. B. Devloo, Sônia M. Gomes, and Frédéric Valentin (2019). A multiscale hybrid method for Darcy’s problems using mixed finite element local solvers, Computer Methods in Applied Mechanics and Engineering 354, 213–244. DOI: 10.1016/j.cma.2019.05.013.