Darcy in two dimensions
Primal, mixed, analytical and residual-enriched Darcy discretizations.
User-defined formulations compose these public numerical owners with
LocalEquations and Equation; their editable providers are listed in the
variational guide.
Physical field records provide optional interpretation of the resulting
coefficients. Named FieldDefinition views require no physical solver object.
Small, auditable triangular finite-element assembly and quadrature kernels.
triangle_quadrature(order: int = 4) -> tuple[FloatArray, FloatArray]
Return barycentric Duffy-product Gauss points and unit-sum weights.
p1_geometry(mesh: TriangleMesh) -> tuple[FloatArray, FloatArray]
Return affine P1 gradients per cell and positive integration areas.
p1_operators(mesh: TriangleMesh, diffusion: Any = 1.0, source: Any = 0.0, reaction: Any = 0.0, advection: Any = (0.0, 0.0), order: int = 4) -> tuple[Any, Any, FloatArray]
Assemble diffusion, mass and load with conventional advective transport.
The bilinear form is (K grad u,grad v)+(beta.grad u,v)+(c u,v).
This operator is not the skew conservative Robin formulation of RAD-MHM.
face_integration(mesh: TriangleMesh, cell: int, fine: TriangleMesh, skeleton: SkeletonSpace) -> tuple[FloatArray, FloatArray]
Integrate skeleton basis against P1 traces and fine boundary flux DOFs.
Quadrature is split at both fine-mesh vertices and skeleton breaks. The
second result maps trace coefficients to integrated outward RT0 boundary
fluxes, ordered by fine.boundary_faces.
boundary_data(skeleton: SkeletonSpace, dirichlet: Any, neumann: dict[int, Any] | None = None, *, neumann_components: dict[int, dict[int, Any]] | None = None, order: int = 5) -> tuple[FloatArray, dict[int, float]]
Integrate primary fields and L2-project full or componentwise boundary fluxes.
neumann_components[face][component] prescribes a Cartesian component
of the multiplier as a scalar field. Other components retain Dirichlet
data. Full and componentwise Neumann declarations cannot share a face.
Components refer to the skeleton's interleaved Cartesian ordering; no
normal/tangential rotation or outward-sign conversion is performed here.
rt0_operators(mesh: TriangleMesh, permeability: Any = 1.0, source: Any = 0.0, order: int = 4) -> tuple[Any, Any, FloatArray]
Assemble RT0 flux mass, integrated divergence and P0 source moments.
A flux DOF is the integral across a globally oriented face. In a cell the
basis associated with side (i,j) is sign*(x-opposite)/(2*area).
rt0_evaluate(mesh: TriangleMesh, flux: FloatArray, barycentric: FloatArray) -> FloatArray
Evaluate RT0 at common (q,3) or cellwise (t,q,3) barycentric points.
Continuous nodal triangle topology and Basix Pk tabulation.
multiindices(degree: int) -> IntArray
Order barycentric lattice nodes by vertices, oriented edges and interior.
reference_values(degree: int, bary: FloatArray) -> FloatArray
Evaluate native equispaced triangular Pk values in PyMHM node order.
bary is a finite real array of shape (points, 3); evaluation is allowed
outside the reference triangle on the unit-sum barycentric hyperplane.
Coordinates and returned values use binary64, with the same node order as
:func:multiindices. Basix tabulates values without allocating derivatives.
reference_basis(degree: int, bary: FloatArray) -> tuple[FloatArray, FloatArray, FloatArray]
Return Basix Pk values and derivatives in the canonical barycentric extension.
Shapes are (points, basis), (points, basis, 3) and
(points, basis, 3, 3). The extension is p(lambda_1,lambda_2):
lambda_0 derivatives are zero, while the other entries are Cartesian
reference derivatives. Contract with all three physical barycentric
gradients to obtain physical derivatives. Input coordinates must sum to
one; off-triangle points on that hyperplane are supported. The node order
remains :func:multiindices; independent off-hyperplane extensions are
not part of the finite-element field.
nodal_space(mesh: TriangleMesh, degree: int) -> tuple[IntArray, FloatArray]
Enumerate continuous nodal unknowns without merging geometrically close nodes.
tabulate(mesh: TriangleMesh, degree: int, bary: FloatArray, *, backend: Literal['portable', 'basix'] = 'basix') -> tuple[IntArray, FloatArray, FloatArray, FloatArray, FloatArray]
Return DOFs, coordinates, values, physical gradients and physical Hessians.
Basix supplies equispaced Pk values and Cartesian derivatives in
:func:multiindices order. Points must sum to one. The affine geometry,
continuous topology and nodal coefficient order are explicit PyMHM data.
backend='portable' is a compatibility spelling for this same Basix
execution; there is no separate polynomial implementation.
element_tabulate(mesh: TriangleMesh, degree: int, bary: FloatArray, *, backend: Literal['portable', 'basix'] = 'basix') -> tuple[IntArray, FloatArray, FloatArray, FloatArray, FloatArray]
Tabulate Pk at separate barycentric quadrature points in each triangle.
bary has shape (cells, points, 3). All returned basis arrays retain
the leading cell axis, so material intersections do not require evaluating
every element's basis at every other element's quadrature points.
backend follows :func:tabulate; tabulation preserves the declared
continuous topology and nodal coefficient order.
trace_coupling(coarse: TriangleMesh, cell: int, fine: TriangleMesh, skeleton: SkeletonSpace, degree: int) -> FloatArray
Integrate signed nodal traces over independently subdivided macrofaces.
scalar_operators(mesh: TriangleMesh, degree: int, *, diffusion: Any = 1.0, source: Any = 0.0, reaction: Any = 0.0, advection: Any = (0.0, 0.0), skew_advection: bool = False, order: int = 6, element_backend: Literal['portable', 'basix'] = 'basix') -> tuple[Any, Any, FloatArray]
Assemble continuous Pk scalar diffusion, reaction and advective operators.
With skew_advection=False the convection term is (beta.grad(u),v).
Its antisymmetric part is selected when skew_advection=True. A variable
conservative velocity also requires half its divergence in reaction.
Basix tabulates local Pk values and derivatives in the same nodal
coefficient order. element_backend='portable' is a compatibility
spelling for this same execution.
Basix tensor Qk tables and scalar operators on Cartesian rectangles.
Rows retain the declared x-fastest nodal ordering. Product quadratures have unit-sum weights; physical cell measures are applied by the operators.
quadrilateral_quadrature(order: int = 4) -> tuple[FloatArray, FloatArray]
Return product Gauss points on [0,1]^2 and unit-sum integration weights.
cardinal_polynomials(degree: int) -> tuple[FloatArray, ...]
cached
Export the reordered Basix interval basis in ascending power coordinates.
This representation conversion supports coefficient archives; evaluation and derivatives use Basix directly. The orthonormal native polynomial j equals sqrt(2j+1) times the shifted conventional Legendre polynomial.
qk_basis(degree: int, points: Any) -> tuple[FloatArray, FloatArray]
Evaluate equidistant Qk basis and reference gradients at unit-square points.
Nodes and basis functions are ordered first in x, then in y. Basix interval factors supply values and derivatives, including at nodal points; their tensor product retains the historical coefficient order.
qk_space(mesh: CartesianMacroMesh, degree: int) -> tuple[IntArray, FloatArray]
Return the shared nodal DOFs and coordinates of a continuous Qk space.
grid_resolves_material(mesh: CartesianMacroMesh, field: CartesianCellField) -> bool
Determine whether every material line coincides with a fine-grid line.
rectangle_intersection_quadrature(origins: FloatArray, spacing: FloatArray, field: CartesianCellField, order: int) -> tuple[FloatArray, FloatArray]
Integrate each rectangle separately on its exact intersections with pixels.
Reference points have shape (cells, points, 2); weights have shape (cells, points) and sum to one on each rectangle. Zero-weight padding combines variable intersection counts without extrapolating the material.
quadrilateral_operators(mesh: CartesianMacroMesh, degree: int = 1, *, permeability: Any = 1.0, source: Any = 0.0, order: int = 4) -> tuple[Any, Any, FloatArray]
Assemble sparse Qk diffusion, mass and load using bounded cell batches.
CartesianCellField coefficients are integrated over exact intersections with material pixels, retaining a uniform-rule fast path on aligned grids. Arbitrary callback discontinuities are not inferred: those require an aligned grid or a declared quadrature-sensitivity study.
quadrilateral_trace_coupling(coarse: CartesianMacroMesh, cell: int, fine: CartesianMacroMesh, skeleton: SkeletonSpace, degree: int) -> FloatArray
Integrate signed Qk traces, splitting at fine and skeleton breakpoints.
Exact scalar energy and face moments in a metric-quadratic simplex basis.
The scalar basis has constant, mean-zero affine and mean-zero metric-quadratic coordinates. These finite-dimensional operators can be composed with arbitrary trial/test trace equations; no global method, gauge or solver is selected.
AnalyticDarcySpace
dataclass
Constant, mean-zero affine and mean-zero metric-quadratic local basis.
evaluate(points: FloatArray) -> tuple[FloatArray, FloatArray]
Evaluate basis and gradients at physical points, preserving metric units.
integrate_rt0(pressure_mean: float, centroid_flux: Any, divergence: float) -> FloatArray
Integrate an RT0 field into the unique potential with its prescribed mean.
The field is q(c) + divergence/d * (x-c). The result contains the
coefficients in this space, with -K grad(p)=q identically. A supplied
classical RT0 pressure mean is preserved; this operation does not add
the variable-source moment correction of the analytical MHM equations.
metric_quadratic_operators(mesh: TriangleMesh | TetraMesh, cell: int, permeability: Any = 1.0, source: Any = 0.0, *, quadrature_order: int = 8) -> tuple[FloatArray, FloatArray, FloatArray, FloatArray, AnalyticDarcySpace]
Build exact harmonic lift equations and all variable-source moments.
Permeability is one constant SPD tensor on this macrocell. The local trial
span is {1, x_i-center_i, (x-center).K^-1.(x-center)/2}. Its gradients
generate exactly the RT0 vector space after multiplication by K. Basis
functions other than the constant have zero physical mean. The source
response in this small polynomial space is used for the global load only;
a full variable-source potential requires a separate local Neumann solve.
Cartesian scalar operators from explicit one-dimensional weighted factors.
interval_nodal_quadrature(count: int, interval: tuple[float, float], degree: int, order: int) -> tuple[FloatArray, FloatArray, FloatArray, FloatArray, np.ndarray, float]
Return points, unit-interval weights, values, derivatives, DOFs and spacing.
Physical points have shape (cells, quadrature); value/derivative tables have
shape (quadrature, degree+1). Derivatives are in reference coordinates and
DOFs share equispaced interval nodes. interval_weighted_operators applies
the physical spacing to mass, stiffness and load. No boundary is prescribed.
interval_weighted_operators(data: tuple, values: FloatArray) -> tuple
Assemble physical weighted mass, stiffness and load in shared nodal order.
data is the explicit interval_nodal_quadrature tuple, and values has its
(cells, quadrature) shape. Signed weights are admissible; PDE coefficient
positivity is a separate condition. Return two CSR matrices and one vector.
require_positive_separated(values: list[tuple[FloatArray, ...]]) -> None
Require sum(a(x)*b(y)) positive at its declared tensor sampling points.
Each pair contains finite sampled x/y factor arrays. Signed terms are admissible. The check uses a bound when sufficient and otherwise bounded temporary blocks. It does not certify the coefficient between these points.
separable_diffusion_operators(mesh: CartesianMacroMesh, degree: int, *, permeability: SeparableField = _UNIT, source: SeparableField = _ZERO, order: int = 6) -> tuple[sparse.csr_matrix, sparse.csr_matrix, FloatArray]
Assemble scalar diffusion, mass and load by exact Kronecker identities.
The quadrature rule uses at least degree+1 Gauss points per coordinate.
Coefficients and source are explicit separated sums; their spatial variation
is integrated with this same rule, without homogenization or low-rank fitting.
Unknowns retain the native ordering ix + (nx*degree+1)*iy.
Finite separated scalar fields without rank fitting or coefficient averaging.
SeparableField
dataclass
Represent sum(a(x)*b(y) for a,b in terms) without approximation.
Each factor is a finite real scalar or a callback accepting a one-dimensional coordinate array. An empty sum represents zero. Signed terms are allowed; a diffusion assembler separately verifies positivity at its Gauss points.
__post_init__() -> None
Normalize term pairs and reject nonscalar, nonfinite constant factors.
__call__(points: FloatArray) -> FloatArray
Evaluate the declared scalar field at physical XY point pairs.
factor_values(value: Any, points: FloatArray) -> FloatArray
Evaluate a real one-coordinate factor with a finite scalar/vector contract.
Physical boundary functionals and nodal data on Cartesian conforming fields.
quadrilateral_boundary_data(mesh: CartesianMacroMesh, degree: int, *, dirichlet: Any = 0.0, neumann: dict[int, Any] | None = None, order: int = 6) -> tuple[FloatArray, dict[int, float]]
Return negative outward-flux load and strongly prescribed nodal values.
Neumann keys identify exterior faces and their physical outflow density. Remaining exterior faces prescribe dirichlet in the executed Qk node order. Intersections share the same physical callback value. The returned load is added to the volume functional; no source, operator or solver is chosen.
Executed metric-quadratic harmonic fields and independent nodal source fields.
AnalyticDarcySolution
dataclass
Analytical trace-driven pressure/RT0 flux plus a separate source potential.
__post_init__() -> None
Archive each executed nodal source basis alongside its coefficient vector.
pressure_update(cell: int, points: FloatArray) -> tuple[FloatArray, FloatArray]
Return p0+p_lambda and its exact RT0 physical flux, without p_f.
evaluate(cell: int, bary: FloatArray) -> tuple[FloatArray, FloatArray]
Evaluate the full pressure and flux on every source-mesh cell at barycentric points.
errors(pressure: Any, flux: Any, order: int = 8) -> tuple[float, float]
Integrate pressure and physical flux errors on the actual source meshes.
Physical evaluation records for independently assembled conforming baselines.
ConformingQuadrilateralSolution
dataclass
One global continuous Qk field, with physical evaluation at arbitrary points.
evaluate(points: FloatArray) -> tuple[FloatArray, FloatArray]
Return pressure and gradient without interpolation or material averaging.
physical_flux(points: FloatArray) -> FloatArray
Evaluate -K grad(p) using the coefficient's declared pointwise convention.
Shared assembly operations for physical normal-flux Darcy formulations.
Element families own tabulation, quadrature, orientation and constant-mode coordinates. These functions assemble supplied bases without changing them.
triangle_pressure_basis(degree: int, bary: FloatArray) -> FloatArray
Tabulate complete cardinal DG pressure, with constant P0 equal to one.
degree is an already validated nonnegative mathematical degree.
Barycentric points have shape (..., 3); the trailing output axis is
the nodal pressure coefficient order used by the simplex FEM tabulator.
triangle_mixed_operators(mesh: TriangleMesh, bary: FloatArray, weights: FloatArray, material: Any, source: Any, flux: FloatArray, divergence: FloatArray, pressure: FloatArray, dofs: IntArray, flux_size: int) -> tuple[Any, Any, FloatArray]
Assemble mass, positive divergence moments and pressure-tested source.
Supplied flux bases already use global-oriented local DOF coordinates.
Arrays begin with (fine_cell, quadrature_point); flux has trailing
(local_flux_dof, 2) and pressure/divergence have a basis axis.
weights are normalized triangle weights, so physical integrals also
multiply by mesh.areas. Material values are scalar or SPD tensors.
DG pressure coefficients are concatenated by fine cell. The returned
divergence has positive sign; the Darcy saddle applies its minus sign.
triangle_pressure_integrals(mesh: TriangleMesh, bary_weights: FloatArray, pressure: FloatArray) -> FloatArray
Integrate each DG pressure basis function in concatenated fine-cell order.
normal_flux_blocks(mass: Any, divergence: Any, source: FloatArray, boundary_dofs: IntArray, trace_mapping: FloatArray) -> tuple[Any, FloatArray, FloatArray]
Build the mixed Neumann saddle, trace coupling and signed source load.
Local coefficients are (flux, pressure, boundary_pressure). The third
equation prescribes exterior flux moments by trace_mapping @ trace;
the second equation is divergence @ flux = source. The first equation
is the weak law K^-1 flux + grad(pressure) = 0 with auxiliary boundary
pressure. Exterior DOFs and mapping rows must have the same order.
Constant-mode kernels, physical pressure means and any scaling remain
caller-owned because nodal and modal pressure coordinates differ.
Physical oriented RT normal-moment pairings on planar macrofaces.
rt_trace_map(coarse: TriangleMesh, cell: int, fine: TriangleMesh, skeleton: SkeletonSpace, degree: int) -> FloatArray
Integrate canonical RT normal moments of the oriented macroface flux density.
Tensor-product Raviart--Thomas elements with independent interior enrichment.
On each rectangle, normal traces have degree k and interior vector functions have RT order s=k+n. The pressure is Q_s. Face DOFs are integral Legendre moments; cell DOFs are coefficients of zero-normal-trace polynomial bubbles. The construction implements the face/interior separation used by the mixed MHM family of Duran et al. (2019), without copying a reference implementation.
tensor_rt_moment_coefficients(degree: int, enriched_degree: int) -> FloatArray
cached
Express declared face lifts and bubble coordinates in native rectangular RT.
tensor_rt_basis(mesh: CartesianMacroMesh, degree: int, enrichment: int, points: FloatArray) -> tuple[FloatArray, FloatArray, FloatArray]
Evaluate globally oriented RT vectors, divergence and modal Q_s pressure.
points are reference coordinates in [0,1]^2, either shared (q,2) or
cell dependent (cells,q,2). Vector/divergence arrays
include a leading cell axis. The pressure basis is cell independent and
ordered with x degree varying fastest. The Piola transform preserves all
normal moments even on anisotropic rectangles.
tensor_rt_dofs(mesh: CartesianMacroMesh, degree: int, enrichment: int) -> IntArray
Map conforming face moments and independent cell bubbles to global DOFs.
tensor_rt_trace_map(mesh: CartesianMacroMesh, cell: int, fine: CartesianMacroMesh, skeleton: SkeletonSpace, degree: int) -> FloatArray
Integrate each oriented macro flux against fine-edge Legendre moments.
tensor_rt_operators(mesh: CartesianMacroMesh, degree: int, enrichment: int, permeability: Any, source: Any, order: int) -> tuple[Any, Any, FloatArray]
Assemble mixed inverse-permeability mass and divergence moments.
Broken pressure and conservative skeleton flux of a Darcy solve.
Primal pressure coefficients are continuous Pk nodal values on each macrocell.
Mixed pressure coefficients
are P0 cell values and the flux is RT0 conforming, with integrated face DOFs.
For the primal formulation, flux stores samples of -K grad(p) at
fine-cell centroids; these are not piecewise-constant coefficients when
the degree exceeds one or K varies. Error norms evaluate the full field
at quadrature points. Raw primal flux is generally not H(div)-conforming.
l2_error(exact: Any, order: int = 5) -> float
Integrate pressure error with independent higher-order quadrature.
flux_l2_error(exact: Any, order: int = 5) -> float
Integrate physical flux error; primal flux is evaluated as -K grad p.
conservation_residuals(order: int | None = None) -> FloatArray
Return macro flux minus source, using assembly quadrature by default.
Supply a higher order to measure source integration error separately.
fine_conservation_residuals() -> tuple[FloatArray, ...]
Return RT0 cellwise mass defects; reject raw primal flux diagnostics.
Locally H(div)-conforming RTm flux and discontinuous nodal Pm pressure.
Mathematical degree m>=0 corresponds to Basix RT element degree m+1. MHM traces prescribe physical normal flux, while classical solutions use a globally conforming fine mesh. The latter have no skeletal multiplier.
l2_error(exact: Any, order: int = 6) -> float
Integrate the discontinuous pressure error without centroidal substitution.
flux_l2_error(exact: Any, order: int = 6) -> float
Integrate the complete physical H(div) vector error.
divergence_l2_error(exact: Any, order: int = 6) -> float
Integrate the full divergence error, distinguishing projection from pointwise balance.
fine_equilibrium_residuals() -> tuple[FloatArray, ...]
Return every DG-Pm moment of div(q)-f on each fine triangle.
fine_conservation_residuals() -> tuple[FloatArray, ...]
Return integrated fine-cell flux/source imbalance via the partition of unity.
conservation_residuals() -> FloatArray
Sum fine-cell physical balances over each local mesh/macrocell.
normal_flux_residuals() -> tuple[FloatArray, ...]
Compare all local exterior RT normal moments with the physical MHM trace.
Conforming BDM(k,n) flux and complete discontinuous pressure from a MHM solve.
Flux stores integrated normal moments and cell moments as specified in
pymhm.fem.hdiv.bdm_family. The default BDM2/P1 uses the equivalent coordinates
in pymhm.fem.hdiv.bdm. Pressure has degree k+n-1 and cardinal coefficients.
The skeletal multiplier is the physical flux q.n in the macro orientation.
l2_error(exact: Any, order: int = 5) -> float
Integrate broken pressure error against an analytical pressure.
flux_l2_error(exact: Any, order: int = 5) -> float
Integrate physical BDM flux error without gradient recovery or smoothing.
divergence_l2_error(exact: Any, order: int = 5) -> float
Integrate divergence error; exact=source measures strong equilibrium error.
fine_equilibrium_residuals() -> tuple[FloatArray, ...]
Return all complete pressure moments of div(q)-f in every fine triangle.
fine_conservation_residuals() -> tuple[FloatArray, ...]
Return integrated div(q)-f over each fine cell, using assembly quadrature.
conservation_residuals() -> FloatArray
Return macro flux minus source integrals, with the signed skeleton flux.
normal_flux_residuals() -> tuple[FloatArray, ...]
Return all boundary normal-moment differences between BDM flux and skeleton.
Repeated physical RT evaluation
RTField(mesh, coefficients, degree) stores an immutable copy of canonical
Raviart–Thomas coefficients and validates it once. evaluate(points, cell_indices)
uses the prescribed incident cell at each physical point; it does not average
traces across interfaces. The flux and divergence use the same Piola map and
moment conventions as rt_evaluate_points. Point and cell-membership checks
remain active for every batch. Serialization preserves immutable coefficient
storage and its numerical precision. Its complete API is documented with the
RT element family.