The half space you do not have to solve for

A boundary element method has one genuine weak spot, and it is elastic stress within about one element of a boundary it solved for. With piecewise-constant density that error is 10–45%, and it does not go away by shrinking the mollification width, because the production value already sits at its optimum there. It is a property of having discretized the boundary at all.

For a homogeneous half space with no topography, that boundary does not exist. The free surface lives in the Mindlin kernel. The half space is laterally and vertically infinite, so there are no sides and no base. There is no material interface. And fault slip is prescribed data rather than an unknown. So in this configuration nothing is solved: no linear system, no collocation, no free term, no boundary density, no box truncation.

What changes, laid out as what each error actually is rather than as a single number on both sides:

Those two rows are not the same quantity, and saying so matters: the left column is an error in a stress a reader would use, the right is the size of a traction that ought to vanish on a surface. The second is the honest statement of this kernel’s accuracy class, and the section below says why it cannot be made zero.

Mollification, and the Mindlin image

Following Cortez, the point force is not softened by hand — it is replaced by a smooth blob of width ε,

and the elastostatic equations are then solved exactly for that body force:

The result is singular nowhere and still satisfies a PDE, which hand-softening does not. A fault therefore has a finite width ε, and stress is defined on the fault, where classical dislocation theory returns a value ordinary integration does not define.

The half space adds Mindlin’s image terms at the mirrored source. Writing $R_\varepsilon^2 = |\mathbf{d}|^2 + \varepsilon^2$ for the mollified distance to the image and $d_3 = x_3 + y_3$ for the summed depth, every image term is rational in

plus a single $\log Q_\varepsilon$. That split by radical is what organises the whole implementation, because the two families integrate over a triangle very differently:

• the R-family closes in closed form, reusing the full-space moment hierarchy on the reflected triangle — the reflection is an isometry, so no new machinery is needed for it.
• the Q-family carries the $(1-2\nu)$ prefactor and is evaluated by adaptive Gauss quadrature under the budget law below. It stays that way on purpose, for a measured reason given at the end of the next section.

Every mollified jump carries an eigenstress

Mollifying a displacement jump smears it over a width ε, and a smeared jump is an inclusion with an eigenstrain. The raw mollified stress therefore diverges like $1/\varepsilon$ on the source surface — the analytic on-fault peak is $\tfrac{3}{4}\,\mu s/\varepsilon$ — and that divergence is the eigenstress the mollification put there, not a failure of the integration:

So stress_matrix subtracts it and returns the elastic stress by default. total_stress_matrix and eigenstress_matrix are the other two views, and total − eigenstress == elastic is an algebraic identity the gates assert at 1 ulp rather than to a tolerance — if it ever needed a tolerance, the three entry points would not be three views of one computation.

There is deliberately no bare strain_matrix. The familiar idiom is strain_matrix then strain_to_stress, and for a mollified kernel that path silently returns the total stress: the single failure this package exists to prevent. The strain readout is named elastic_strain_matrix, so the choice is visible at the call site.

What that buys, and what it does not

What it buys. On-fault stress is at machine precision — 1 × 10⁻¹⁵ at every realistic collocation point, including the top element of a fault that breaks the surface, which is the configuration that motivated the whole adaptive quadrature rule. That figure is agreement between the shipped kernel and an independent reference; it is not a statement about the free surface, and the two are worth keeping apart.

What it does not. The free-surface traction condition is satisfied to O(ε²), not exactly: 3.3 × 10⁻³ at ε/h = 0.1, 3.3 × 10⁻⁵ at ε/h = 0.01, with the order in ε measured at 1.90–1.96. Classical Mindlin gives $t_{3k} = 0$ at $z = 0$ identically; this kernel gives it only to second order, and the residual is worst as $\nu \to \tfrac{1}{2}$ — roughly twenty times larger at ν = 0.49 than at ν = 0.25. That is structural rather than incidental: the image terms carrying the $(1-2\nu)$ prefactor are exactly the $1/Q_\varepsilon$ and $\log Q_\varepsilon$ family whose mollification is the approximate part.

An exact free surface was derived, measured, and foreclosed. The route was a convolution rule applied to Papkovich–Neuber potentials, and it is exact only for harmonic potentials; for merely biharmonic ones it is O(ε⁴). Mindlin’s potentials are biharmonic in the source variable — verified to 10⁻¹⁶ — so the rule improves them to fourth order, but the traction takes two derivatives, and differentiating an error whose spatial scale is ε twice costs ε². O(ε⁴)/ε² is O(ε²) again. The arithmetic is consistent and no better implementation of that route escapes it, because the Cortez blob has algebraic tails: some of the blob’s mass sits above $z = 0$, and its image lands back inside the body.

There is a second recorded negative result, because it decides why the Q-family is quadrature. The Q-family can be rationalised — $1/Q^q = (R + d_3)^q / P^q$ — and the resulting closed form was built and is correct, verified at 10⁻⁴³ to 10⁻⁶¹ in sixty-digit arithmetic. It is also unusable in float64 everywhere in the physical domain. Because the body is $z \le 0$, the quantity $d_3 = x_3 + y_3$ is never positive, and that is precisely the cancelling sign: $R - |d_3| = P/(R + |d_3|)$ is a difference of nearly equal numbers. The best relative residual anywhere is 2.6 × 10⁻⁷, against the shipped quadrature’s 10⁻¹⁵. A closed form eight orders worse than the quadrature it would replace is not a speed-up; it is a regression with a speed-up attached.

Quadrature where it must be, closed form where it can be

The Q-family’s quadrature order is not a knob a caller tunes. It follows a measured budget law,

with $L$ the triangle’s scale, $\delta$ the observer’s distance to the image triangle, and $p$ the patch order, because a higher-order shape function raises the integrand’s polynomial degree. $C$ measured 5.8 to 8.4 for a relative accuracy of 10⁻⁹, and the shipped default is 10 — chosen with headroom over the maximum, not the mean, because starving this rule is silent. An earlier version used a flat 16, picked from a buried element where the error happens to be ε-independent, and that left an on-fault collocation point at 7 × 10⁻⁴ and a readout 0.03 h below a surface trace at order one.

Two guards matter more than the law itself:

• The floor of 8 is load-bearing. The law under-predicts far away — it asks for one point at sixteen triangle-lengths — so on a realistic matrix the floor is carrying 99.7% of all quadrature points, and with it the whole of the far-field accuracy. A flat floor of 4 fails at 1.9 × 10⁻³.
• The ceiling of 192 is reported, not silently accepted. Past it the Q-family would need a closed form rather than more points, and the cost is quadratic. The paragraph above explains why that closed form is not coming.

Beyond ten triangle-lengths the closed-form R-family is itself the worse choice — it loses digits as $(R/L)^{4\text{–}5}$ — so the producer switches to quadrature there. On a compact mesh that branch is never taken, which is worth knowing before timing anything: the far-field fraction of a mesh whose diameter is smaller than ten elements is exactly zero.