Skip to content

Robin MH and three-field MH²M

Native two- and three-dimensional local maps, physical boundary conventions and independent trace spaces.

The Robin multiplier is \((q-p\sigma)\cdot n\), while MH²M's local conormal is \(A\nabla p\cdot n=-q\cdot n\) and its global trace represents pressure. MH²M uses a boundary-mean-zero local pressure complement. The reconstruction and gauge conventions are documented in the source API below; the original formulations are given by Barrenechea, Gomes and Paredes (2024) and de Barros, Madureira and Valentin (2026).

Barrenechea et al. (2024) prove their Robin finite-element results in two dimensions and indicate the three-dimensional extension. De Barros et al. (2026, version 3) define their functional maps for dimensions two and three, but their sufficient polynomial-space Fortin constructions are two-dimensional. The implemented tetrahedral paths have separate analytical/native verification; passing an algebraic rank check does not transfer those mesh-uniform stability theorems.

All API families

Three-field multiscale Darcy method of de Barros, Madureira and Valentin.

The pressure trace Gamma is globally continuous, while Lambda has independent outward conormal coefficients on each macrocell. Equations (25)--(29) of the 2026 preprint define two local Neumann inverses and a global pressure-trace system. The paper's lambda=A grad(p).n is minus the physical Darcy flux q.n.

solve_mh2m(mesh: TriangleMesh | PolygonMesh, *, permeability: Any = 1.0, source: Any = 0.0, dirichlet: Any = 0.0, neumann: dict[int, Any] | None = None, pressure_trace: PressureTraceSpace | None = None, flux_space: SkeletonSpace | None = None, degree: int = 1, local_refinement: int = 4, local_meshes: tuple[TriangleMesh, ...] | None = None, quadrature_order: int = 6, mean_pressure: float = 0.0, solver: str = 'scipy', local_solver: str = 'scipy') -> MH2MSolution

Solve the 2D three-field MH2M, equations (25)--(29), with independent spaces.

Gamma is continuous across macro vertices; Lambda is discontinuous and is duplicated on each incident macrocell. Default Gamma is P1, Lambda is P0. Macroelements may be triangles or simple polygons, including nonconvex polygons triangulated by the shared local-mesh builder. Each polygon side retains its independent conormal space and continuous pressure endpoints. Pressure degree and local refinement specify Vh independently. Insufficient Lambda/Vh injectivity or a rank-deficient global trace system is rejected. This numerical rank check is not a mesh-uniform inf-sup proof. Sufficient mesh conditions (M1),(M2) and the exceptions of Remark 16 remain relevant. The optional local_meshes supplies one validated conforming triangular partition per macrocell. It replaces the uniform local_refinement geometry without changing Gamma, Lambda, or either Neumann map.

neumann maps exterior faces to outward physical q.n. Other exterior faces interpolate dirichlet at Gamma nodes. Pure Neumann data require compatibility; mean_pressure then fixes the volume average of the full reconstructed pressure. Nonhomogeneous and mixed data are the consistent extension of the paper's homogeneous Dirichlet examples.

Multiscale Hybrid diffusion with coercive Robin local problems (2024).

The multiplier is lambda=(-K grad(p)-p sigma).n, where sigma=nu*(x-origin)/2. It is not the physical Darcy normal flux. The local operators and the Dirichlet condensed operator are symmetric positive definite under the stated ellipticity and trace-compatibility conditions. Physical Neumann conditions use a symmetric boundary-pressure extension; its global system is indefinite.

solve_mh(mesh: TriangleMesh | PolygonMesh, *, permeability: Any = 1.0, source: Any = 0.0, dirichlet: Any = 0.0, neumann: dict[int, Any] | None = None, mean_pressure: float = 0.0, skeleton: SkeletonSpace | None = None, degree: int = 2, local_refinement: int = 4, quadrature_order: int = 6, robin_parameter: float | None = None, origin: Any = None, ellipticity_lower_bound: float | None = None, solver: str = 'scipy', local_solver: str = 'scipy', local_refinement_precision: Literal['double', 'extended'] = 'double', refinement_precision: Literal['double', 'extended'] = 'double', backend: Literal['serial', 'thread', 'process'] = 'serial', workers: int | None = None) -> MHSolution

Solve diffusion with the 2024 Multiscale Hybrid Robin local formulation.

Define sigma=nu(x-origin)/2, with positive robin_parameter=nu. The default origin is the lower corner of the bounding box. Let C_sigma=max_vertex |x-origin|/2. We require nu<=K_min/(4C_sigma**2), the sufficient coercivity estimate using the actual sigma norm. The default is half this bound. This sharpens the paper's conservative substitution of the bounding-domain diameter; it does not change its bilinear form. Literal scalar/tensor and Cartesian materials supply K_min automatically; other callbacks require a certified ellipticity_lower_bound.

Local Pk spaces and signed skeletal spaces are independent. A fine mesh resolving each skeletal segment with two boundary edges gives the mesh condition used by the article, with the corresponding compatible degrees. An algebraically singular trace choice is rejected by the shared solver. Neumann keys identify exterior faces and prescribe the outward physical flux q.n, not the Robin multiplier. Remaining exterior faces prescribe pressure weakly. Auxiliary boundary pressure moments preserve the Robin local operators; this mixed-boundary extension has a symmetric indefinite global system. Pure Neumann problems require integral compatibility and mean_pressure fixes the volume average of the reconstructed field. At nu=0 the elliptic condensation is undefined; use solve_darcy for MHM. Assembly and condensation occur together inside the selected workers. The explicit refinement-precision options reuse the shared residual-checked iterative refinement. Extended accumulation requires a wider long-double type; neither option changes the local or global residual tolerance.

Three-dimensional Robin MH with physical mixed boundary conditions.

solve_mh_3d(mesh: TetraMesh, *, permeability: Any = 1.0, source: Any = 0.0, dirichlet: Any = 0.0, neumann: dict[int, Any] | None = None, mean_pressure: float = 0.0, skeleton: TriangularSkeleton | None = None, degree: int = 2, local_refinement: int = 2, quadrature_order: int = 5, robin_parameter: float | None = None, origin: Any = None, ellipticity_lower_bound: float | None = None, solver: str = 'scipy', local_solver: str = 'scipy', local_refinement_precision: Literal['double', 'extended'] = 'double', refinement_precision: Literal['double', 'extended'] = 'double', backend: Literal['serial', 'thread', 'process'] = 'serial', workers: int | None = None) -> MH3DSolution

Solve tetrahedral Robin MH using sigma=nu(x-origin)/3 and physical q.n data.

This dimensional extension retains the 2024 bilinear form and uses div(sigma)=nu. With C=max_vertex|x-origin|/3, nu<=K_min/(4C*2) gives the same sufficient Young-inequality coercivity bound in three dimensions. The default is half that bound. Material callbacks require a certified ellipticity lower bound. Local Pk and triangular skeletal modes are independently chosen; the local refinement must align their face partitions. Rank checks do not establish a uniform inf-sup constant for arbitrary degree/subdivision pairs.

Neumann values are physical outward fluxes, never Robin multipliers. Pure Neumann data require integral compatibility and a prescribed volume mean. Other exterior faces prescribe pressure weakly. Assembly uses tetrahedral Gaussian integration; discontinuous coefficients require adequate local geometric resolution rather than an implicit enrichment of the trial space.

Three-field MH²M on tetrahedra with independently resolved interface spaces.

solve_mh2m_3d(mesh: TetraMesh, *, permeability: Any = 1.0, source: Any = 0.0, dirichlet: Any = 0.0, neumann: dict[int, Any] | None = None, pressure_trace: PressureTraceSpace3D | None = None, flux_space: TriangularSkeleton | None = None, degree: int = 2, local_refinement: int = 2, quadrature_order: int = 5, mean_pressure: float = 0.0, solver: str = 'scipy', local_solver: str = 'scipy') -> MH2M3DSolution

Solve the three-dimensional MH²M using equations (25)--(29) of version 3.

Gamma is continuous nodal Pk on uniformly refined triangular macrofaces. Lambda is independently subdivided and uses broken Bernstein polynomials; each incident macrocell has its own outward conormal copy. All partitions must align the fine tetrahedral boundary. A singular local Neumann map or pressure-trace system is rejected; this check is not a uniform inf-sup proof. The article's functional local maps include dimensions two and three, but its sufficient polynomial Fortin constructions in Section 6 use two-dimensional edges. Those constructions do not certify arbitrary triangular-face degree and subdivision pairs here.

Neumann data prescribe outward physical q.n. Remaining exterior faces interpolate pressure at Gamma nodes. Pure Neumann data require integral compatibility and mean_pressure fixes the full volume-mean reconstruction. Local operators use positive tetrahedral Gaussian integration. Discontinuous material requires geometrically adequate local resolution.

Continuous pressure traces and physical boundary quadrature on tetrahedra.

PressureTraceSpace3D dataclass

Continuous nodal Pk functions on a uniformly triangulated tetrahedral skeleton.

Each original triangular face has subdivisions**2 subtriangles. Polynomial degrees and subdivision counts are independent of the broken conormal space. Nodes on shared macroedges and vertices are identified by exact integer barycentric weights, never by rounded physical coordinates. The uniform face degree and subdivisions ensure conforming edge traces.

size: int property

Return the global continuous pressure-trace dimension before boundary elimination.

__post_init__() -> None

Construct topological node identities on each canonical subtriangle.

cell_dofs(cell: int) -> IntArray

Return the sorted pressure-trace nodes incident to one macrocell.

evaluate(face: int, bary: FloatArray) -> FloatArray

Evaluate all face-local nodal functions at original-face barycentric points.

boundary_rules(coarse: TetraMesh, cell: int, fine: TetraMesh, degree: int, order: int) -> tuple[tuple[int, IntArray, FloatArray, FloatArray, FloatArray, FloatArray], ...]

Tabulate physical fine-boundary facets and their original macroface coordinates.

Each tuple contains macroface index, incident volume DOFs, nodal basis, physical XYZ points, physical weights, and macroface barycentric coordinates. Fine tetrahedra must conform to the macrocell boundary.

broken_face_basis(skeleton: TriangularSkeleton, face: int, bary: Any) -> FloatArray

Evaluate DG or macroface-C0 conormal modes in the skeleton's coefficient order.