Hybrid operators and multiscale constructions

Field queries can supply explicit cell owners for independent one-sided values. This evaluates only the requested point/cell pairs. TrianglePointLocator provides reusable geometric searches for planar triangular fields.

User-defined variational equations, local condensation, recursive reconstruction and distinct multiscale formulations. The variational guide describes the four blocks, global additions and current compilation limits.

All API families

Mesh, interface and variational contexts

Portable interface layouts and transport of user-defined variational blocks.

An interface binding maps selected global coefficients into the declared local trial/test trace bases. Test functionals travel through the transpose map. The maps may encode signs, permutations or dense basis changes; they never select a physical coupling sign. Numerical condensation remains owned by equations.

TraceBinding dataclass

Independent trial/test embeddings into a square global interface layout.

trial_map @ global[dofs] gives local trial trace coefficients, and test_map @ global[test_dofs] gives local test coefficients. Each map has one column per selected global coordinate. Rectangular and rank-deficient maps permit projections/restrictions; require_injective=True requests an injectivity check when the formulation needs it. Omitted maps are identities; omitted test coordinates/maps reuse the trial binding. basis_id names the executed local basis supplied by an adapter. Its identity, maps and selected coordinates enter basis_digest together. The digest identifies declared representation data, not PDE consistency.

trial_size: int property

Return the local trial basis width before transport.

test_size: int property

Return the local test basis width before transport.

__post_init__() -> None

Freeze owned coordinate maps and validate independent basis transports.

__reduce__() -> tuple[Any, tuple[Any, ...]]

Preserve immutable coefficient maps and basis identity across spawn/replay.

InterfaceSpace

Bases: Protocol

Describe global interface coordinates without requiring inheritance.

A custom implementation owns its basis convention and local support. It can use arbitrary numbering and dense changes of trace basis. Geometry, basis evaluation and native form binding are separate optional capabilities.

size: int property

Return the complete global interface coefficient dimension.

binding(cell: int) -> TraceBinding

Return independent trace maps for the selected macrocell.

ComponentTraceSpace dataclass

Interleave scalar trace coordinates for a declared vector-valued interface.

base owns the scalar basis, geometry and any continuity constraints. Coordinate components*i+j is component j of scalar coordinate i. This wrapper shares exactly the same vertices/edges as the base and introduces no normal or tangential rotation. Use bind_interface to declare value or canonical-normal transport, or TraceBinding for a different frame.

mesh: Any property

Return the unchanged scalar interface's macro mesh.

size: int property

Return the component-interleaved global coefficient dimension.

__post_init__() -> None

Validate a scalar mesh-associated layout without loading a FEM backend.

dofs(face: int) -> np.ndarray

Return one face's vector coefficients, including any shared scalar nodes.

cell_dofs(cell: int) -> np.ndarray

Preserve the scalar cell-coordinate order with interleaved components.

BoundInterface dataclass

Adapt an existing face space to the structural interface-space contract.

convention='value' preserves scalar/vector face values. 'normal' transports coefficients expressed in canonical face normals to outward macrocell normals, reusing the mesh's existing incidence signs. The adapter never infers a physical sign in B or C. Existing operators that already include orientation must use coordinates='global' when bound.

mesh: Any property

Return the macro mesh associated with the underlying trace basis.

size: int property

Return the complete scalar/component-interleaved interface dimension.

__post_init__() -> None

Require a mesh-associated trace layout and a supported convention.

binding(cell: int) -> TraceBinding

Derive local coefficient transport from the existing topology owner.

MeshHierarchy dataclass

Associate macro cells with worker-created local meshes.

local_meshes is a callable cell -> mesh or one portable mesh per macrocell. The callable runs only when a local context requests its mesh, so native objects can remain worker-owned. items selects macrocell indices and their deterministic assembly order; its default selects all. Native objects, factories and meshes retain their documented execution requirements: process execution needs picklable portable descriptions.

__init__(macro: Any, local_meshes: Callable[[int], Any] | Sequence[Any], items: Sequence[int] | None = None) -> None

Check topology association without eagerly constructing local meshes.

local_mesh(cell: int) -> Any

Construct or retrieve a local mesh belonging to this hierarchy.

validate_trace_binding(binding: TraceBinding, size: int) -> TraceBinding

Reject stale or out-of-range layouts before assembling local equations.

bind_interface(space: Any, *, convention: Literal['value', 'normal'] = 'value') -> BoundInterface

Bind a built-in face space using its declared value or normal convention.

Custom spaces implement InterfaceSpace directly and return their own TraceBinding; no inheritance, registry or native FEM dependency is needed. Tangential/component rotations require an explicitly declared custom map.

bind_local_equations(binding: TraceBinding, *, a: Any, L: Any, b: Any, c: Any, d: Any = 0, g: Any = 0, kernel: Any = None, moments: Any = None, coarse_basis: Any = None, left_kernel: Any = None, test_basis: Any = None, test_moments: Any = None, metadata: Any = None, coordinates: Literal['local', 'global'] = 'local', compiler: FormCompiler = compile_form) -> LocalEquations

Transport trace pairings independently before the existing condensation.

For maps R and S, local blocks become B R, S.T C, S.T D R and S.T g. Local volume coefficients, retained modes and their moment functionals are unchanged. coordinates='global' declares that supplied pairings already use selected global coefficients and skips this transport, preventing double orientation for existing signed coupling owners. Native forms compile in their own local integration spaces through compiler.

Mesh-associated variational contexts lowering to the existing hybrid algebra.

Mathematical signs, spaces, modes and gauges are declarations. Contexts own their coordinate binding; condensation, reduction and solving keep their current owners.

GlobalContext dataclass

Bound interface and declared retained coordinates for global forms.

Coordinates are interface first, then retained modes in hierarchy item order. The counts describe user-declared physical modes; no local PDE is solved or nullspace guessed to determine this layout. Physical load signs remain explicit.

trace_size: int property

Return the number of declared global interface coordinates.

size: int property

Return the square global dimension, including declared retained modes.

__post_init__() -> None

Freeze declared dimensions independently of a mutable custom provider.

trace_load(value: Any) -> FloatArray

Embed an interface functional, leaving retained compatibility loads zero.

interface_equation(forms: Callable[[Any, Any, Any], Any]) -> Equation

Assemble user-written global interface UFL forms in the bound coordinates.

forms(trial, test, measure) returns Equation or its (a, L) pair. The built-in adapter integrates supported edge/triangle polynomial spaces; custom spaces may declare the same capability. Retained-mode rows and columns are padded by zero, without inventing a coupling or boundary sign.

retained_dofs(cell: int) -> np.ndarray

Expose a cell's retained coordinates for explicitly assembled global forms.

boundary_data(value: Any, neumann: Mapping[int, Any] | None = None, *, order: int = 6) -> tuple[FloatArray, dict[int, float]]

Integrate weak value data and project fixed normal-density coefficients.

This returns an unsigned interface functional and fixed coefficients; callers choose its sign in their global equation. Custom spaces may implement their own boundary_data(value, neumann, order=...) capability. No Dirichlet/Neumann role is inferred from a local differential operator.

fix_faces(faces: Any, value: float = 0.0) -> dict[int, float]

Prescribe a constant face value in a declared built-in scalar/vector basis.

This prescribes the interface variable itself, without interpreting a PDE boundary condition. Nonconstant or custom data use a declared projection.

LocalContext dataclass

One worker-owned macrocell's mesh, spaces and trace coordinate binding.

Local native spaces are created lazily. The context itself never travels back to the coordinator; equation metadata must contain portable evaluation data. Custom local solvers may use only mesh/binding and assembled numerical blocks.

mesh: Any property

Return the lazily constructed portable local mesh.

macro: Any property

Return the hierarchy's macro mesh for declared material/geometry data.

binding: TraceBinding property

Return this context's single executed coefficient-map snapshot.

__post_init__() -> None

Check selected entity and trace association without creating native objects.

native_space(element: Any = None, *, degree: int = 1, shape: tuple[int, ...] = ()) -> Any

Bind a user element, or equispaced Lagrange default, to this local mesh.

Native optional libraries load only on this call. Returned NativeSpace exposes its FunctionSpace and coordinate conversions for advanced use. Its resources belong to this context and are released by close().

trace_pairings(expression: Callable[[Any, Any], Any], *, axis: Literal['columns', 'rows'] = 'columns', interface: InterfaceSpace | None = None) -> Any

Bind a user-written boundary form to each supported local trace basis.

expression(trace_value, face_measure) returns its UFL linear form. The trace value has the declared scalar/vector shape; geometry orientation is applied once by equations(), not hidden in the expression. Built-in edge and triangular polynomial traces have a native adapter; custom spaces may supply a trace_pairings(context, expression, axis=...) capability explicitly. interface selects an independent local boundary basis when a formulation has local trace unknowns as well as its global interface. Returned forms use unsigned local basis coordinates; equations() transports only the primary global pairings, while local operator blocks retain their basis.

interface_pairing(test_interface: InterfaceSpace, *, order: int = 6) -> FloatArray

Pair an independent local boundary basis with the primary interface basis.

The built-in operation is the unsigned L2 boundary mass matrix, with local test coordinates in rows and local primary trial coordinates in columns. Face numbering, union partitions and shared-vertex accumulation belong to the trace owner. Users choose its signs and block placement. A custom test interface may define interface_pairing(context, order=...). Primary global maps are applied later by equations(), exactly once.

equations(*, coordinates: Literal['local', 'global'] = 'local', **forms: Any) -> LocalEquations

Bind declared forms without exposing global DOF numbering to the provider.

Local trace coordinates are transported through the interface maps. coordinates='global' explicitly retains already oriented existing blocks. No coupling transpose, retained mode, moment or physical sign is inferred.

field(name: str, space: Any = None, *, component: int | None = None, evaluator: Callable[..., Any] | None = None, gradient_evaluator: Callable[..., Any] | None = None, reconstruction: Any = None, trace_reconstruction: Any = None, trace_dofs: Any = None, offset: Any = None, mesh: Any = None, basis_id: str = '') -> None

Declare a named native or custom field without transferring live resources.

space is an optional NativeSpace; component selects a mixed subfield. Custom fields provide an evaluator and basis_id, or expose coefficients only. reconstruction maps the local solution into field data. The optional trace_reconstruction maps global trace[trace_dofs]; offset supplies an explicit affine lifting, for example essential boundary data.

nested(problem: MultiscaleProblem[Any], *, boundary_dofs: Any = None, trace_map: Any = None, **options: Any) -> NestedEquations

Bind a child problem to the parent interface through an exact restriction.

Automatic restriction reuses the existing planar normal-trace owner for bound SkeletonSpaces. Supplied boundary_dofs/trace_map declare a complete map from selected parent GLOBAL coefficients into child boundary trace coefficients; it is not oriented again. Child boundary/gauge, injectivity and parallel restrictions remain those of NestedEquations. No projection of an unrepresentable parent trace is silently introduced.

close() -> None

Release native spaces while retaining independently owned numerical data.

__enter__() -> LocalContext

Own native resources during a manually requested local context.

__exit__(exc_type: type[BaseException] | None, exc_value: BaseException | None, traceback: TracebackType | None) -> None

Release spaces on every exit; bound providers already manage this lifetime.

BoundProblem dataclass

Bases: MultiscaleProblem[int]

An existing multiscale problem with mesh-associated variational contexts.

local_context(cell: int) -> LocalContext

Create an explicit context for diagnostics or a manual provider invocation.

The caller owns this context and must close any native resources it creates. Normal assemble/solve calls manage that lifetime automatically in workers.

bind_problem(hierarchy: MeshHierarchy, interface: InterfaceSpace, local_provider: Callable[[LocalContext], Any], *, global_equation: Equation | Callable[[GlobalContext], Equation] = _EMPTY_EQUATION, retained: int | tuple[int, ...] = 0, fixed: Mapping[int, float] | Callable[[GlobalContext], Mapping[int, float]] | None = None, constraints: tuple[tuple[FloatArray, float], ...] = (), compiler: FormCompiler = compile_form) -> BoundProblem

Bind meshes/spaces and user forms to the existing assemble/solve pipeline.

retained declares the number of retained modes per selected cell; a scalar count repeats it, while a tuple follows hierarchy item order. The provider receives LocalContext rather than an integer. A global callback can use the declared layout to express additional interface/retained terms. Fully manual MultiscaleProblem and custom InterfaceSpace implementations remain available.

Named field views carrying the executed coefficient and evaluation contract.

FieldDefinition dataclass

Portable description of one physical field reconstructed from a local vector.

A native nodal descriptor stores its executed basis/map. A custom evaluator has signature evaluator(mesh, coefficients, points, , cells=None); its declared basis_id identifies the coordinate convention. reconstruction maps local coefficients to the evaluator's coefficients. Without an evaluator or descriptor, the field supports coefficient access but no spatial evaluation. The reconstructed coefficients are Rlocal + T*global_trace[trace_dofs] + offset. Every trace orientation and basis change is explicit in T; no incident normal sign or local binding is inferred. All maps and the affine offset form part of the archived coefficient and basis-digest contract.

basis_matrix: FloatArray property

Return an owned read-only copy of a declared single polynomial matrix.

Nodal/native and modal evaluators declare this coordinate matrix. A custom or composite Piola basis may require several matrices and then raises TypeError. Archive this entire definition with its digest to retain cell maps, orientation and reconstruction as well as the matrix.

basis_digest: str property

Identify the executed descriptor and any additional reconstruction map.

__post_init__() -> None

Own a reconstruction map and reject ambiguous evaluation conventions.

__reduce__() -> tuple[Any, tuple[Any, ...]]

Restore an owned immutable reconstruction through the validating constructor.

DiscreteField dataclass

One reconstructed field associated with its mesh and executed basis.

Coefficients are independently owned and read-only. Separate macrocell fields preserve one-sided interface values; evaluation does not average or smooth.

mesh: Any property

Return the declared physical mesh, without constructing native resources.

basis_digest: str property

Return the field's executed basis/reconstruction identity.

portable_coefficients: FloatArray property

Delegate coefficient access in the descriptor's recorded portable order.

__post_init__() -> None

Freeze the reconstructed field's finite real coefficients.

__reduce__() -> tuple[Any, tuple[Any, ...]]

Restore immutable coefficients through the same checked field constructor.

evaluate(points: Any, *, cells: Any = None) -> FloatArray

Delegate spatial evaluation with optional explicit one-sided cell owners.

gradient(points: Any, *, cells: Any = None) -> FloatArray

Delegate physical gradients, retaining explicit one-sided cell owners.

values_and_gradient(points: Any, *, cells: Any = None) -> tuple[FloatArray, FloatArray]

Delegate values and physical gradients with a shared native sampling pass.

evaluate_field(field: DiscreteField, points: Any, *, cells: Any = None) -> FloatArray

Evaluate a declared native or custom field; reject absent capabilities.

evaluate_field_gradient(field: DiscreteField, points: Any, *, cells: Any = None) -> FloatArray

Evaluate physical derivatives in the executed native or custom field basis.

Scalar gradients have axes (point, coordinate); vector gradients use (point, component, coordinate). These are raw derivatives, not a conservative H(div) flux reconstruction. Custom fields declare gradient_evaluator with the same arguments as their value evaluator; absent capabilities fail.

portable_field_coefficients(field: DiscreteField) -> FloatArray

Return coefficients through the executed descriptor's portable coordinate map.

The result is independently owned and read-only, with mixed components selected by their stored map. This only changes ordering; it does not regenerate the basis or discard its digest. These are nodal values when the executed basis is the corresponding nodal basis. A custom field without this declared map has no automatic portable coefficient conversion.

evaluate_field_and_gradient(field: DiscreteField, points: Any, *, cells: Any = None) -> tuple[FloatArray, FloatArray]

Sample values and raw physical derivatives without duplicate native geometry.

Native descriptors share cell location, geometry pullback and executed basis tabulation. Custom fields delegate to their independent value and gradient callbacks. Output axes and one-sided conventions match evaluate_field and evaluate_field_gradient; this does not construct a conservative flux.

solution_field(solution: Any, name: str, *, recursive: bool = False) -> tuple[DiscreteField, ...]

Recover named fields in deterministic macrocell order.

With recursive=True, a cell without its own definition delegates to its recovered child solution. Cells declaring the field are visited once; their children are not also returned. Each reconstruction uses its own level's trace vector, retaining one-sided values and executed orientations.

local_trace(solution: Any, index: int, *, test: bool = False) -> FloatArray

Recover selected global coefficients in their executed local trace basis.

index is a position in deterministic assembly item order, like solution.fields. The trial map is used by default; test=True uses the independent test map. Both maps already declare geometric orientation and basis changes. This operation does not infer a physical boundary sign or average adjacent fields. Explicit equations without a stored binding have no local-trace capability.

Validated one-sided ownership for large physical point samples on triangular meshes.

TrianglePointLocator

Locate points through centroid candidates and actual barycentric containment.

Candidates are a search acceleration, never a nearest-cell interpolation. A vectorized search tests candidates nearest centroids; unresolved points use all cells, including skew or strongly graded meshes. Among valid candidates the first cKDTree result is chosen. Shared interfaces therefore select one incident side deterministically for the same mesh and query; coordinates(..., cells=...) selects an explicit independent side. Containment uses physical orientation predicates with the declared nonnegative, dimensionless barycentric tolerance. Roundoff-ambiguous signs are resolved exactly for the represented binary coordinates; zero tolerance therefore retains vertices and edges while excluding points strictly outside. Coordinates are neither moved nor averaged. Returned barycentric coordinates retain their original binary64 affine tabulation and can have roundoff-size negative entries at a geometrically valid edge. A conservative mesh bounding box excludes strictly exterior points, and bounded batches keep exhaustive fallback storage independent of query size.

__init__(mesh: TriangleMesh, *, candidates: int = 4, tolerance: float = 1e-10) -> None

Cache the actual affine maps; validate the geometric query convention.

locate(points: Any, *, allow_outside: bool = False) -> IntArray

Return one cell per point, optionally marking exterior points with -1.

By default every point must lie in the mesh within the declared barycentric tolerance. allow_outside=True preserves the query order and returns -1 only for points outside every triangle; this supports masks when sampling independent subdomains. Exterior indices must be filtered before evaluating a cell field.

coordinates(points: Any, *, cells: Any = None) -> tuple[IntArray, FloatArray]

Return owners and barycentric triples, optionally preserving declared incident sides.

Portable named nodal fields using the recorded polynomial basis and topology.

Evaluation uses archived Basix coefficient matrices directly. Mesh interfaces keep independent one-sided values through an optional explicit cell owner.

nodal_field(name: str, mesh: TriangleMesh | TetraMesh | CartesianMacroMesh, degree: int, *, components: int = 1, discontinuous: bool = False, reconstruction: Any = None, trace_reconstruction: Any = None, trace_dofs: Any = None, offset: Any = None) -> FieldDefinition

Declare a scalar/vector Pk/Qk field with portable evaluation and gradients.

Components are interleaved at each node. Continuous simplex coordinates use the package's topological nodal order; DG coordinates concatenate that local order by cell, with degree zero representing one cell value. Cartesian Qk coordinates use x-fastest nodes and the actual archived interval factors used by Qk assembly; no independent tensor basis is regenerated. reconstruction maps the user's complete local unknown vector into this field's coordinates, including block selection and scaling. The executed basis matrix and DOF map are retained by the picklable evaluators and identified by basis_id.

Broken scalar/vector polynomial fields in explicitly ordered modal coordinates.

modal_field(name: str, mesh: Any, powers: tuple[tuple[int, ...], ...], *, convention: str = 'monomial', components: int = 1, reconstruction: Any = None) -> FieldDefinition

Declare cellwise polynomial coefficients with an explicit exponent order.

monomial uses reference coordinate powers. legendre uses conventional L_j(2*x-1), normalized to one at x=1, independently in each direction. This covers complete simplex, tensor Qk and restricted complete tensor Pk spaces through the supplied multiindices. Coefficients interleave physical components within a polynomial, then concatenate cells. The actual interval coordinate matrix is archived and reused for evaluation and physical gradients.

Portable physical Piola fields with literal executed polynomial and moment matrices.

ReferenceVectorBasis dataclass

Executed reference values/divergence before orientation and physical Piola.

Native mode archives Basix's coefficient_matrix and the actual moment map. Monomial mode archives explicitly ordered component powers and the actual moment matrix. Tensor Bernstein mode additionally retains its executed triangle and interval matrices for each vector component. Evaluation never recomputes a dual basis or arbitrarily oriented nullspace.

digest: str property

Identify every executed polynomial/moment matrix and its coordinate convention.

__post_init__() -> None

Archive literal arrays through construction and pickle reconstruction.

__reduce__() -> Any

Replay archived arrays instead of native interpolation and moment solves.

__call__(points: FloatArray) -> tuple[FloatArray, FloatArray]

Return (q,basis,reference_component) values and (q,basis) divergence.

piola_field(name: str, mesh: Any, *, reference_basis: ReferenceVectorBasis, dofs: Any, orientation: Any, components: int = 1, reconstruction: Any = None, divergence: bool = False, trace_reconstruction: Any = None, trace_dofs: Any = None, offset: Any = None) -> FieldDefinition

Declare a physical vector/tensor from explicitly supplied moment and Piola maps.

The reference basis, global DOF map and per-cell orientation matrix are independent mathematical inputs. Components interleave per scalar H(div) DOF; one component returns a physical vector, d components row-wise tensor. divergence=True declares its scalar/vector divergence using the very same archived basis and maps. No continuity or physical PDE is inferred from successful algebraic evaluation. Names such as flux/stress are caller data.

hdiv_field(name: str, mesh: Any, family: Any = 'RT', *, degree: int = 0, enrichment: int = 0, components: int = 1, reconstruction: Any = None, divergence: bool = False) -> FieldDefinition

Declare a built-in moment H(div) field in its actual executed coordinate basis.

family is BDMFamily, HDiv3DFamily, or 'RT', 'tensor-RT', 'mapped-RT'. Family-specific degree conventions remain those of their volume owners. Literal native/moment/orientation matrices are archived once, and replay uses them for both physical evaluation and divergence. Pressure and broken kinematic polynomial fields have separate scalar/nodal definitions.

Physical trace integration and projection independent of a PDE or hybrid method.

trace_quadrature(space: Any, face: int, *, order: int = 6) -> tuple[IntArray, FloatArray, FloatArray, FloatArray]

Return face coordinates, physical points/weights and the executed value basis.

Shapes are (n,), (q,dimension), (q,), (q,n,value_components). Basis values use canonical face coordinates; no incidence sign is inserted. The order floor integrates each polynomial Gram matrix exactly. Nonpolynomial data can require a larger order. Scalar edge/triangle, Cartesian-component and tangential wrappers retain their basis owner's coefficients and continuity. A custom space may implement trace_quadrature(face, order=...) with this same contract, including vector or tensor values flattened by component.

trace_linear_form(space: Any, value: Any, *, faces: Iterable[int] | None = None, order: int = 6) -> FloatArray

Assemble integral(value dot test) in the declared global trace basis.

value is a scalar/vector constant or a callable evaluated at physical points with output (q,value_components). A scalar constant broadcasts to every component; scalar fields may return (q,). Signs and selected physical boundaries are caller choices. Shared coordinates accumulate contributions.

trace_bilinear_form(space: Any, coefficient: Any = 1.0, *, faces: Iterable[int] | None = None, order: int = 6) -> Any

Assemble integral(test dot coefficient trial) as a sparse trace operator.

A scalar coefficient gives an isotropic scalar/vector pairing. A constant or physical-point callable may instead supply a matrix in value components, with shape (m,m) or (q,m,m). This supports mass, reaction, Robin, impedance and tensor pairings without choosing a physical boundary convention.

project_trace(space: Any, value: Any, *, faces: Iterable[int] | None = None, order: int = 6) -> dict[int, float]

L2-project physical values onto selected trace coefficients with shared nodes.

Solve the executed Gram matrix on exactly the union of selected face DOFs. The returned coefficient map can prescribe interface values in any problem. Normal densities, Robin pseudo-tractions and physical scalar/vector traces remain distinct choices of the user-supplied value and space.

trace_boundary_data(space: Any, value: Any, prescribed: Mapping[int, Any] | None = None, *, order: int = 6) -> tuple[FloatArray, dict[int, float]]

Integrate exterior value data and project explicitly prescribed trace data.

Faces in prescribed supply fixed trace coefficients; every other exterior face supplies the unsigned functional integral(value dot test). No PDE, outward sign or physical variable is inferred. Projection is joint across prescribed faces when their coefficient coordinates are shared. Face-dependent prescribed functions use physical scalar/vector values in the owner's actual basis. Interior faces cannot be prescribed here.

Unsigned physical pairings between independently declared trace spaces.

evaluate_trace_basis(space: Any, face: int, points: FloatArray) -> tuple[IntArray, FloatArray]

Evaluate literal face coordinates as physical scalar/vector basis values.

Return shapes (n,) and (q,n,components), with no incidence normal sign. Cartesian-component and tangential wrappers retain their physical frame. Custom owners may supply trace_evaluation(face, points) with this contract. Points are physical; independently subdivided basis owners locate their own pieces. No projection, common-node identification or continuity is inferred.

interface_pairing(test_space: Any, trial_space: Any, cell: int, *, order: int = 6) -> FloatArray

Integrate unsigned test dot trial on a macrocell's complete boundary.

Rows and columns follow each owner's literal cell_dofs(cell). Edge integration uses the union of independent breakpoints. Built-in triangular partitions must be nested, including explicitly declared nonuniform refinements. Their containment is checked before the finer owner supplies the common rule. Nonnested 3D partitions require a custom common_trace_quadrature(other,face,order) owner returning physical points/weights resolving both bases. Physical components are contracted explicitly. Outward signs, coefficients and PDE choices are caller data; this same mass pairing supports diffusion, elasticity and curl traces.

Scalar boundary quadrature and assembly independent of a physical formulation.

The rules expose physical points, positive surface weights, outward normals and an executed basis/map. Coefficients may depend on points and normals, allowing Robin, impedance, reaction and weighted boundary mass terms to use one owner.

BoundaryRule dataclass

One physical facet rule in the supplied scalar coefficient convention.

basis has axes (point, local_dof), dofs selects the global local-field coordinates, and normal points out of the domain whose boundary is integrated. A caller can also supply rules for a different mesh or basis.

boundary_quadrature(mesh: TriangleMesh | TetraMesh | CartesianMacroMesh, degree: int, *, order: int = 6) -> tuple[int, tuple[BoundaryRule, ...]]

Tabulate an exterior scalar Pk/Qk nodal boundary with positive physical weights.

order counts one-dimensional Gauss points, with the degree+1 floor for exact polynomial mass integration. Only exterior facets are returned. Mesh normals use their canonical outward boundary orientation. The returned size is the complete field dimension, including nodes with zero boundary trace.

boundary_bilinear(size: int, rules: tuple[BoundaryRule, ...], coefficient: Any = 1.0) -> sparse.csc_matrix

Assemble integral_boundary coefficient*u*v in explicitly supplied coordinates.

A callback has signature coefficient(points, outward_normal) and returns scalar point values. Negative coefficients are allowed. Rules can represent a selected boundary subset or a custom basis; they define all geometry and orientation. No essential condition, PDE or solver is selected.

boundary_functional(size: int, rules: tuple[BoundaryRule, ...], datum: Any = 0.0) -> FloatArray

Assemble integral_boundary datum*v with the same coordinate/quadrature contract.

A callback has signature datum(points, outward_normal). This operation interprets neither the data as a physical flux nor its sign in an equation; both conventions belong to the formulation written by the user.

Physical normal and representation checks for explicitly declared trace fields.

normal_component_zero(samples: Any, normal: Any, *, rtol: float = 64 * np.finfo(float).eps, pointwise: bool = False) -> bool

Check sampled normal components against their declared cancellation scale.

samples has shape (points, dimension), and normal has shape (dimension,). The comparison is abs(samples @ normal) <= rtol * scale. The scale is the maximum uncancelled absolute dot product, or its value at each point when pointwise=True. No absolute amplitude floor is added. This algebraic check does not certify a callback between sampling points.

boundary_tangent_2d(space: SkeletonSpace, vector: Any, order: int = 8, faces: Iterable[int] | None = None) -> bool

Check sampled tangency on each declared edge partition in two dimensions.

The owner's partition quadrature uses order at least eight. Each face uses the maximum uncancelled dot-product scale and relative tolerance 64*eps, without an absolute amplitude floor. This is a sampled compatibility check for a constant local kernel, not a continuum certificate for a callback.

boundary_tangent_3d(space: TriangularSkeleton | PolygonalSkeleton3D, vector: Any, order: int = 8, *, faces: Iterable[int] | None = None) -> bool

Check tangency on the canonical triangulation of original 3D faces.

Triangular faces use one positive rule of order at least eight. Polygonal faces integrate through their geometry owner's triangles. Physical points use anchored affine evaluation to preserve coordinates constant on a face. Each face uses normal_component_zero with its maximum uncancelled scale and relative tolerance 64 * eps. This preserves the declared sampled compatibility convention; no unsampled continuum condition is certified.

boundary_normal_zero(space: Any, vector: Any, *, order: int = 8, faces: Iterable[int] | None = None) -> bool

Check sampled vector dot normal zero relative to its uncancelled terms.

Exterior face normals are the mesh's canonical outward normals. Quadrature resolves declared trace partitions. This is a floating-point compatibility check, not a certified continuum tangency statement for arbitrary callbacks. Exact zero data remain zero without absolute amplitude thresholds.

trace_represents(space: Any, face: int, value: Any, *, order: int = 8) -> bool

Check a scalar/vector field belongs to the declared trace at integration nodes.

A weighted least-squares solve uses the literal physical value basis and positive quadrature, avoiding a squared Gram condition number. Unresolved basis rank returns false. Residuals are checked relative to uncancelled input and projected magnitudes. The declaration does not certify unsampled values. This operation is independent of a PDE, material, gauge and solver family.

Explicit equations and shared hierarchy

User-defined local and global variational blocks, independent of a physical model.

Rows are test coordinates and columns are trial coordinates. The two equations are A u + B lambda = f and C u + D lambda = g. In particular, the sign of C is declared by the user rather than inferred from a method or equation. UFL forms are assembled by an optional native adapter; arrays are already assembled forms. Neither representation selects a PDE or a stabilization.

Equation dataclass

A bilinear form a and linear right-hand side L in declared coordinates.

A literal zero denotes a zero form with dimensions supplied by the problem. A UFL form keeps its trial/test spaces, mesh, coefficients and quadrature. An array or sparse matrix denotes an already assembled form. Global forms act on trace coordinates followed by the ordered retained cell coordinates.

LinearForms dataclass

Linear forms representing columns B or rows C in an explicit trace basis.

Every entry represents one trace basis function, including its orientation and physical units. An empty sequence requires the problem's dimensions. This representation avoids requiring cross-mesh UFL assembly: each linear form lives on its local integration mesh. It also accepts numerical vectors.

__post_init__() -> None

Freeze the sequence and reject an ambiguous pairing direction.

FormCompiler

Bases: Protocol

Assemble a real rank-one or rank-two form without changing its coordinates.

__call__(form: Any, shape: tuple[int, ...] | None = None) -> Any

Return an owned coefficient array or sparse operator of the declared shape.

LocalEquations dataclass

Define both local equations and their contribution to the global balance.

a, L, b express a(u,v)+b(lambda,v)=L(v). c, d, g express c(u,mu)+d(lambda,mu)=g(mu). dofs and test_dofs map the local trial/test trace bases to the global coordinates. They can differ in order and length. Basis functions include normal orientation; their signs are not inferred. moments and test_moments are column pairings with explicit retained trial/test modes. No kernel is discovered automatically and no stability or consistency theorem is implied.

Scalar, vector and mixed unknowns use this same contract. A method with several local fields can use a mixed UFL space or a block matrix. Native objects must remain worker-owned, with optional reuse between provider invocations and explicit close() cleanup at worker shutdown, when using spawn workers; metadata must then contain only picklable data. A 0-by-0 a declares no eliminated local coordinates: d and g then contribute the complete face equation without a local factorization.

__post_init__() -> None

Freeze explicit maps without assembling forms or inspecting native spaces.

CompiledLocalEquations dataclass

Numerical local contract and its uncondensed direct global terms.

The signed balance C is stored as test_coupling=-C.T in the shared condensation owner. matrix and load are D and g in the union trace numbering, and are added once during ordered global reduction.

columns(*forms: Any) -> LinearForms

Declare one local linear form for each global trial coordinate.

rows(*forms: Any) -> LinearForms

Declare one local linear form for each global test coordinate.

compile_form(form: Any, shape: tuple[int, ...] | None = None) -> Any

Assemble arrays, linear pairings or native UFL without importing unused FEM.

Sparse matrices remain sparse. Literal bilinear zeros use CSC storage; literal linear zeros use dense vectors. Finite wider real array precision is retained. A nonzero scalar is not silently broadcast to a vector or operator. UFL assembly currently requires a real serial DOLFINx mesh. Other backends can supply a callable compiler with this same form/shape contract.

compile_local_equations(equations: LocalEquations, compiler: FormCompiler = compile_form) -> CompiledLocalEquations

Compile four variational blocks and reuse the checked local condensation contract.

Trial and test trace maps embed into their ordered union. This permits rectangular B/C/D pairings while the complete global problem remains square. Retained compatibility equations preserve their explicit test basis and physical moments. Already assembled matrices do not pass through UFL.

Assembly, solution and recursive reconstruction of user-defined hybrid equations.

All numerical elimination delegates to the common constrained condensation and ordered contribution owners. A recursive local operator is itself another MultiscaleProblem; the same contracts apply at every scale. No physical model or method name selects the operator.

MultiscaleProblem dataclass

Bases: Generic[Item]

A global variational equation coupled to independently declared local equations.

global_equation supplies additional bilinear and linear global forms; each local provider supplies its own explicit C u + D lambda = g balance. Global coordinates are trace first, then each cell's retained modes in item order. coarse_sizes declares those counts, including zero. A global form can contain facet terms between different cells, Robin terms or other user-defined couplings. Its full reduced dimensions must match.

The provider runs once per item, inside the selected worker, and returns LocalEquations or a compiled equation record. Scalar/vector/mixed fields and recursive local operators share this API. Boundary coefficients and physical gauge rows are explicit; the package does not choose them from a PDE name.

__post_init__() -> None

Validate coordinates and form/compiler contracts without consuming the items.

from_global(equation: Equation, size: int, *, fixed: Mapping[int, float] | None = None, constraints: tuple[tuple[FloatArray, float], ...] = (), compiler: FormCompiler = compile_form) -> MultiscaleProblem[Any] classmethod

Declare a global equation with no local fields or condensation.

All size coordinates belong to the stated global trial/test basis. Essential coefficient values and physical integral constraints retain the ordinary solve contract. This supports classical reference forms or any explicitly assembled global operator through the same API.

NestedEquations dataclass

Declare a child hybrid problem with a parent boundary-trace constraint.

If the child's reduced equation is A x=f, E selects its boundary trace coordinates and M is trace_map, the local equations are -A x + E mu=-f and E.T x-M lambda=0. The global pairing is +M.T mu in the declared DSL balance C. The original physical moment is -M.T mu; negating that equation matches the shared reduced saddle convention. The local field contains child coordinates followed by boundary reactions. This is a generic equality-constrained operator, independent of a physical PDE. Explicit kernel/left_kernel are in child coordinates; their boundary reactions are lifted by the shared nesting owner.

Children leave boundary data and physical gauges to their parent. Supplied moments have one row per child coordinate and are extended by zero reaction weights. M includes every restriction and orientation; it is not inferred.

MultiscaleSolution dataclass

Bases: HybridSolution

Hybrid coefficients and recursively recovered child solutions.

children[cell] is None for a leaf local discretization. Otherwise its coordinates are exactly that cell's reconstructed field; its own fields and children recover the finer scales without another global solve.

field(name: str, *, recursive: bool = False) -> Any

Recover named views; recursive=True descends into undeclared child cells.

local_trace(index: int, *, test: bool = False) -> FloatArray

Delegate trace recovery in the executed local basis, in assembly item order.

MultiscaleSystem

Bases: HybridSystem

Assembled variational equations with application data and recursive reconstruction.

Assembly stores numerical operators, executed local bases and child systems, without retaining native UFL meshes/forms. The inherited physical moment utility can build gauge rows from local field integrals before solving.

trace_bindings: tuple[Any, ...] property

Expose executed interface maps for persistence and custom-basis replay.

with_rhs(rhs: Any, *, load_scale: Any = None) -> MultiscaleSystem

Delegate a compatible load update to :func:with_global_load.

solve(**options: Any) -> MultiscaleSolution

Delegate solution and hierarchical recovery to :func:solve_multiscale_system.

reconstruct(coordinates: Any, *, rhs: Any = None) -> MultiscaleSolution

Delegate recovery of executed coordinates to :func:reconstruct_multiscale.

with_global_load(system: MultiscaleSystem, rhs: Any, *, load_scale: Any = None) -> MultiscaleSystem

Reuse executed local responses for a new global balance right-hand side.

Retained rows enforce the existing local source's compatibility and must remain unchanged. Changing a volume source requires new local equations or explicitly cached local factors. Child operators and basis matrices are reused literally, and all local fields remain reconstructible.

solve_multiscale_system(system: MultiscaleSystem, **options: Any) -> MultiscaleSolution

Solve the declared boundary/gauge problem and recover every finer field.

Keyword overrides use the common HybridSystem solver contract. Omitting fixed and constraints uses those declared by MultiscaleProblem. Numerical criteria, original equations and basis coordinates are the same as in the shared hybrid solver.

reconstruct_multiscale(system: MultiscaleSystem, coordinates: Any, *, rhs: Any = None) -> MultiscaleSolution

Recover finer fields from executed hybrid coordinates without solving again.

rhs is the actual load at this scale, including parent coupling. This method checks the represented global equations before descending. Parent loads must not change retained leaf compatibility equations.

leaf_moment(system: MultiscaleSystem, weights: Callable[[Any], Any]) -> tuple[FloatArray, float]

Represent a physical finest-scale moment as row @ coordinates + offset.

weights(metadata) supplies one integration vector in each leaf's executed coefficient basis. Vector fields use that basis's component ordering. Child source offsets and boundary-reaction padding are included recursively. A physical gauge with prescribed moment m uses (row, m-offset) as its global constraint; no field is solved or sampled.

with_global_equation(system: MultiscaleSystem, equation: Equation, *, compiler: FormCompiler = compile_form) -> MultiscaleSystem

Add global forms without rebuilding any local operator or executed basis.

The forms act on the existing reduced coordinates. They may depend on the executed local response maps, for example a jump of reconstructed fields. This operation adds both the declared operator and load once, preserving child systems and coefficient interpretation. It returns a new system; it does not infer jump signs, penalties or stability conditions.

assemble(problem: MultiscaleProblem[Item], *, execution: ExecutionConfig = _DEFAULT_EXECUTION, solvers: SolverConfig = _DEFAULT_SOLVERS) -> MultiscaleSystem

Assemble user-defined local/global equations in serial or bounded parallel batches.

Workers build and eliminate local equations. Only the coordinator sums shared global indices, in item order. The additional global equation is assembled once after those contributions. Recursive children run serially within their owning worker, so one resource budget controls the hierarchy. Live optional native resources never cross the worker boundary.

solve(problem: MultiscaleProblem[Item] | MultiscaleSystem, *, execution: ExecutionConfig = _DEFAULT_EXECUTION, solvers: SolverConfig = _DEFAULT_SOLVERS, **options: Any) -> MultiscaleSolution

Solve a problem description or an already assembled variational hierarchy.

execution and solvers configure assembly of a problem description. An assembled system retains those settings; keyword options override its global solve, boundary coefficients and physical constraints through the same owner.

Monolithic assembly of explicitly declared leaf A/B/C/D coefficient equations.

OriginalBlocks dataclass

Original local fields followed by global trace coordinates, without elimination.

local_offsets delimits each original local vector. trace_offset starts the global trace vector. Prescribed trace indices must therefore be shifted by this offset before a caller applies essential data. No gauge, symmetry, invertibility or physical interpretation is inferred.

assemble_original_blocks(system: MultiscaleSystem) -> OriginalBlocks

Assemble the declared uncondensed equations without invoking a physical solver.

The local rows are A u+B lambda=L. Global rows sum C u+D lambda=g and the user's additional global Equation. The signed C is taken literally from the compiled test coupling, independently of B. Matrix contributions are scattered in original local/global coordinate order. The source-dependent rows come from the supplied system, including updates made through OfflineMultiscaleSystem.

This operation currently accepts leaf problems with no retained modes. Recursive hierarchies and retained-mode gauges require a separately declared original layout and are rejected instead of being flattened. It builds a larger sparse operator for direct comparisons; it does not establish stability or constitute an independent physical discretization.

Operator-independent reconstruction from physical coefficient moments.

energy_reconstruction(matrix: Any, moments: FloatArray, refinement_precision: Literal['double', 'extended'] = 'double') -> tuple[FloatArray, FloatArray]

Construct R with A R + C mu = 0 and C.T R = I.

matrix is the finite real square coefficient operator A. moments contains a nonempty set of independent real moment columns C in the original coefficient basis; they retain their physical units. The augmented A/C operator must be invertible. Moment independence alone does not prove that condition, stability of a discretization, or a coercivity estimate.

Return the reconstruction R and the represented reduction R.T A R. For a symmetric positive semidefinite A coercive on the moment-zero space, R minimizes its energy under the declared moments. For a nonsymmetric A, the same original saddle defines a constrained Galerkin reconstruction; its reduction keeps the full bilinear operator, including antisymmetry. Neither operator nor reduction is projected onto its symmetric part. Wider arithmetic follows solve_moment_system without changing rows, moment units or the original 1e-10 residual criterion.

Shared numerical operations and fixed hybrid contracts

Explicit local forms and the scoped algebraic global hybrid form.

Local expressions are opaque to the portable core. A supplied compiler owns their finite-element interpretation; the global form uses the existing hybrid trace/retained coordinates and does not compile arbitrary skeleton UFL forms.

LocalForm dataclass

Describe a(u,v) + sum_j lambda_j b_j(v) = L(v) on one cell.

a and L are compiler-owned expressions, for example UFL forms. trace_forms are the signed linear forms b_j in exactly the order of distinct nonnegative global trace_dofs. Outward/global-normal signs, spaces, quadrature and essential boundary elimination are explicit inputs. No trace orientation or boundary convention is inferred.

kernel declares literal nullspace coefficients; coarse_basis instead retains arbitrary modes without a nullspace claim. They are mutually exclusive. moment_forms has one linear physical moment per retained mode; omitting it requests the compiler's coefficient-orthogonality default. The assembled LocalProblem checks operator/nullspace compatibility and moment pairing. Those algebraic checks do not establish inf-sup stability. Forms are stored as tuples and numerical maps/bases are copied read-only.

__post_init__() -> None

Freeze the declared order and validate maps and literal mode arrays.

GlobalForm dataclass

Declare the sum of local condensed hybrid forms in explicit coordinates.

Trace coordinates occupy [0, trace_size). Retained coordinates follow in cell order with widths coarse_sizes. For each cell, its response's global_contribution supplies the bilinear block and linear load; the shared hybrid assembler sums these records on their declared indices. An empty coarse_sizes tuple declares a global form with no local cells; its complete operator must be supplied explicitly by its owning formulation. This is an algebraic global form, not an arbitrary UFL skeleton compiler.

boundary_load is a trace vector subtracted from the assembled RHS with the HybridSystem boundary convention. fixed_trace prescribes trace coefficients, for example oriented normal fluxes. Each constraints entry is (row, target) imposing row.T @ (lambda,c) = target in the full reduced coordinates. Supply physical pressure/rigid-motion gauges; neither a gauge nor a stability condition is inferred. Local field moments can be converted to these rows by HybridSystem.mean_constraint. Wider floating coefficients are preserved; records are copied on creation.

__post_init__() -> None

Validate the ordered coordinate partition, loads and physical rows.

compile_local_forms(forms: LocalForm, compiler: Callable[[LocalForm], LocalProblem]) -> LocalProblem

Compile one record while preserving its trace map and literal basis.

A compiler is an ordinary callable; no inheritance or optional backend is required. It must return a validated LocalProblem. The caller owns any compiler-native resources and must release them before returning arrays. Numerical operator and integration accuracy remain the compiler's concern. A changed map, undeclared mode or equivalent rotated basis is rejected.

Ordered local providers and assembly of explicitly declared hybrid forms.

Workers own local construction and solves. Only the coordinator accumulates shared trace entries, in macrocell order and without intermediate batch sums. The module requires no optional FEM, MPI or accelerator imports.

LocalLinearSolver

Bases: Protocol

Solve all columns of one constrained local system in its declared basis.

A callable receives owned copies of the augmented sparse operator and source/trace/retained-mode right-hand sides. Return a finite real array of the same shape as rhs, including moment multipliers. External codes or learned response models can implement this contract without inheritance. They own their native resources and may implement close() to release resident solver state when local execution ends. The original constrained operator has an independent numerical rank check, and each column is checked before decoding a local response.

__call__(matrix: Any, rhs: FloatArray) -> FloatArray

Return coefficients satisfying matrix @ result = rhs.

SolverConfig dataclass

Choose local and global algebra without changing a variational form.

Named solvers use the existing PyMHM adapters. A callable local_solver implements :class:LocalLinearSolver and is checked at the unchanged relative residual tolerance 1e-10, separately for each right-hand side. Explicit precision choices affect native correction arithmetic, not the operator, quadrature, gauge or acceptance criterion. A custom solver owns its precision policy; use the default local precision setting with it.

__post_init__() -> None

Reject ambiguous solver and precision settings before executing work.

HybridProblem dataclass

Bases: Generic[Item]

Combine a global hybrid form with an independent local provider.

items supplies one specification per macrocell, in the order declared by global_form.coarse_sizes. The provider returns LocalProblem or LocalAssembly with oriented trace maps, actual bases and physical moments. It is called exactly once per item per assembly. An iterator is consumed once; use a reusable iterable for repeated assemblies. Native resources remain worker-owned and must not cross a process boundary. A provider may retain reusable workspaces between items and implement close() for worker-lifetime cleanup. Threaded providers must keep mutable native workspaces separate across concurrent calls. Providers and items must be picklable for the spawn backend. contribution can supply explicit direct global terms in addition to a cell's Schur block. By default only the coordinator invokes a custom callback, in cell order. contribution_execution="worker" declares an independent, picklable callback that runs after the local solve on that worker; it must not write shared state or require coordinator resources. The default Schur contribution always runs on its local worker. Only the ordered coordinator reduction writes shared global entries. Set require_local_trace_coverage=False when a separate declared global form owns the remaining trace coordinates; retained intervals are checked.

__post_init__() -> None

Check the global form and provider without consuming local items.

assemble_hybrid(problem: HybridProblem[Item], *, execution: ExecutionConfig = _DEFAULT_EXECUTION, solvers: SolverConfig = _DEFAULT_SOLVERS) -> HybridSystem

Assemble a hybrid form using ordered serial work or bounded parallel jobs.

Serial execution solves and adds each cell before requesting the next item. Thread/process execution completes one batch before consuming the next; contributions are reduced individually in input order by the coordinator. With execution.pipeline=True, a rolling window replaces batch barriers; other workers continue while the coordinator scatters the next ordered block. The default Schur reduction and explicitly worker-owned callbacks run with the local solve, before that handoff. Shared faces therefore have no concurrent writes or batch-dependent sums. The declared trace and retained layout is checked against every response. Local responses are retained for reconstruction and physical mean rows; this API bounds in-flight work, not total response/global matrix storage.

solve_hybrid(problem: HybridProblem[Item], *, execution: ExecutionConfig = _DEFAULT_EXECUTION, solvers: SolverConfig = _DEFAULT_SOLVERS) -> HybridSolution

Assemble, solve and reconstruct with the declared boundaries and gauges.

fixed_trace prescribes oriented trace coefficients. constraints holds physical rows and targets in global trace/coarse coordinates. The unchanged global solver verifies the original equations after eliminating prescribed values and enforcing gauges. To derive mean rows from local physical weights, call :func:assemble_hybrid and system.mean_constraint before system.solve instead.

Validated local equations, cached responses and coefficient data contracts.

LocalProblem dataclass

A local variational problem A u + B lambda = f.

kernel contains a basis Z of the right nullspace of A; left_kernel defaults to the same basis. Distinct test and trial kernels are allowed. Alternatively, coarse_basis retains an arbitrary basis Z, including nearly null modes, in the global system. It is mutually exclusive with kernel; no nullspace claim is made for this general basis. constraints is C in C.T w = 0 and must pair nonsingularly with the retained Z; by default C=Z. Physical L2 orthogonality uses C=M Z. Signed normal orientations must already be included in B. trace_dofs maps its columns to the global skeleton. test_coupling is C in the global equation C.T u = g and defaults to B. test_basis and test_constraints specify the corresponding test-side retained modes and compatibility lifting. This permits Petrov formulations with different left/right kernels and independent B/C couplings. A declared kernel is tested at the matrix scale. Its nonzero floating-point actions A Z and A.T W are retained in the original equations through E=Z-R A Z; only exactly vanishing represented actions use the kernel shortcut. No boundary pressure constraints are applied locally.

__init__(matrix: Any, coupling: Any, load: Any, trace_dofs: Any, kernel: Any = None, constraints: Any = None, *, coarse_basis: Any = None, test_coupling: Any = None, left_kernel: Any = None, test_basis: Any = None, test_constraints: Any = None) -> None

Check finite square blocks and the declared kernel, including a zero local space.

condense(solver: str = 'scipy', *, refinement_precision: Literal['double', 'extended'] = 'double') -> LocalResponse

Delegate constrained source/trace elimination to :func:condense_local.

condensation_system() -> tuple[Any, FloatArray]

Return the constrained operator and source/trace/retained right sides.

condensation_matrix() -> Any

Return the scaled constraint operator without allocating response columns.

reconstruct(trace: Any, coarse: Any, *, retained_basis: Any = None, solver: str = 'scipy', factorization: LinearFactorization | None = None, refinement_precision: Literal['double', 'extended'] = 'double') -> FloatArray

Reconstruct physical fields through :func:reconstruct_local.

response_from_solution(solution: Any) -> LocalResponse

Decode finite real responses, preserving explicitly retained wider precision.

condensed_load(source_response: Any, *, load: Any = None, corrected_retained: bool | None = None) -> FloatArray

Project executed source responses into the original reduced Petrov right side.

with_load(load: Any, *, preserve_precision: bool = False) -> LocalProblem

Copy this contract with a different finite source vector.

The default retains the canonical double input convention. Explicit preserve_precision=True keeps real wider defect-source digits for original-equation refinement, including augmented/reduced right sides. Matrices, kernels, moment constraints and trace orientation are unchanged.

LocalAssembly dataclass

A local variational problem and optional application reconstruction data.

metadata can hold a NumPy mesh, DOF map or other application data that should return with the condensed response. Process workers require picklable inputs and metadata: do not capture or return live FEniCS, PETSc, MPI or CUDA objects. Construct native resources inside the worker, extract numerical arrays into LocalProblem, and release resources inside that worker.

__post_init__() -> None

Reject a factory payload that does not contain a validated local problem.

LocalResponse dataclass

Cached source/trace lifts and reconstructed retained modes of one cell.

retained_basis: FloatArray property

Return Z for an exact represented kernel, otherwise E=Z-R A Z.

reconstruct(trace: FloatArray, coarse: FloatArray) -> FloatArray

Evaluate cached lifts and retained modes through :func:reconstruct_response.

kernel_roundoff_bound(field: FloatArray) -> FloatArray

Bound declared-kernel evaluation through the shared reconstruction owner.

global_load() -> FloatArray

Return source-dependent reduced loads without rebuilding matrix blocks.

global_contribution(coarse_dofs: Any) -> tuple[IntArray, FloatArray, FloatArray]

Return the oriented cell block and load; see local_global_contribution.

HybridSolution dataclass

Skeleton coefficients, retained coarse amplitudes and reconstructed local fields.

residual is the verified relative compatibility defect after any componentwise declared-kernel evaluation allowance. Global solves also preserve its unadjusted value in raw_residual and its absolute norm in raw_residual_norm. Manually constructed solutions may omit these records.

Constrained local elimination and portable local assembly workers.

condense_local(problem: LocalProblem, solver: str = 'scipy', *, refinement_precision: Literal['double', 'extended'] = 'double') -> LocalResponse

Factor once for source, trace and retained responses, with explicit precision.

problem declares A u+B lambda=f and the oriented local trace columns. The response stores R f, R B and the executed retained basis E=Z-R A Z; an exactly represented kernel uses Z directly. Constraints, Petrov test data and solver residual criteria remain those of the supplied problem. Direct factors can retain corrected lifts in extended precision without changing their residual criterion. AMG local projection currently accepts only the default double-precision accumulation. AMG refines the projected original equations to relative tolerance 1e-12 while each pinned correction solve keeps its 1e-10 residual criterion. This reserves two decimal orders for source/trace combinations; their defects add with the reconstruction coefficients. This reserve does not replace the physical reconstructed-field residual criterion, checked independently. AmgX prepares the pinned operator once, reusing its hierarchy for every column and projected correction, and closes it on success or failure. Compatibility and gauge moments use the shared accurate accumulation, including contractions across retained modes. All four permitted corrections are checked against the original projected equations before acceptance.

local_condensation_system(problem: LocalProblem) -> tuple[Any, FloatArray]

Return the constrained operator and source/trace/retained-mode right sides.

The augmented operator is [[A, C_left], [C_right.T, 0]] with constraint columns scaled down only when they exceed the operator's largest entry. This change of auxiliary coordinates prevents constraints from dominating small-unit operators without magnifying constraints in mixed high-contrast saddles. It leaves the physical complement unchanged. Its inverse is never explicitly formed. Backends may factor this matrix once, retain the factorization, and solve new source columns online.

local_condensation_matrix(problem: LocalProblem) -> Any

Return the augmented operator without allocating response right-hand sides.

This is the same scaled left/right constraint operator returned by condensation_system. It permits one caller-owned factorization for several direct reconstructions, including nonsymmetric Petrov data.

local_response_from_solution(problem: LocalProblem, solution: Any) -> LocalResponse

Decode finite real responses, preserving explicitly retained wider precision.

Physical reconstruction and reduced loads in executed retained coordinates.

reconstruct_local(problem: LocalProblem, trace: Any, coarse: Any, *, retained_basis: Any = None, solver: str = 'scipy', factorization: LinearFactorization | None = None, refinement_precision: Literal['double', 'extended'] = 'double') -> FloatArray

Reconstruct local fields without storing every trace response.

trace uses local coupling-column order and coarse uses retained trial coordinates. Vectors produce one field; matching nonempty column arrays produce several fields with one factorization. With R the volumetric block of the augmented inverse, the field is R(f-B trace-A Z coarse)+Z coarse for an arbitrary retained Z. This includes nonsymmetric Petrov problems and an empty retained basis. retained_basis supplies an executed E=Z-R A Z for persisted coarse coordinates; then the field is R(f-B trace)+E coarse. Its numerical basis and orientation belong to the caller's archived coefficient contract. Only combined source columns are solved. Explicit extended precision preserves forcing and correction digits; backend factors stay double precision. A supplied factorization must match the augmented operator exactly and remains caller-owned. The linear residual uses the combined augmented RHS; full physical compatibility remains a global original-equation check.

local_condensed_load(problem: LocalProblem, source_response: Any, *, load: Any = None, corrected_retained: bool | None = None) -> FloatArray

Project executed source responses to the original reduced Petrov RHS.

source_response is R times the supplied load (the stored source by default). A vector gives one reduced RHS; nonempty column arrays give several. Real wider source/load digits are preserved. The trace part is test_coupling.T R load. General or rounded retained modes use (A.T W).T R load-W.T load; exact represented kernels use -W.T load. corrected_retained can preserve an executed LocalResponse branch. It does not change the augmented operator or test/trial conventions.

reconstruct_response(response: LocalResponse, trace: FloatArray, coarse: FloatArray) -> FloatArray

Recover one local coefficient vector from trace and retained amplitudes.

trace follows the local coupling columns, and coarse follows the executed retained basis columns. The field is source-lifts@trace+E@coarse; wider real operands follow NumPy's existing arithmetic promotion. This cached-response operation takes vectors; reconstruct_local also accepts matching nonempty RHS column arrays.

kernel_roundoff_bound(response: LocalResponse, field: FloatArray) -> FloatArray

Bound floating-point evaluation of retained kernel equations.

For declared left/right kernels only, sum the componentwise row bound gamma_(nnz_i+1) * (abs(A) @ abs(u))_i against abs(W) and include the final test projection's dot-product bound. Here gamma uses the stored double-precision operator, even for wider reconstruction. No operator entry, represented kernel action or physical residual is replaced by zero. General retained modes receive no kernel allowance.

Ordered sparse reduction of local blocks with shared trace indices.

local_global_contribution(response: LocalResponse, coarse_dofs: Any, *, direct_matrix: Any = None, direct_load: Any = None) -> tuple[IntArray, FloatArray, FloatArray]

Return indices, matrix and RHS of this cell's reduced Petrov equations.

This is the common assembly contract for serial, MPI and accelerator backends. The trace rows enforce C.T u=g; retained test rows enforce W.T (A u+B lambda-f)=0. The negative signs preserve the symmetric saddle convention whenever C=B and test/trial data coincide. direct_matrix and direct_load add explicitly declared local D/g terms to the trace rows, in the same local trace ordering. Their defaults are zero. Wider real floating digits are preserved. The same operation therefore owns four-block variational forms in serial, CPU workers and MPI.

assemble_hybrid_contributions(contributions: Iterable[tuple[IntArray, FloatArray, FloatArray]], *, trace_size: int, kernel_offsets: IntArray, boundary_load: Any = None, require_local_trace_coverage: bool = True) -> tuple[Any, FloatArray, FloatArray]

Reduce ordered cell contributions to a CSC operator, RHS and load scale.

kernel_offsets starts at trace_size and bounds each cell's retained global coordinates. Exactly one contribution per interval is consumed in iterable order, without retaining a tuple of contributions. Shared trace entries are added in that order; wider real matrix and load digits survive. The COO matrix entries remain stored until the final CSC conversion. load_scale sums absolute loads before cancellation. The weak Dirichlet boundary_load is subtracted once from the trace RHS and its absolute value added to that scale. No trace orientation, gauge or local basis is inferred. The caller supplies the executed, consistently numbered blocks. require_local_trace_coverage=False permits coordinates belonging only to an explicitly assembled global form; retained cell intervals are always checked. The complete global solver still rejects a singular operator. A single kernel offset and no contributions represent a global-only form; uncovered trace coordinates require require_local_trace_coverage=False.

Global hybrid layout, boundary elimination, gauges and physical compatibility.

HybridSystem

Assemble the global MHM equations independently of the local FEM backend.

The unknown is (lambda, c). Dirichlet traces enter through boundary_load as integrals against oriented skeleton basis functions. Prescribed Neumann flux coefficients are passed to solve as fixed DOFs. Arbitrary linear mean constraints can remove global pressure or rigid-motion nullspaces.

__init__(problems: list[LocalProblem] | tuple[LocalProblem, ...], *, boundary_load: Any = None, local_solver: str = 'scipy', local_refinement_precision: Literal['double', 'extended'] = 'double', backend: Literal['serial', 'thread', 'process'] = 'serial', workers: int | None = None, native_threads: int | None = 1, batch_size: int | None = None) -> None

Condense independent subdomains and assemble their sparse contributions.

from_local_factory(factory: Callable[[Any], LocalProblem | LocalAssembly], items: Iterable[Any], *, boundary_load: Any = None, local_solver: str = 'scipy', local_refinement_precision: Literal['double', 'extended'] = 'double', backend: Literal['serial', 'thread', 'process'] = 'serial', workers: int | None = None, native_threads: int | None = 1, batch_size: int | None = None) -> HybridSystem classmethod

Assemble and condense each independent local problem inside its worker.

factory(item) returns LocalProblem or LocalAssembly. The factory and one condensation execute together, exactly once per item; global assembly and solution remain in the parent. Input order defines response order and local_metadata, including None for a bare problem. Exceptions propagate without a serial fallback.

Serial and thread backends accept closures. The portable spawn backend requires picklable factories, items and metadata, and an executable script protected by if __name__ == '__main__':. Native FEM, PETSc, MPI and CUDA resources must be created and released inside a worker, never sent to or returned from a worker. An optional factory.close() or solver close() runs at the worker's lifetime boundary, including local failures. native_threads defaults to one thread; backend library thread safety is the factory's concern. batch_size limits submitted local jobs while preserving input order. Responses are retained before global assembly because this constructor infers the global layout. For immediate serial contribution with an explicitly declared layout, use assemble_hybrid.

from_responses(responses: Iterable[LocalResponse], *, boundary_load: Any = None, metadata: Iterable[Any] | None = None) -> HybridSystem classmethod

Assemble cached local responses without repeating any factorization.

from_contributions(contributions: Iterable[tuple[IntArray, FloatArray, FloatArray]], *, trace_size: int, coarse_sizes: Iterable[int], boundary_load: Any = None) -> HybridSystem classmethod

Assemble compact cell contributions without retaining local lifts.

Each tuple is the output of LocalResponse.global_contribution with coarse indices numbered after the trace_size skeleton unknowns. coarse_sizes declares their ordered cell partition. solve returns trace/coarse coefficients and an empty field tuple; callers reconstruct one local response at a time using the executed, persisted local basis. Physical mean rows must be supplied explicitly to solve: this compact representation cannot derive them or certify kernel evaluation roundoff.

mean_constraint(local_weights: list[FloatArray] | tuple[FloatArray, ...], value: float = 0.0) -> tuple[FloatArray, float]

Express a prescribed integral of reconstructed fields as r.T x=b.

For incompressible flow, weights integrate only the pressure components; for scalar pure Neumann diffusion, they integrate the pressure field.

with_rhs(rhs: Any, *, load_scale: Any = None) -> HybridSystem

Reuse the executed condensed matrix with a new explicit reduced load.

This compact online system retains the original trace numbering and ordered retained partitions, and returns no reconstructed fields. load_scale is the sum of absolute contributions before cancellation; omitting it uses abs(rhs). Wider real floating loads are preserved. Physical moment rows must be supplied explicitly to solve.

solve(*, solver: str = 'scipy', fixed: dict[int, float] | None = None, constraints: list[tuple[FloatArray, float]] | None = None, factorization: LinearFactorization | None = None, refinement_precision: Literal['double', 'extended'] = 'double', rtol: float | None = None) -> HybridSolution

Solve with boundary elimination and physical checks; see solve_hybrid_system.

solve_hybrid_system(system: HybridSystem, *, solver: str = 'scipy', fixed: dict[int, float] | None = None, constraints: list[tuple[FloatArray, float]] | None = None, factorization: LinearFactorization | None = None, refinement_precision: Literal['double', 'extended'] = 'double', rtol: float | None = None) -> HybridSolution

Solve the saddle system, enforce gauges and verify original equations.

A singular skeleton typically indicates missing gauges or a trace space richer than the local response space. Such systems are never repaired by an undocumented diagonal perturbation or least-squares solution. For a nonsymmetric global operator, the symmetric gauge augmentation requires constraint rows pairing with both left and right nullspaces. A supplied factorization must match the complete constrained matrix exactly. It is reused without transferring ownership, allowing repeated sources and boundary values through an offline/online preparation. Explicit refinement_precision="extended" retains correction digits in trace, coarse coefficients and reconstructed fields. Backend factors remain double precision. rtol selects the global linear-solver relative residual tolerance; None uses 1e-10 for a fresh solve or the tolerance of a supplied factorization. An explicit tolerance must match that prepared factorization. Physical compatibility and local reconstruction criteria remain unchanged. Extended arithmetic requires a wider NumPy long-double type. The original-equation residual uses the absolute local load contributions before assembly, and the net prescribed action before elimination. This accounts for cancellation roundoff without applying a relative tolerance to cancelled large prescribed fluxes or to artificial gauge equations. Prescribed elimination also admits its componentwise floating-point bound gamma_(m+1) * abs(A_fixed) @ abs(x_fixed) per physical row, where m is its number of terms and gamma uses double precision. A gauged declared kernel can additionally use the componentwise evaluation bound of its retained physical rows. This bound is consulted only after the usual compatibility test fails; it changes neither the represented operator nor the raw residual or previously accepted result. General retained modes receive no kernel allowance.

hybrid_mean_constraint(system: HybridSystem, local_weights: list[FloatArray] | tuple[FloatArray, ...], value: float = 0.0) -> tuple[FloatArray, float]

Express a prescribed integral of reconstructed fields as r.T x=b.

For incompressible flow, weights integrate only the pressure components; for scalar pure Neumann diffusion, they integrate the pressure field.

Optional defect correction of the original uncondensed hybrid equations.

HybridRefinement dataclass

Corrected physical fields and full-equation residual history in input units.

residual_norms includes the initial state and every accepted correction. It combines the original local rows, unconstrained weak-trace rows and supplied physical moments. rhs_norm uses these same rows and units; this algebraic norm is not a finite-element error estimate. Persist solution.fields with their executed evaluation basis: trace and coarse coefficients alone cannot replay the added defect-source responses. solution.residual checks free weak-trace and retained-test equations on these corrected fields; artificial gauge rows are excluded from that diagnostic. Gauge multipliers include all accepted corrections.

HybridRefinementCase dataclass

Original reduced system, executed coordinates and physical data of one case.

moments contains local physical weights and original integral targets. Compact systems additionally require moment_rows: the executed reduced rows expressing these same integrals, without subtracting source moments. The number of moments must equal the initial gauge multiplier count.

HybridRefinementLocal dataclass

One original local operator and the retained evaluation basis actually used.

Every streamed case uses the same problem and retained_basis. An optional trace_maps entry (indices, injection) for each case maps its oriented trace coefficients to the columns of problem.coupling. Omission uses problem.trace_dofs with the identity injection. The local test and trial couplings remain those of the original problem.

HybridRefinementStore

Bases: Protocol

Caller-owned bounded storage for executed fields and correction records.

Arrays have shape (local_field_size, number_of_cases) and the selected correction dtype. read_fields/write_fields access the current physical fields. write_record/read_record preserve records named initial (step zero), load (the measured original defect), forcing (the applied defect, zero for already accepted cases), source (its executed augmented inverse response), and correction (the actual physical delta). Replaying fields adds the recorded deltas to initial in step order and the executed dtype. Persist each case's returned coordinate increments and the executed local bases as well. Storage must preserve every digit; each write is read back and verified before the next operation.

read_fields(cell: int) -> FloatArray

Read the current physical fields of one cell.

write_fields(cell: int, fields: FloatArray) -> None

Store the current physical fields without narrowing their precision.

write_record(name: str, step: int, cell: int, values: FloatArray) -> None

Persist one executed array in the declared precision.

read_record(name: str, step: int, cell: int) -> FloatArray

Restore an executed array without recomputing its conditional solve.

HybridStreamRefinement dataclass

One streamed case's original residual history and executed coordinate deltas.

solution.fields is empty; physical fields remain in the supplied store. increments contains the actual trace, retained and gauge corrections in execution order. local_contract_digests covers original operators, data, evaluation bases and oriented trace maps, and is checked on every pass. The residual norm has the same physical-row convention as HybridRefinement.

refine_hybrid_stream(cases: Sequence[HybridRefinementCase], local_factory: Callable[[int], HybridRefinementLocal], store: HybridRefinementStore, *, max_steps: int = 2, min_steps: int = 0, rtol: float = 1e-10, atol: float = 0.0, solver: str = 'scipy', local_solver: str = 'scipy', refinement_precision: Literal['double', 'extended'] = 'double') -> tuple[HybridStreamRefinement, ...]

Correct original hybrid equations cellwise with a fixed condensed operator.

Cases share original local matrices and data, but may have different oriented trace spaces, boundary loads, prescribed coefficients and physical moments. Each bounded pass calls local_factory(cell) in cell order; actual local bases and all operator/data values must stay identical between passes. Local factors solve batched defect sources and are released before the next cell. The executed initial field, measured load, inverse response and actual physical increment are stored; no conditional solve is repeated for replay.

Acceptance uses the original volume rows, free weak trace rows and physical moments with max(atol, rtol*rhs_norm). Artificial gauge equations are excluded. rhs_norm includes prescribed physical trace action. The same reduced matrix is used for corrections, and prescribed trace increments are zero. No tolerance, operator entry or approximation space is modified. Compact solutions and per-case coordinate increments are returned; fields remain in store. Rejected states raise LinearSolveError and must not be published as accepted acquisitions. min_steps optionally requires further original-equation corrections even when the residual is already accepted. It is at most max_steps; zero preserves residual-based stopping. It does not certify field accuracy or change the original residual criterion.

refine_hybrid(system: HybridSystem, solution: HybridSolution, *, boundary_load: Any = None, fixed: dict[int, float] | None = None, moments: Sequence[tuple[Sequence[FloatArray], float]] = (), max_steps: int = 2, min_steps: int = 0, rtol: float = 1e-10, atol: float = 0.0, solver: str = 'scipy', local_solver: str = 'scipy', refinement_precision: Literal['double', 'extended'] = 'double') -> HybridRefinement

Refine A_i u_i+B_i lambda=f_i and the original weak trace equations.

Supply the original boundary_load (zero by default), prescribed trace coefficients fixed, and physical moments as in the initial solve. A moment is (local_weights, target); do not pass reduced gauge rows. Distinct test/trial couplings and multiple retained modes are supported.

Each correction retains the original lifts and condensed operator. Local direct factors solve only the new defect source and are released cellwise; no full uncondensed matrix is assembled. The default solution is untouched unless this function is explicitly called. max_steps=0 checks it only. Acceptance requires the full original residual to be at most max(atol, rtol*rhs_norm). Failure raises LinearSolveError with the achieved residual, including stagnation at a precision floor. Extended correction storage is explicit and requires a wider NumPy type. min_steps optionally requires original-equation corrections after the residual is already accepted, up to max_steps. Zero preserves the default stopping convention. A small residual or a requested correction count alone does not certify physical field accuracy.

Reusable MHM operators for repeated sources and fixed material/geometry.

An offline object keeps the local augmented factorizations, harmonic lifts and one constrained global factorization. Online calls may change source vectors, boundary moments and prescribed trace values, but not matrices, spaces, fixed DOF locations or the definition of physical mean constraints.

LocalFactorCache

Reuse direct factors only for exactly identical constrained local matrices.

Sources, trace couplings and retained responses may differ between cells. Keys contain the canonical CSC shape, index arrays and every numerical byte; no tolerance, geometric equivalence or approximate hashing is used. The cache owns all factors until close and is local to one process/thread of use.

size: int property

Return the number of distinct retained factorizations.

__init__(*, solver: str = 'scipy') -> None

Create an empty explicit-lifetime cache for one direct solver backend.

condense(problem: LocalProblem) -> LocalResponse

Condense one validated local problem using an exact-matrix cache entry.

close() -> None

Release all factors exactly once and invalidate further use.

__enter__() -> LocalFactorCache

Enter an open exact-matrix cache context.

__exit__(exc_type: type[BaseException] | None, exc_value: BaseException | None, traceback: TracebackType | None) -> None

Release cached CPU or accelerator resources on every exit path.

OfflineLocalProblem

Keep a local augmented factorization and source-independent responses.

The input LocalProblem validates the mathematical kernel/basis contract. New loads are expressed in that same test space. This object owns native solver resources and must be closed; concurrent calls are unsupported.

__init__(problem: LocalProblem, *, solver: str = 'scipy') -> None

Factor once and build the harmonic and retained-mode responses.

response(load: Any) -> LocalResponse

Solve only the new source response, reusing all homogeneous lifts.

close() -> None

Release the reusable local factorization exactly once.

__enter__() -> OfflineLocalProblem

Enter an open local-factorization context.

__exit__(exc_type: type[BaseException] | None, exc_value: BaseException | None, traceback: TracebackType | None) -> None

Release native resources even if an online solve fails.

OfflineHybridSystem

Prepare a fixed MHM discretization for multiple source/boundary queries.

moments contains (local_weights, target) pairs, with one physical weight vector per local problem. Their rows are fixed; their right-hand sides account for each new source lifting. Online targets may change prescribed integrals. fixed values can change, but their keys cannot.

This is an algebraic offline/online organization of the supplied method; it does not change that method into MsHHO or project arbitrary sources onto the polynomial source space of a different formulation.

original_residuals: tuple[float, float] | None property

Return the last checked original residual and physical RHS Euclidean norms.

The residual includes original local rows, free weak trace equations and physical moments. The RHS accounts for prescribed trace values. An unchecked or failed online query clears this record. These norms measure algebraic satisfaction, not approximation error or stability.

__init__(problems: Iterable[LocalProblem], *, boundary_load: Any = None, fixed: dict[int, float] | None = None, moments: Sequence[tuple[Sequence[FloatArray], float]] | None = None, local_solver: str = 'scipy', solver: str = 'scipy') -> None

Prepare local lifts and the constrained global factorization once.

solve(loads: Sequence[Any], *, boundary_load: Any = None, fixed: dict[int, float] | None = None, targets: Sequence[float] | None = None, check_original: bool = False) -> HybridSolution

Solve one new source query without rebuilding or refactoring operators.

check_original=True checks the executed solution against the original physical equations with the shared refine_hybrid(max_steps=0) criterion, residual <= 1e-10 * rhs_norm. It performs no corrections, returns the same solution, and records the norms in original_residuals.

solve_many(loads: Iterable[Sequence[Any]], **kwargs: Any) -> tuple[HybridSolution, ...]

Evaluate source cases in order while retaining all offline resources.

close() -> None

Release all local and global native resources; repeated calls are safe.

__enter__() -> OfflineHybridSystem

Enter an open offline/online context.

__exit__(exc_type: type[BaseException] | None, exc_value: BaseException | None, traceback: TracebackType | None) -> None

Release resources after the last query or an exception.

condense_cached(problems: Iterable[LocalProblem], *, solver: str = 'scipy') -> tuple[LocalResponse, ...]

Condense ordered cells using exact-matrix reuse within this single call.

Repeated variational sources with unchanged operators and executed bases.

OfflineMultiscaleSystem

Reuse local factors for arbitrary new volume and interface functionals.

The input system supplies A/B/C/D, retained trial/test modes, coefficient transports and field definitions. Only L, g and the additional global functional may change. Loads use the compiled test-coordinate order; balance_loads uses the compiled union of local trace coordinates. No physical method, sign, source projection or time scheme is selected.

This object owns local factor resources and must be closed. Updated systems preserve their operator and harmonic lifts literally and are independently reconstructible after close. A global factor can be reused through the ordinary solve factorization option. Source-dependent gauge targets must be recomputed from the updated system's mean_constraint or leaf_moment. A recursive hierarchy updates its child sources before assembling its parent; this leaf-source cache does not flatten that tree.

__init__(system: MultiscaleSystem, *, solver: str = 'scipy') -> None

Factor the declared local constrained operators with explicit ownership.

with_loads(local_loads: Sequence[Any], *, global_load: Any = None, balance_loads: Sequence[Any] | None = None) -> MultiscaleSystem

Return a reconstructed-source system without rebuilding any operator.

Omitted global/balance loads retain the original additional functional and g. Explicit global_load replaces that additional functional in the full trace-plus-retained layout, rather than the condensed RHS. Local compatibility contributions are recomputed from the new L.

solve(local_loads: Sequence[Any], *, global_load: Any = None, balance_loads: Sequence[Any] | None = None, **options: Any) -> MultiscaleSolution

Solve an updated source through the existing generic solver contract.

close() -> None

Release owned local factors once; previously returned systems remain usable.

__enter__() -> OfflineMultiscaleSystem

Enter the explicit lifetime of an open factor cache.

__exit__(exc_type: type[BaseException] | None, exc_value: BaseException | None, traceback: TracebackType | None) -> None

Release resources on normal completion or an unsuccessful source query.

Recursive local MHM through a second, explicit trace constraint.

The inner hybrid unknowns and their boundary reactions form the local field of the next level. Inner volume fields are reconstructed only after the outer solve. No inverse, penalty, or implicit boundary gauge is introduced.

NestedSolution dataclass

Recovered inner hybrid coordinates, leaf fields and boundary reactions.

NestedLocalProblem dataclass

One local problem whose interior discretization is another MHM system.

Use problem in the parent :class:HybridSystem and pass the resulting local field to :meth:reconstruct. The trace map sends outer coefficients to selected inner trace coordinates; it includes every orientation and basis transformation. A kernel is declared in inner hybrid coordinates and lifted to the boundary reactions, so retained constants or rigid modes survive each condensation level.

reconstruct(field: Any) -> NestedSolution

Recover all inner local fields and check unconstrained inner equations.

moment(local_weights: list[FloatArray] | tuple[FloatArray, ...]) -> tuple[FloatArray, float]

Express a physical leaf moment as weights @ nested_field + offset.

The offset is essential when the inner source responses have nonzero moments. Subtract it from a parent physical gauge's prescribed value.

nested_trace_map(outer: SkeletonSpace, cell: int, inner: SkeletonSpace) -> tuple[IntArray, FloatArray]

Restrict oriented parent face polynomials exactly to inner boundary faces.

The inner mesh must partition the selected parent cell. Its boundary edges must lie on individual straight parent edges, and every restricted outer basis function must belong to the inner face space. Extra material cuts or smaller inner degrees cannot silently project an unrepresentable trace. Endpoint projections agree with parent breakpoints within the declared coordinate-roundoff envelope. Corrections smaller than half the child edge retain its ordering; genuine interior material cuts remain unrepresentable. Components retain the standard component-interleaved ordering.

nest_hybrid_system(inner: HybridSystem, boundary_dofs: Any, trace_map: Any, trace_dofs: Any, *, kernel: Any = None, left_kernel: Any = None, constraints: Any = None, test_constraints: Any = None) -> NestedLocalProblem

Promote an inner hybrid discretization to a parent local problem.

With A x=b the inner equations, E selecting its boundary trace and M=trace_map, the local equations are [[-A,E],[E.T,0]] [x,mu] + [0,-M] lambda = [-b,0]. Thus E.T x=M lambda and mu is minus the inner potential moment. The parent continuity pairing -M.T mu has the original physical sign.

kernel and left_kernel are arrays with one row per inner hybrid unknown, and the same number of independent columns. Their reaction parts are calculated from A and A.T. They must vanish in the selected trace coordinates and satisfy the interior homogeneous equations. Constraints have the same inner shape and are extended by zero boundary reactions. The caller supplies the physical constant/rigid modes; numerical nullspace detection does not replace this formulation contract. Distinct left and right kernels are supported when the inner method is Petrov--Galerkin.

Exact restriction of prepared local responses to nested skeletal subspaces.

moment_complement(moments: Any, *, pivots: Any = None) -> FloatArray

Construct coordinates annihilating independent declared moments.

moments is a real (n, m) matrix: its columns represent linear functionals on an n-coefficient vector. A vector declares one moment. The returned (n, n - m) matrix Z satisfies moments.T @ Z = 0. Its free-coordinate rows form the identity, in ascending coordinate order; pivot rows enforce the moments. This convention fixes scale and orientation without choosing an arbitrary numerical nullspace basis.

pivots optionally declares the m constrained coordinates. Otherwise the largest absolute entry is selected for one moment, and column-pivoted QR of moments.T selects coordinates for several moments. Near-ties in QR can depend on the numerical platform: archive the executed matrix with persisted coefficients, or supply fixed admissible pivots. No physical kernel or integration measure is inferred from these algebraic moments.

restrict_response(response: LocalResponse, injection: Any, trace_dofs: Any, *, test_injection: Any = None) -> LocalResponse

Reuse harmonic lifts on an exactly embedded trial/test trace subspace.

injection maps the new trace coefficients to the prepared trial space. test_injection defaults to the same map, or independently restricts the Petrov test space. Both maps must have full column rank. Their geometrical meaning is supplied by the caller; this operation does not project or fit unresolved traces. The local operator, source and retained spaces are fixed.

Multiscale HHO from constrained local energy minimization on triangles.

The cell and face unknowns are integral moments, an invertible change of coordinates from the polynomial coefficients in Chaumont-Frelet et al. (2022). Local reconstructions minimize diffusion energy subject to these moments. The cell moments are statically condensed, leaving the represented face system. The exact diffusion form is symmetric; its assembled operator is retained without a further projection onto its symmetric part. Both the projected-source formulation (4.6) and reconstructed-source variant (5.1) are available. The face-only case m=-1 requires the latter variant.

solve_mshho(mesh: TriangleMesh, *, permeability: Any = 1.0, source: Any = 0.0, dirichlet: Any = 0.0, skeleton: SkeletonSpace | None = None, neumann: dict[int, Any] | None = None, mean_pressure: float = 0.0, cell_degree: int = 0, degree: int = 2, local_refinement: int = 4, source_variant: Literal['projected', 'reconstructed'] = 'projected', quadrature_order: int = 6, solver: str = 'scipy', local_refinement_precision: Literal['double', 'extended'] = 'double') -> MsHHOSolution

Solve multiscale HHO with local continuous Pk energy reconstructions.

cell_degree=m and skeleton specify the coarse moment spaces; degree controls the local numerical approximation only. For m=-1, select source_variant='reconstructed' as required by Remark 5.4. Exterior Dirichlet data are imposed through face moments. neumann maps selected exterior faces to physical outward flux -A grad(u) dot n; remaining faces have Dirichlet data. All-Neumann problems require compatible total source/flux and fix the volume mean through mean_pressure. This finite-dimensional realization uses the same local variational spaces as primal MHM, allowing exact discrete equivalence checks for polynomial/projected sources. Theorem 5.1 concerns ideal local PDE solves; finite-space agreement requires matched spaces and represented source loads. The ideal projected-source flux is H(div) conforming (Lemma 4.2), whereas the reconstructed-source flux need not be (Remark 5.3). Finite local raw gradients require their own conformity check. Theorem 6.3 gives an energy estimate with piecewise H^(m+1) source and H^(ell+1) A grad(p), regular matching macrocell submeshes and material-dependent constants; it supplies no pressure L2 rate by itself. local_refinement_precision='extended' retains reconstruction, represented Galerkin operators and moment-coordinate digits through defect corrections. Factors remain binary64 and the original residual criterion stays 1e-10; platforms without a wider long-double type raise SolverUnavailableError.

Three-dimensional MsHHO on tetrahedral and convex polyhedral macro partitions.

solve_mshho_3d(mesh: TetraMesh | PolyhedralMesh, *, permeability: Any = 1.0, source: Any = 0.0, dirichlet: Any = 0.0, neumann: dict[int, Any] | None = None, skeleton: TriangularSkeleton | PolygonalSkeleton3D | None = None, cell_degree: int = 0, degree: int = 2, local_refinement: int = 2, source_variant: Literal['projected', 'reconstructed'] = 'projected', quadrature_order: int = 5, mean_pressure: float = 0.0, solver: str = 'scipy', local_refinement_precision: Literal['double', 'extended'] = 'double') -> MsHHO3DSolution

Solve MsHHO with tetrahedral Pk energy lifts and polynomial face moments.

Polyhedral faces retain one polynomial on each original polygon. Tetrahedral faces additionally permit independent aligned subdivisions. neumann prescribes physical outward Darcy flux; other exterior faces prescribe pressure. Pure Neumann data require compatibility and fix the volume mean. The face-only case m=-1 requires the reconstructed-source variant. Local numerical spaces must represent all independent volume and face moments. Theorem 5.1 concerns ideal local PDE solves, rather than arbitrary finite tetrahedral lifts. Theorem 6.3 bounds energy with piecewise H^(m+1) source, H^(ell+1) A grad(p), regular matching macrocell submeshes and material-dependent constants. Pressure L2 rates require separate evidence. Projected-source ideal flux conformity (Lemma 4.2) does not apply to reconstructed-source loads (Remark 5.3) or to generic finite local raw gradients without a conformity check. local_refinement_precision='extended' uses the same preserved moment arithmetic as the two-dimensional and polygonal solvers, without changing local spaces, represented operators or the original 1e-10 residual check.

Oriented geometry operations for polygons and star-shaped polyhedral cells.

clip_polygon(polygon: ArrayLike, axis: int, level: float, lower: bool) -> FloatArray

Intersect an ordered convex planar polygon with a closed Cartesian half-plane.

polygon contains finite real vertices with shape (n,2), without a repeated closing vertex. axis is 0 or 1. lower=True retains coordinate >= level; lower=False retains coordinate <= level. Output vertices follow the input's clockwise or counterclockwise orientation. An empty intersection has shape (0,2); point/segment intersections and repeated boundary vertices are retained for the caller's measure checks.

Computation uses binary64 and strict comparisons of represented coordinates. An intersection's clipped coordinate is set exactly to level. No geometric tolerance, vertex merging or grid snapping is inferred; applications requiring those conventions declare them before clipping. Convexity and cyclic ordering are caller hypotheses.

oriented_cell_faces(points: FloatArray, faces: list[IntArray], cells: list[IntArray]) -> tuple[IntArray, ...]

Orient each connected closed cell boundary through opposing directed edge incidences.

The global orientation is selected by its signed surface-volume integral; neither a convex hull nor an assumed interior centroid is used. The caller subsequently checks that the surface is embedded by its positive cone mesh.

kernel_center(vertices: FloatArray, origins: FloatArray, outward: FloatArray) -> tuple[FloatArray, float, bool]

Find a strictly interior ball in the intersection of all inward face half-spaces.

Convex cells retain the arithmetic vertex center. Nonconvex cells maximize the ball radius by a four-variable linear program in cell-scaled coordinates. Canonically ordered inequalities make degenerate optimal-center selection independent of face-list order. The returned radius is verified against every original face, independently of the optimizer's feasibility tolerance. An empty or unresolved kernel is rejected rather than filled by a convex hull.

validate_disjoint_cones(points: FloatArray, cells: IntArray) -> None

Reject interior overlap of tetrahedral cones using the convex separating-axis theorem.

Face normals and pairwise edge cross products are a complete axis set for two tetrahedra. Coordinates are translated/scaled before projection; touching faces, edges and the common kernel center do not constitute overlap.

triangle_in_face(corners: FloatArray, triangles: FloatArray, normal: FloatArray) -> bool

Test a fine boundary triangle against the original face's canonical triangulation.

Conforming local refinements lie within a single canonical triangle even when the original polygon is nonconvex. This preserves individual coplanar faces and does not accept a triangle spanning a reentrant notch.

Source modules

Shared real array aliases and exact discrete parameter validation.

positive_int(value: Any, name: str, minimum: int = 1) -> int

Validate a discrete parameter without silently truncating real numbers.

real_array(value: Any, name: str) -> FloatArray

Copy a finite real array into binary64 without discarding imaginary parts.

dyadic_refinement(value: Any, name: str = 'refinement') -> int

Validate a positive power of two for nested simplex subdivisions.

This is a geometric alignment condition, not a polynomial degree or mathematical stability certificate. The returned integer is unchanged.

Shared solves of represented moment-coordinate equations.

solve_moment_system(matrix: Any, rhs: Any, *, solver: str = 'scipy', refinement_precision: Literal['double', 'extended'] = 'double') -> FloatArray

Solve matrix @ result = rhs without discarding represented wider digits.

matrix is a finite real square operator; rhs is a real vector or a matrix of independent right-hand sides. The named factorization validates solvability and owns its native resources within this call. No symmetry, boundary convention or physical meaning of a coefficient is inferred.

Factors use binary64. refinement_precision='extended' retains the original operator and forcing in wider storage for two additional defect corrections and the unchanged 1e-10 check of every original RHS column. It requires a genuinely wider platform long-double representation and returns wider coefficients; it never substitutes binary64 silently. Local reconstructions, condensed operators and gauged moment equations in all supported spatial dimensions use this same operation.

Dimension-independent physical Neumann extension of the Robin MH method.

RobinBoundaryFace dataclass

Integrated boundary data in one scalar Robin multiplier basis.

The mass matrix pairs multiplier densities and auxiliary pressure traces. Coefficient is sigma.n with the outward exterior normal. Load integrates the prescribed physical flux q.n; constant represents one in this basis.

solve_robin_neumann(system: HybridSystem, trace_size: int, faces: tuple[RobinBoundaryFace, ...], volume_weights: tuple[FloatArray, ...], volume: float, mean_pressure: float, pure_neumann: bool, solver: str, precision: Literal['double', 'extended']) -> tuple[HybridSolution, dict[int, FloatArray], Any, FloatArray]

Solve symmetric boundary-pressure equations, checking physical compatibility.

Rows are S lambda + J rho=t-g_D and J.T lambda+R rho=h_N. Local coercive inverses are unchanged. The optional mean fixes the full reconstructed volume pressure, not an average of trace coefficients.

Petrov-Galerkin MHM with residual face enrichment for two-dimensional Darcy.

Equations (27)--(34) of Fernando, Martins, Pereira and Valentin (2023) use the conormal A grad(p).n. Here all stored multipliers use physical Darcy flux q.n=-A grad(p).n. Only the enriched multiplier is macro conservative.

solve_pgmhm(mesh: TriangleMesh | PolygonMesh, *, stabilization_parameter: float, permeability: Any = 1.0, source: Any = 0.0, dirichlet: Any = 0.0, neumann: dict[int, Any] | None = None, skeleton: SkeletonSpace | None = None, degree: int = 2, local_refinement: int = 1, local_meshes: tuple[TriangleMesh, ...] | None = None, quadrature_order: int = 6, ellipticity_lower_bound: float | None = None, mean_pressure: float = 0.0, solver: str = 'scipy', local_solver: str = 'scipy', local_refinement_precision: Literal['double', 'extended'] = 'double', backend: Literal['serial', 'thread', 'process'] = 'serial', workers: int | None = None) -> PGMHMSolution

Solve equations (27)--(34), including the conservative residual enrichment.

stabilization_parameter=alpha is explicit: the theorem requires alpha sufficiently small but does not provide a universal numerical threshold. The weight is alphaK_min/(2H_E), using the full macroface length, not the lengths of refined skeleton segments. The local degree must satisfy k>=ell+2. The stated uniform theorem assumes ell>=1; ell=0 is included among the article's numerical experiments. Rank/residual checks establish neither a uniform inf-sup bound nor parameter-independent conditioning.

Dirichlet faces contribute jumps against their prescribed trace. Physical Neumann faces are fixed essentially and receive no jump penalty or flux enrichment. Compatible pure Neumann data use the physical mean-pressure gauge; the local enrichment has zero volume mean. The boundary projection in equation (32) is evaluated on the common fine-trace partition so that its Pk moments agree with the Dirichlet functional in equation (29). local_refinement_precision="extended" retains residual-correction digits in both the local response and its residual enrichment. It leaves the operator and residual tolerances unchanged and requires a platform with a wider long-double type, as in :func:pymhm.linalg.linear.solve_linear.