Local finite-element operators
Reusable local volume and trace operators, quadrature and reference-space operations. The documentation below is generated from their source docstrings. User-written variational problems combine these operations with their own coefficients and coupling equations.
Source modules
BDM-family forms in the declared normal and interior coefficient convention.
bdm_pressure_basis(family: BDMFamily, bary: FloatArray) -> FloatArray
Return a cardinal complete pressure space with partition of unity.
bdm_basis(family: BDMFamily, mesh: TriangleMesh, bary: FloatArray) -> tuple[FloatArray, FloatArray]
Tabulate BDM fluxes in their declared normal/interior moment coordinates.
bdm_evaluate(family: BDMFamily, mesh: TriangleMesh, flux: FloatArray, bary: FloatArray) -> tuple[FloatArray, FloatArray]
Evaluate vectors in the same executed coordinate convention as assembly.
bdm_trace_map(family: BDMFamily, mesh: TriangleMesh, cell: int, fine: TriangleMesh, skeleton: SkeletonSpace) -> FloatArray
Map the physical multiplier into the family's oriented boundary moments.
Mapped H(div) volume integration and declared normal-density trace moments.
HexSkeleton
dataclass
Tensor-Qk reference flux-density traces with aligned uniform face subdivisions.
face_size: int
property
Return the number of reference flux-density modes on each macroface.
size: int
property
Return the global number of scalar normal-flux trace modes.
__post_init__() -> None
Validate scalar polynomial degree and positive subdivision count.
cell_dofs(cell: int) -> IntArray
Return the six oriented macroface blocks in local side order.
evaluate(points: FloatArray) -> FloatArray
Evaluate reference flux densities, normalized by segment moments.
mapped_rt_operators(mesh: HexMesh, degree: int, permeability: Any, source: Any, order: QuadratureOrder) -> tuple
Assemble mixed operators with bounded quadrature batches and physical material values.
quadrature_slices(mesh: HexMesh, degree: int, count: int) -> Iterator[slice]
Partition one unchanged rule to bound simultaneous physical-basis storage.
mapped_rt_trace_mapping(mesh: HexMesh, cell: int, fine: HexMesh, reference: FloatArray, skeleton: HexSkeleton, degree: int) -> FloatArray
Integrate fine normal moments of each oriented mapped macroface trace.
mapped_boundary_data(skeleton: HexSkeleton, dirichlet: Any, neumann: dict[int, Any], order: int) -> tuple[FloatArray, dict[int, float]]
Integrate weak pressure data and project outward physical normal fluxes.
Basix reference bases in declared H(div) candidate coordinates.
Only polynomial tabulation is delegated here. Physical moment definitions, normal restrictions, bubble coordinates and Piola maps belong to their callers. Native handles are cached in each process and never passed to spawn workers.
vector_tabulation(family: str, cell: str, degree: int, points: FloatArray) -> tuple[FloatArray, FloatArray]
Tabulate a native RT/BDM basis and its Cartesian reference divergence.
degree follows Basix: mathematical RT order m uses degree m+1; BDM
uses its vector polynomial degree. Values have axes (point, DOF, component)
and divergence (point, DOF). Moment numbering is supplied by the caller.
RT inverse-material mass, DG divergence and scalar load integration.
pressure_basis(degree: int, bary: FloatArray) -> FloatArray
Tabulate discontinuous nodal Pm pressure, using a constant for RT0.
rt_operators(mesh: TriangleMesh, degree: int, permeability: Any = 1.0, source: Any = 0.0, order: int = 5) -> tuple[Any, Any, FloatArray, FloatArray]
Assemble inverse-diffusion mass, DG-Pm divergence/load and pressure integral weights.
Explicit degree-dependent Gaussian point-count policies for nodal forms.
nodal_quadrature_order(degree: int, requested: int, *, increment: int = 2) -> int
Return max(requested,degree+increment) as a declared Gaussian point count.
This policy is a count in each tensor/Duffy coordinate, not a polynomial exactness degree or a certificate for nonpolynomial coefficients. The caller declares the nonnegative increment and independently verifies quadrature adequacy for its form, material, source and error norm.
Geometric compatibility of independently declared simplicial integration pieces.
simplex_partition_refines(fine: Any, coarse: Any, *, tolerance: float = 1e-12) -> bool
Check that each fine simplex lies in one coarse simplex of a partition.
Arrays have axes (piece,vertex,ambient_coordinate); codimension is allowed. Vertices suffice because the coarse pieces are convex simplices. This tests the common-integration hypothesis, rather than equality of partition counts. The nonnegative tolerance applies to barycentric coordinates and geometric residual divided by the coarse piece diameter. Partition coverage and nonoverlap remain the declaring geometry owner's separate responsibility.
Original Pk/Qk Helmholtz forms, material cuts and exact realification.
complex_values(value: Any, points: FloatArray) -> Any
Evaluate scalar complex data, rejecting unintended broadcasting and NaNs.
positive_values(value: Any, points: FloatArray) -> FloatArray
Evaluate scalar acoustic material, accepting isotropic planar tensor data.
stretch_values(value: Any, points: FloatArray) -> Any
Evaluate outgoing PML stretches with positive real and nonnegative imaginary parts.
real_vector(value: Any) -> FloatArray
Interleave the real and imaginary coordinates of a complex vector.
complex_vector(value: Any) -> Any
Decode the component-interleaved real coordinates without a sign change.
real_matrix(matrix: Any) -> Any
Represent complex multiplication by [[Re,-Im],[Im,Re]] per scalar entry.
acoustic_space(mesh: Any, degree: int) -> tuple[IntArray, FloatArray]
Return continuous local nodal coordinates and element maps for Pk or Qk.
acoustic_quadrature(mesh: Any, degree: int, material: Any, order: int) -> tuple[IntArray, FloatArray, FloatArray, FloatArray, FloatArray, Any]
Tabulate physical local fields over all intersections with material pixels.
Return element DOFs, physical points, physical weights, basis values, gradients and material samples. Separate terms can use separate material partitions, so density and modulus need not share a pixel grid.
volume_forms(mesh: Any, degree: int, density: Any, modulus: Any, source: Any, order: int, pml_stretch: Any = None) -> tuple[Any, Any, Any]
Assemble diffusion with rho^-1, mass with kappa^-1 and complex forcing.
acoustic_edge_pieces(macro: Any, fine: Any, cell: int, face: int, degree: int) -> Iterator[tuple[FloatArray, IntArray]]
Yield parameter endpoints and local edge DOFs without crossing a fine edge.
edge_rule(macro: Any, fine: Any, cell: int, face: int, degree: int, space: Any, order: int, materials: tuple[Any, ...]) -> Iterator[tuple[IntArray, FloatArray, FloatArray, FloatArray, FloatArray]]
Integrate on the common fine, skeletal and explicit material partition.
material_trace(material: Any, points: FloatArray, interior: FloatArray) -> FloatArray
Evaluate acoustic material on the incident side of a boundary face.
Diffusion Neumann maps in an explicitly declared boundary-mean complement.
The supplied local stiffness has a constant kernel. Boundary coupling represents outward conormal moments, while volume moments are kept for physical gauges. The returned data records the executed zero-mean basis instead of recomputing it.
NeumannMaps
dataclass
Local operators and Eq. (29) lifts in declared nodal/Legendre coordinates.
conormal_lift and conormal_source use the paper's lambda convention.
neumann_energy acts on zero-boundary-average Lambda coordinates.
Boundary averages, rather than volume averages, define the local complement.
neumann_maps(fine: TriangleMesh | TetraMesh, trace_ids: IntArray, stiffness: Any, mass: Any, load: FloatArray, coupling: FloatArray, pairing: FloatArray, constant: FloatArray, solver: str) -> NeumannMaps
Construct the dimension-independent Neumann inverses and Eq. (29) lifts.
Coupling has the unsigned outward-conormal convention. Boundary integrals supply its physical mean; mass supplies the volume mean used only by gauges.
Affine Robin boundary forms and physical nodal volume moments in two dimensions.
These integrate supplied finite-element bases; they do not select a hybrid method, solve a PDE, or infer a multiplier convention. Robin sigma is explicitly parameter*(x-origin)/2 and normals point outward from the supplied fine mesh.
nodal_volume_moments(fine: TriangleMesh, degree: int) -> FloatArray
Integrate every local nodal basis function for the physical pressure mean.
robin_boundary_operator(fine: TriangleMesh, degree: int, parameter: float, origin: FloatArray, order: int) -> Any
Integrate sigma.nuv using the general boundary owner, sigma=parameter*(x-origin)/2.
Affine Robin coefficient adapter for the general boundary bilinear owner.
The user selects sigma=parameter*(x-origin)/3. Every fine boundary integral uses the outward macrocell normal and the declared nodal coefficient basis.
tetra_robin_boundary_operator(mesh: TetraMesh, cell: int, fine: TetraMesh, degree: int, order: int, nu: float, origin: FloatArray) -> Any
Integrate sigma.nuv by supplied fine facet rules, sigma=nu*(x-origin)/3.
Topological tetrahedral Pk nodes and inherited Basix derivatives.
tetra_indices(degree: int) -> IntArray
cached
Return equispaced barycentric nodal multi-indices for any positive degree.
Vertices precede edge interiors, face interiors and cell interiors. Node identity depends on vertex indices and integer barycentric weights, never on rounded physical coordinates.
continuous_tetra_nodes(mesh: Any, degree: int) -> tuple[IntArray, FloatArray]
Enumerate continuous Pk DOFs through exact shared vertex/edge/face identities.
tetra_polynomials(degree: int, bary: FloatArray) -> tuple[FloatArray, FloatArray, FloatArray]
Return native Pk values and canonical barycentric gradients/Hessians.
Basix tabulates equispaced Pk in :func:tetra_indices order. Derivative
axes have length four, with lambda_0 entries zero and Cartesian reference
derivatives in lambda_1 through lambda_3. This declares the extension
p(lambda_1,lambda_2,lambda_3); contraction with physical barycentric
gradients gives the affine physical derivatives. Finite real coordinates
must sum to one, including extrapolation points. Positive degree support
is inherited from Basix; equispaced interpolation retains its conditioning
limitations at high degree.
tetra_values_gradients(degree: int, bary: FloatArray) -> tuple[FloatArray, FloatArray]
Return native Pk values/canonical gradients without allocating Hessians.
The node order, unit-sum point contract and padded four-coordinate
derivative convention follow :func:tetra_polynomials.
Conservative scalar volume forms and explicit strong-residual stabilization.
ScalarTransportOperators
dataclass
Integrated A/load and physical nodal moments, with exact sampled constant-mode flags.
triangle_transport_operators(fine: TriangleMesh, *, degree: int, diffusion: Any = 1.0, velocity: Any = (0.0, 0.0), velocity_divergence: Any = 0.0, diffusion_divergence: Any = (0.0, 0.0), reaction: Any = 0.0, source: Any = 0.0, stabilization: str = 'galerkin', order: int = 6, unusual_parameters: UnusualParameters | None = None) -> ScalarTransportOperators
Integrate conservative skew transport and Galerkin, SUPG or UNUSUAL residuals.
Trial strong residual is -div(K grad(u))+beta.grad(u)+(c+div(beta))*u. SUPG tests against tau beta.grad(v); UNUSUAL subtracts tau times the complete reaction-diffusion residual pair with zero advection. The weak skew operator has c+div(beta)/2 and requires that sampled energy coefficient nonnegative. Matrices contain only volume and an explicitly supplied numerical-normal advection callback. Boundary enforcement, trace maps and kernels are caller declarations. Material-cut quadrature does not resolve interfaces for strong residual methods; those interfaces must fit the local approximation mesh.
Coefficient and inverse-inequality contracts for scalar UNUSUAL stabilization.
UnusualParameters
dataclass
Declare elementwise bounds for the reaction--diffusion residual method.
diffusion_lower bounds the smallest eigenvalue of A from below and
reaction_upper bounds the nonnegative reaction from above throughout
each fine cell. Scalars are global bounds; callables receive fine-cell
centroids and must return bounds valid on the whole corresponding cell.
Smooth coefficient callbacks require these declarations. Constant and
resolved Cartesian materials admit bounds obtained directly from their data.
inverse_constant optionally prescribes m in Eq. (15) of Martins,
Pacazuca and Valentin (CILAMCE 2025). It must satisfy 0<m<=1/3 and the discrete
elementwise coercivity check in :func:unusual_scale. Without an explicit
value, a conservative operator-inverse value up to 1/3 is computed. Thus constant-A
P1 uses exactly the published m=1/3. Higher-degree and variable-A constants
depend on the physical element and operator, not only the degree.
__post_init__() -> None
Reject inverse parameters outside the stated admissible interval.
unusual_scale(basis: FloatArray, gradient: FloatArray, diffusion_residual: FloatArray, weights: FloatArray, diameters: FloatArray, tensor: FloatArray, reaction: FloatArray, centers: FloatArray, parameters: UnusualParameters) -> FloatArray
Compute elementwise Eq. (15), including its zero-reaction limit.
Arrays have leading axes (cell, quadrature point); gradients additionally have (basis, dimension), and diffusion_residual is div(A grad(phi)). Bounds are checked at quadrature points. The automatic inverse parameter obeys mh²||div(A grad(v))||² <= A_min²*||grad(v)||² on the local polynomial space, with both norms defined by the supplied positive integration rule. Adequate quadrature and valid bounds between samples remain caller obligations. An explicit m is accepted only when the negative residual is strictly smaller than the physical diffusion--reaction energy for every polynomial on each cell (modulo constants when reaction vanishes). This generalized eigenvalue check admits stable explicit choices beyond the conservative automatic inverse bound without claiming those bounds are identical.
Eq. (15) is evaluated as mh²/[max(sigma_maxmh²,2A_min)+2*A_min], avoiding division by zero when sigma vanishes. This is the continuous diffusion limit of that equation, with no artificial reaction floor.
Polynomial and oscillatory skeleton spaces for frequency-domain Helmholtz.
PolynomialNeumannTrace
dataclass
Prescribe physical outward flux q.n by declared Legendre coefficients.
Coefficients multiply L_j(2*t-1), with t increasing from the first to the
second vertex in mesh.faces[face]. They are complex scalar data on one
complete straight macroface, independent of its subdivisions. No normal
sign is inferred from the coefficient vector: the values already denote
physical outward q.n on that exterior face.
This declaration preserves its polynomial degree under trace refinement and exact nested restriction. A generic callback has no declared degree and is instead projected using the executed boundary quadrature. Nonzero coefficients, however small, are never removed when determining degree.
__post_init__() -> None
Require finite one-dimensional data and trim only exact trailing zeros.
evaluate(parameter: Any) -> np.ndarray
Evaluate the declared polynomial at oriented macroface coordinates.
coefficients_on(space: FaceSpace) -> np.ndarray
Transfer declared data exactly in degree to every target face segment.
Affine pullbacks are computed as polynomial compositions. Coordinates above the declared degree are structurally zero, including when lower coefficients have rounding error. No tolerance discards a mode. Nonconstant data require ordinary polynomial trace spaces of adequate degree; constants use the actual basis's constant representation.
OscillatoryFaceSpace
dataclass
Bases: FaceSpace
Real coordinates for the complex exponential space of Section 6 (2020).
On each segment, an even degree 2k uses a constant and k cosine/sine
pairs; an odd degree 2k+1 additionally contains the linear polynomial.
Frequencies are omegalengthcos(n*pi/(2k+2)). A deterministic weighted
QR changes coordinates within this exact span. Its executed matrix is
available as transforms for persisted coefficient-vector contracts.
continuous is unavailable because the published space is discontinuous.
__post_init__() -> None
Validate the frequency and compute independent segment basis transforms.
from_executed_transforms(breaks: tuple[float, ...], degrees: tuple[int, ...], transforms: Sequence[Any], *, wave_number_length: float, continuous: bool = False) -> OscillatoryFaceSpace
classmethod
Restore owned binary64 segment coordinate matrices without recomputing QR.
Matrices multiply the declared raw columns on the right. They must be finite, square, nonsingular real binary64 arrays of size degree+1; their input bytes are copied and protected from later mutation. Face orientation is the direction of increasing parameter, including reversed evaluation at 1-t. Evaluation and the representation of one use these exact executed matrices. Callers verify their persisted digest and associate the same partition, degree, frequency and mesh orientation with the archived coefficient vector before restoration. This reconstruction does not reassert the current QR's conditioning test or the approximation hypotheses for an arbitrary supplied basis.
evaluate(parameter: Any) -> FloatArray
Evaluate the executed real basis at oriented face parameters.
constant_coefficients() -> FloatArray
Represent one in the executed segment bases without a fit.
helmholtz_skeleton(mesh: Any, omega: float, *, degree: int = 1, subdivisions: int = 1, oscillatory: bool = False) -> SkeletonSpace
Build two-component real/imaginary coordinates for a complex scalar trace.
The oscillatory frequencies use the physical edge length and the supplied angular frequency. They correspond to the homogeneous unit-wave-speed experiments in the article; no effective heterogeneous speed is inferred.
Physical face moment rules on triangular and polygonal macrofaces.
face_moment_rule_3d(mesh: Any, skeleton: Any, face: int, order: int) -> tuple[FloatArray, FloatArray, FloatArray]
Return physical moment rules, integrating the degree-2ell face mass exactly.
Physical normal moments of already assembled vector boundary functionals.
boundary_normal_integral(skeleton: Any, boundary: FloatArray, *, absolute: bool = False) -> float
Apply the hydrostatic trace to the assembled displacement boundary moments.
This uses exactly the same quadrature as the global equations; a separate
integration rule would make the finite-compressibility identity inconsistent.
Compensated accumulation reduces cancellation between opposite boundaries
on platforms with or without a wider native real type.
absolute sums uncancelled moments for scale-invariant compatibility.
require_balanced_moment(moment: float, absolute_terms: float, *, rtol: float = 1e-10) -> None
Require a signed integral vanish relative to its uncancelled physical terms.
This condition is independent of amplitude and physical units. Both inputs must come from the same integration rule. It supports flux compatibility, force balance and moment constraints without an absolute zero threshold.
require_compatible_displacement_flux(flux: float, absolute_moments: float) -> None
Require volume displacement balance in incompressible elasticity.
Globally continuous scalar pressure traces on planar macro edges.
Vertices shared by different macrofaces share coordinates. This differs from facewise C0 multiplier spaces, whose separate macrofaces remain independent.
PressureTraceSpace
dataclass
Globally continuous nodal polynomials on triangular or polygonal edges.
faces may prescribe different positive degrees and segment partitions.
Endpoints of different macrofaces share the mesh vertex unknown. Interior
edge nodes belong to that face only. This is Gamma, not the flux skeleton.
size: int
property
Return the number of shared pressure-trace unknowns before boundary elimination.
__post_init__() -> None
Assign shared vertices and private face-interior interpolation nodes.
uniform(mesh: TriangleMesh | PolygonMesh, degree: int = 1, segments: int = 1) -> PressureTraceSpace
classmethod
Use the same continuous nodal degree and subdivision on every face.
cell_dofs(cell: int) -> IntArray
Return unique Gamma indices touching one macrocell, in increasing order.
Componentwise physical data in an explicitly signed vector triangular trace.
vector_boundary_components_3d(skeleton: TriangularSkeleton, dirichlet: Any, traction: dict[int, Any], components: dict[int, dict[int, Any]], order: int) -> tuple[FloatArray, dict[int, float]]
Project physical outward traction and assemble all remaining velocity moments.