 |
BELFEM 0.9.0
Berkeley Lab Finite Element Framework
|
|
BELFEM 0.9.0
Berkeley Lab Finite Element Framework
|
Date: 2026-08-24 Purpose: Derive a dimensionally consistent Coulomb-gauge penalty for the h-formulation, establish where it can and cannot act in BELFEM's edge-element spaces, and record what the quench-deck conditioning numbers actually measure. Module: fem/maxwell, fem/interpolation/nedelec Status: Theory plus measured evidence. The gradient operator G is implemented and tested on every edge-function class except EF_LINE3; the penalty is implemented in the h-kernels behind the deck key nonlinear magnetic { coulomb gauge penalty { chi } }. It is opt-in since 2026-09-01: an absent block means chi = 0 and the term is not assembled; a deck switches it on by stating a positive chi (the value in cl_IWG_Maxwell.cpp is only the fallback for paths that never reach Controller::set_params). It was on by default at chi = 1e-4 from 2026-08-27 to 2026-09-01; on the tape decks that value was invisible to the conditioning estimate, which is why the default went back to off.
A Coulomb-gauge penalty on the h-formulation is meant to regularize the curl null space of the discrete operator. The results below say where that works and where it cannot:
Conventions follow Messe et al. 2023 (Eqs. 2–8) and the code in src/fem/maxwell/matrices/mt_maxwell_h.cpp (anchors are greppable tokens, not line numbers).
| symbol | meaning | code counterpart |
|---|---|---|
| h | magnetic field strength (edge dofs = circulations along edges, unit A) | dof vector q |
| b = μh | magnetic flux density | none |
| ρ | electric resistivity (scalar; HTS power law ρ(|j|,T,B,θ)) | mx->compute_rho |
| μ | magnetic permeability (μ0 in HTS/metal; μ(|h|) in ferro-conductors) | mx->compute_mu |
| e | electric field | boundary term only |
| j = ∇×h | current density | mx->compute_j |
| δh, v | test function (same edge space) | none |
| w_k | edge basis function k | column k of E |
| E | interpolation operator, d×n, h(x) = E·q | Calculator::E(k) / EF_*::E |
| C | curl operator, 3×n (3D) or 1×n (2D) | Calculator::C(k) / EF_*::C |
| G | gradient operator (§8) | Calculator::G(k) / EF_*::G |
| M | mass matrix ∫ μ EᵀE dV | aMatrices->M() in h_picard |
| K | stiffness ∫ ρ CᵀC dV | aMatrices->K() |
| Δt | timestep (mDeltaTime) | h in BDF comments |
| α, βᵢ | variable-step BDF coefficients | mAlpha, mBeta |
| qhist | β-weighted dof history | collect_qhist |
| χ | dimensionless penalty gain | deck key chi, penalty( 2 ) |
| γ | dimensional penalty coefficient (§6.2) | none |
| ρ* | reference resistivity (element, block or global) | save_resistivity |
| n | number of edge dofs per element | mNumDofs |
| d | spatial dimension | none |
The h-formulation weak form, Messe et al. 2023 Eq. (6):
∫_Ω δhᵀ ∂(μh)/∂t dV + ∫_Ω (∇×δh)ᵀ ρ (∇×h) dV + ∮_Γ δhᵀ (n×e) dS = 0.
As implemented (mt_maxwell_h.cpp, kernels h_picard, h_newton_mu0, h_newton_mu):
M += EᵀE·(μ·w·dV) K += CᵀC·(ρ·w·dV)
Timestep assembly (cl_IWG_Timestep.cpp, bdf1…bdf5; cl_TimestepMatrices.cpp assemble_dJdx): the element fixed-point system is
A = α·M + Δt·K, b = Δt·f + M·qhist,
and the Newton tangent adds dJdx = α·dMdX_times_x − dMdX_times_h + Δt·dKdX_times_x − Δt·dFdX into a second global matrix.
Solve structure (cl_FEM_DofMgr_SolverData.cpp, SolverAlgorithm::NewtonRaphson branch): Picard solves A·q_new = b directly with relaxation ω; Newton computes the residual r = A·x − b at the current iterate x, solves (A + dJdx)·Δ = r, and updates x ← x − ω·Δ. Two global matrices exist: the fixed-point matrix (used for residual evaluation) and the tangent (used for the increment solve).
Assembly placement consequence used below: anything added into the element K-position is automatically Δt-scaled into A and enters the residual through A·x; anything in the M-position is automatically α-scaled into A and contracted with qhist on the RHS. A penalty that is linear in h with frozen coefficients needs no dJdx channel and no explicit RHS term.
Taking the divergence of Faraday's law gives ∂(∇·b)/∂t = −∇·(∇×e) = 0: the continuous evolution preserves ∇·b = 0 exactly if the initial condition satisfies it (Monk 2003, §1.2 makes the same point for the Maxwell system: the divergence conditions are consequences of the curl equations plus initial data).
Discretely: the space of gradients of the scalar (nodal) space is a subspace of the Nédélec space, ∇S_h ⊂ X_h, and for a simply connected domain with connected boundary the null space of the discrete curl is exactly N_h(curl) = ∇S_h (Monk 2003, §7.2.1, where that hypothesis is stated). On a multiply connected conductor the kernel is ∇S_h plus discrete harmonic (cohomology) fields; §4.3 shows the argument below covers those as well. Testing the discrete h-formulation with v = ∇ξ_h kills the K term (∇×∇ξ = 0), leaving d/dt ∫ μ h·∇ξ_h dV = 0. The Galerkin scheme conserves the weak divergence of b against all discrete gradients, step by step, up to solver tolerance. This is the discrete analog of the continuum statement, and it is the reason the h-formulation "does not require any choice of gauge" in the conductor (Arsenault et al. 2021, introduction; the same paper notes that first-order curl elements satisfy ∇·B = 0 only element-locally, with the deficit in the normal jumps, and that transient runs keep the property when the initial condition is solenoidal).
Consequence: a penalty cannot restore a lost conservation property, because none is lost. Its only legitimate purpose here is spectral: to move the curl-null-space eigenvalues (currently set by α·M alone) closer to the curl-range eigenvalues (set by Δt·K), reducing the condition number. That reframing matters, because it makes the success criterion measurable: the penalty must raise the gradient-mode diagonal without perturbing the physical (curl-range) modes.
The classical construction (Bíró & Preis 1989, for the A-formulation) adds γ(∇·u, ∇·v) to a curl-curl form to ellipticize it. Monk 2003 §7.4 ("The ellipticized Maxwell system") analyzes exactly this and proves two things:
Edge elements exist precisely to avoid that pathology (Monk, Ch. 5 introduction). Monk §7.4 then lists the edge-element-compatible substitutes: a weak discrete divergence ∇_h·, which requires a mass-matrix inversion made local by mass lumping (Eqs. 7.43–7.44), or explicit divergence enforcement via a Lagrange multiplier. The natural gauge for edge spaces is the tree–cotree gauge (Dular et al. 1997 §III: "natural gauge condition in the space of edge finite elements, which requires the construction of a tree").
The operative question for BELFEM is therefore not how the Coulomb penalty is defined, but whether an element-local penalty built from ∇h does anything useful in these edge-element spaces. §6–§10 answer that quantitatively.
Members of ∇S_h are gradients of continuous piecewise-linear scalars: element-wise constant vector fields. The cohomology/cut fields that join the kernel on multiply connected conductors are element-wise constant too: On a single tet, every Whitney ∇w is antisymmetric. A Whitney-1 field with zero curl therefore has no antisymmetric gradient part, so ∇h = 0 on that element and the field is constant there. Tangential continuity glues these per-element constants into a field that fails to be a global single-valued gradient around a non-contractible loop, but inside every element ∇h ≡ 0 still holds. BELFEM's h-φ cuts are moreover largely φ-side jumps plus multiplier rows (Messe et al. 2023, §2), untouched by any volume term.
Any element integral whose integrand is built from the pointwise gradient ∇h (the full Jacobian, its trace ∇·h, ∇·(μh) with element-wise constant μ, or any quadratic form in these) evaluates to zero on such fields. Hence on TET4/TRI3 and the thin-shell kernels (§12) no element-local ∇-based penalty can raise the curl-null-mode eigenvalues at all, cohomology modes included. The discrete divergence content of an edge-element field lives in the normal jumps across element faces (Monk §5.5.1: "there is a singular contribution to the divergence at the faces in the mesh"), reachable only by face/jump terms, a weak (mass-inverted) divergence, or graph-based gauging.
Scope caveat: this element-wise-constant argument is a simplex/Whitney property. The lowest-order hex kernel contains non-constant curl-free fields (∇(xy) = (y,x,0) ∈ Q_{0,1,1}×Q_{1,0,1}×Q_{1,1,0}, curl-free, ∇h ≠ 0), so on HEX8 meshes a full G′G penalty does reach part of the kernel. §10 gives what that costs, and §11 gives why it still does not move the condition number at usable χ.
Order-of-magnitude entries for mesh size h_m (E ~ 1/h_m, C ~ 1/h_m², dV ~ h_m^d, shown for d = 3):
Their ratio r = Δt·ρ/(α·μ·h_m²) is the gap a gauge could close. Concretely, with μ0 = 1.26e-6 and α = 1:
| h_m | Δt | ρ | r |
|---|---|---|---|
| 1 mm | 1e-4 s | 1e-6 Ω·m (Hastelloy-class) | ~80 |
| 1 mm | 1e-4 s | 2e-9 Ω·m (Cu at 77 K) | ~0.16 |
| 1 μm | 1e-4 s | 1e-6 Ω·m | ~8e7 |
| 1 μm | 1e-2 s | 1e-6 Ω·m | ~8e9 |
Nothing in this range produces 17 decades. A perfect gauge, whether tree–cotree or an ideal weak divergence, closes at most this factor. That bounds what any gauging work can be expected to achieve on the condition number.
Note the direction of the Δt dependence: if λ_max ~ Δt·λ_K and λ_min ~ α·λ_M, then κ ∝ Δt. Collapsing the timestep during a quench shrinks this particular split rather than widening it.
A near-critical quench run (I ≈ 1.1·Ic, 97 completed steps) with the conditioning diagnostic on in the linear magnetic block, κ supplied by MUMPS from its own error analysis. (Measured before 2026-08-30; the MUMPS numbers are now behind their own key, mumps error analysis — under the current deck contract compute conditioning alone yields the eigenvalue ratio instead, a different quantity.)
| quantity | value |
|---|---|
| range | 1.1e17 … 7.4e19 |
| median on steps converging in ≤ 6 iterates | 5.9e17 (n = 85) |
| median on steps needing ≥ 12 iterates | 3.0e17 (n = 5) |
Two results matter more than the magnitude:
The run actually stopped on a residual plateau at −90…−96 dB. Measurements excluded timestep, linear solver, linear precision, and conditioning as causes; loosening the nonlinear tolerance from 1e-11 to 1e-9 was what got past it. A term that only improves κ therefore has no measured problem to solve on that deck. What a gauging design must answer first is whether that residual plateau lives in the gradient null space of the curl-curl operator, a question settled by projecting a converged residual onto the discrete-gradient / cotree subspace, not by tuning a penalty.
κ itself needs careful interpretation: on the MUMPS path the reported number is COND1 (RINFOG(10)), a componentwise 1-norm condition estimate of a nonsymmetric mixed system — computed on the original matrix with its actual right-hand side, not the equilibrated one: MUMPS scales for the factorization (ICNTL(8) = 77), but its error analysis hands the condition estimator an identity weight vector (MUMPS 5.9.1 dsol_driver.F, "Notice that D is always the identity"), so no scaling enters the estimate. (An earlier version of this paragraph said "already-scaled"; corrected 2026-08-30.) The companion quantities from the same factorization (COND2 and the forward and backward errors) say more about whether a solve is actually losing the solution than COND1 does.
Target functional (Coulomb gauge on b): Π = ½ ∫ γ |∇·(μh)|² dV, γ > 0 to be determined. First variation against the edge test space:
P(h, v) = ∫ γ (∇·(μh)) (∇·(μv)) dV.
Properties: (i) symmetric positive semi-definite; (ii) consistent: the exact solution has ∇·(μh) = 0, so the penalty vanishes on it and adds no modeling error in the continuum limit; (iii) linear in h for h-independent μ and γ.
It belongs in the stationary part of the equation (the K-channel), because the constraint ∇·b = 0 is an algebraic condition on the state at t^{n+1}, not on the rate. The rate variant is examined in §6.4.
Bookkeeping in SI with edge dofs in A (E ~ 1/m, G ~ 1/m², dV ~ m³):
| object | entry units | ×dof (A) |
|---|---|---|
| M = ∫μEᵀE | μ·m = V·s/A | V·s |
| K = ∫ρCᵀC | ρ/m = Ω = V/A | V |
| P = ∫γμ²Gdivᵀ Gdiv | γ·μ²/m | none |
Matching P to the K-position requires γ·μ²/m = V/A, i.e.
γ = χ · ρ* / μ², χ dimensionless, ρ* a reference resistivity,
giving P(h,v) = χ·ρ*·∫ (∇·h)(∇·v) dV for element-wise constant μ: the μ² cancels pointwise, so nothing is gained by carrying it. χ then has the clean meaning penalty resistance per physical resistance.
The ρ* fork. The choice of averaging scope is a genuine design fork, not a detail:
The implemented kernels use the live pointwise ρ, namely the first branch. A floored coefficient must not be combined with the full G′G form: on TET4 that is an exact resistivity injection Δρ = (χ/2)·ρ_floor on every curl mode, which turns subcritical HTS into a poor metal and silently undoes the unfloored power law. A floored coefficient with a divergence-only operator avoids that, but the divergence row is identically zero on most of BELFEM's element families (§10), so it is close to the zero operator on the decks of interest.
A natural-looking construction adds, per element and Gauss point, α·[ μ²(G′G) + μ E′E dH/dt ]·w·dV with α = χ·μ/Δt and χ dimensionless. Both terms fail, independently:
Term 1, α·μ²·G′G. Units: (χμ/Δt)·μ²/m = χ·μ³/(Δt·m) = χ·V³s²/(A³m⁴) per entry. The K-position needs V/A; the M-position needs V·s/A. Neither matches, so χ would have to carry units A²m⁴/(V²s²). With dimensionless χ this is not a valid SI formula.
Term 2, α·μ·E′E·dH/dt. Two independent problems. Dimensions: entry·rate = χμ²·m·A/(Δt·s), again neither V nor consistent with Term 1. Structure: E′E is the mass Gram; it contains no derivative of h at all. A term ∝ E′E·dH/dt is an artificial addition to the permeability, not a divergence penalty of any kind: it acts identically on curl-range and curl-null modes, so it cannot selectively fix the gradient subspace, while it does shift the physical eddy-current time constant.
Where the form comes from. Using the identity ∫|∇u|² = ∫|∇×u|² + ∫|∇·u|² + boundary terms, one may try to build a divergence penalty as "full-gradient penalty minus curl-curl content", then eliminate the curl-curl part through the strong equation ∇×(ρ∇×h) = −∂(μh)/∂t, which turns it into a mass-rate term. That produces exactly the shape above. It fails twice more: the substitution is valid only at the converged solution, so inside a Newton loop it changes the tangent and the transient; and on simplices the discrete identity degenerates, because ∫∇wₐ:∇w_b = ½∫curlₐ·curl_b holds exactly on TET4. "G′G minus curl content" is then −½·C′C, a negative curl stiffness carrying no divergence information.
Choice of channel. In BELFEM's assembly path, the channel determines whether a contribution belongs in the matrix or the RHS: a coefficient times E′E·dH/dt is a mass-channel term and splits automatically (α-part into the matrix, qhist part into the RHS); a stationary constraint G′G·h is K-channel, matrix only, Δt-scaled and residual-consistent through r = A·x − b. Nothing should ever be added to the RHS by hand; the BDF contract owns that split.
Penalizing ∂(∇·b)/∂t (G-terms contracted with dH/dt) enters the M-channel and only prevents divergence growth: any divergence content present in the initial state or introduced by a restart stays frozen. Combined with §4.1 (the weak divergence is already conserved) the rate form has no target left.
For P = χ·ρ*·∫(∇·h)(∇·v) dV, or the full-gradient variant, with ρ* frozen during the nonlinear loop:
G is the full Jacobian of the interpolated field, not only its divergence.
At integration point k, G(k) is the (d·d) × n matrix with
G( i + d·j , e ) = ∂(w_e)_j / ∂x_i i, j ∈ {0…d−1}, e ∈ {0…n−1}
(column-major vec of the d×d tensor ∇h with convention (∇h)_{ij} = ∂h_j/∂x_i). Shapes: 9×n in 3D, 4×n in 2D. The binding layout contract, including the explicit curl-tie rows, lives at the mGrad member in src/fem/interpolation/nedelec/cl_EF_EdgeFunction.hpp; Calculator::G( aIndex ) is the assembly-side accessor.
Derived quantities:
Why full G and not a 1×n divergence row:
Cost: everything G needs (the per-point inverse Jacobians and reference derivatives) is already computed for E and C; G adds O(d²·n) multiplies per Gauss point, the same order as one E evaluation. EF_LINE3 deliberately does not implement G: on a 1D manifold element, the ambient gradient is not defined from the element alone. Its G() is a permanent stub that always fails.
In h-φ, edge dofs exist only in conducting domains (air/vacuum and non-conducting iron carry φ; Messe et al. 2023 Eqs. 6–8). Therefore:
| family | dofs | geometry map | ∇·(basis) element-wise |
|---|---|---|---|
| TRI3 | 3 | affine | ≡ 0 (Whitney, 2D) |
| TET4 | 6 | affine | ≡ 0 (Monk §5.5.1 states it for the Whitney element) |
| TRI6 | 8 | straight or curved | ≠ 0 for the higher-order members (div ∇(quadratic) = const ≠ 0); the curved path adds mapping terms |
| TET10 | 20 | straight or curved | ≠ 0 for higher-order members; curved adds mapping terms |
| HEX8 | 12 | trilinear | box ≡ 0, sheared affine ≠ 0, general trilinear ≠ 0. But the full ∇h of kernel modes is ≠ 0 even on boxes (the ∇(xy) mode, §4.3) |
| PENTA6TS | 6 | flat prism, thickness τ | ≡ 0 (∇w antisymmetric-tangential; n·w = 0) |
| QUAD4TS | 2 | 2D thin shell | ≡ 0 only for orthogonal stacking: the basis is s·F(η)·∇ξ, so ∇·h = F′(η)·(∇η·∇ξ) ≠ 0 when the stacking is sheared. Kernel fields (equal top/bottom dofs, hence constant h) stay blind regardless |
| HEX8TB | 4 | side-connector wall, imposed orthonormal cuboid frame | ≡ 0 identically: the basis is s·F(η,ζ)·∇ξ and ∇η, ∇ζ ⊥ ∇ξ by construction |
| HEX8TS | 8 | thin-shell machinery, per-point mid-surface nablas | ≡ 0 on rectangular mid-surfaces; nonzero under sheared stacking, as for QUAD4TS |
| LINE3 | 2 | 1D manifold | divergence not meaningful on the curve |
Mechanism for the affine cases. Under the covariant Piola map h = J⁻ᵀĥ with constant J, the physical divergence is ∇·h = (J⁻¹J⁻ᵀ) : ∇̂ĥ: the contraction of the symmetric inverse metric with the reference Jacobian of the basis. Whitney simplex functions have antisymmetric constant ∇̂ĥ, so the divergence is exactly zero for every affine map. The hex reference functions have symmetric off-diagonal ∇̂ĥ content, so their divergence is zero only when the metric is diagonal (a box) and nonzero under shear. On hex meshes the penalty magnitude is therefore partly a measure of element shape, not of physics.
Implication. On TET4/TRI3/PENTA6TS/HEX8TB meshes the div-penalty matrix is identically the zero matrix; on QUAD4TS and HEX8TS it is zero for orthogonal stacking and shear-dependent otherwise. Per §4.3 the full G′G penalty cannot touch the curl null space on any of them, so a nonzero basis divergence under shear buys nothing. On mixed meshes the penalty switches on and off by element type, which is itself a conditioning hazard.
The assembled system is A = α·M + Δt·K. Two candidate scalings place the penalty on different sides of that sum, and adaptive stepping makes the difference operational.
Scaling the penalty as χ·μ/Δt on a K-channel term. The Δt from assembly cancels the 1/Δt: the assembled penalty becomes Δt-independent, i.e. proportional to the mass block. Then the penalty-to-K ratio goes as 1/Δt, so at large Δt the penalty fades relative to the physics, and at collapsed quench timesteps it grows, giving maximum distortion at maximum fragility. The gap it is meant to close, r = Δt·ρ/(α·μ·h_m²), is itself Δt-dependent, so a Δt-free penalty cannot track it.
The consistent scaling γ = χ·ρ*/μ² (§6.2), K-channel. Penalty and physical stiffness carry the same Δt factor. The penalty-to-K ratio is χ·(ρ*/ρ_local)·(div-content / curl-content), Δt-invariant across the adaptive range, and the gradient-mode diagonal lifts to the same Δt·ρ*/h_m² scaling as the curl modes.
What χ can and cannot do. It is easy to conflate two thresholds:
A null result at χ = 1e-4 is therefore expected on every element family, including the hexes where the term is not vacuous. A K-channel penalty can visibly move a near-1/ε condition number only if it fills exact zeros on the small-end subspace, which it does nowhere on BELFEM's current meshes. Combined with §5.1 (the kernel-versus-range gap is at most 1e2…1e8) and §5.3 (κ does not predict solver difficulty on the measured deck), χ tuning is not the lever it appears to be.
If χ is ever needed, validate it by measuring the quantity it is supposed to improve with and without the penalty on the same deck, not by argument alone.
In 2D (TRI3/TRI6; mE is 2×n, mC is 1×n scalar):
G is 4×n: rows (∂x h_x, ∂y h_x, ∂x h_y, ∂y h_y); ∇·h = row(0) + row(3); curl h (scalar, out-of-plane) = row(2) − row(1).
Penalty reduced form: P = χ·ρ*·∫ (∂x h_x + ∂y h_y)(∂x v_x + ∂y v_y) dA. The effectiveness verdict is identical to the 3D simplices: TRI3 divergence ≡ 0 element-wise and the curl null space is element-wise constant, so the penalty is vacuous on the 2D workhorse meshes. TRI6 carries nonzero divergence in its higher-order members only.
An element-local ∇h penalty should not be used as the conditioning strategy for BELFEM's simplicial or thin-shell quench systems: on those families it is either exactly zero or exactly an artificial resistivity, and on the families where it is non-vacuous it cannot reach the mass scale at any consistent χ (§11). The implemented term is retained, defaulted off, for hex-dominant meshes and experiments.
If conditioning or the residual plateau is to be addressed at the formulation level, the edge-element-compatible options are, in increasing order of invasiveness:
| option | mechanism | reference | remark |
|---|---|---|---|
| a. Measure first | attribute the number: which block owns ‖A‖, what the forward error is, whether the residual plateau lives in the gradient/cotree subspace | none | prerequisite to everything below; costs one instrumented run |
| b. Scaling / equilibration | MC64 and diagonal scaling on the assembled matrix | already load-bearing on both solver paths | zero formulation risk; on the MUMPS path the reported κ is already post-scaling |
| c. Weak-divergence penalty | γ(∇_h·(μh), ∇_h·(μv)) with a lumped-mass weak divergence | Monk §7.4, Eqs. 7.43–7.44 | acts on the true (jump) divergence; needs face-adjacent assembly, not just G. Blind to discrete harmonics: Monk's Lemma 7.26 splits E_h with ∇_h·(εẼ_{h,0}) = 0, so cohomology modes survive it |
| d. Tree–cotree gauge | eliminate the curl-null subspace in conductors | Dular et al. 1997 §III (tree gauge of the source field h_s; gauging the eddy-current unknown itself is the Albanese–Rubinacci-style use); Denis et al. 2026 | exact null-space removal, no scale threshold; interacts directly with cut/cohomology generators, and enriched meshes carry twin edges that share both end nodes, so a tree keyed on edge endpoints is wrong there |
| e. Full-G Tikhonov, hex/high-order meshes only | χρ*(∇h, ∇v) | §6.2 | inconsistent: perturbs the physics at O(χ·ρ*/ρ_local), fatal in SC zones unless support-restricted to quenched elements |
The G-operator is prerequisite for (c) and (e), useful for the diagnostics in (a) under every option, and independently valuable for postprocessing.
The dimensional algebra, the null-space arguments, and the per-family divergence mechanism are derivations. They are given inline above and reproducible from the element constructions in src/fem/interpolation/nedelec/. Four were also checked numerically at the mathematical level before implementation: the Whitney divergence identity on TET4, the Gram identity ∫∇wₐ:∇w_b = ½∫curlₐ·curl_b on TET4 (max deviation 7.1e-15), the vanishing of ∇h on reconstructed discrete-gradient fields (max 6.9e-16), and the box-versus-sheared HEX8 divergence contrast. The compiled implementations are gated separately in tests/fem/test_EdgeFunctions.cpp and tests/fem/test_InterfaceOrientation.cpp: finite differences of G against E, the antisymmetric part of G against the compiled C, the divergence trace against §10, and per-family kernel and control modes.
The conditioning statements in §5.3 are measurements from an instrumented quench run, not estimates. The scale bounds in §5.1 and §11 are arithmetic from those same material and mesh parameters.
Full citations and DOIs: doc/literature_references.md.
1.18.0