Two things classical dislocation theory cannot do

Classical elastic dislocation theory represents a fault as a surface of zero thickness across which displacement jumps. It is enormously useful and it has two defects that are not numerical but structural.

The first is that the stress on the fault itself is not defined. The displacement-discontinuity kernel is singular as the observation point approaches the source, and the integral over the element does not converge there. Codes return a number anyway — by taking a one-sided limit, by evaluating at an offset, by discarding the self term — and those numbers disagree with each other. For a quantity as consequential as the shear stress resolved on a fault that is about to slip, that is an uncomfortable place to be.

The second is that a fault of zero thickness is not a fault. Real fault zones are damaged, finite-width volumes. A representation that insists on zero width is not an approximation of a thin thing; it is a different object, and the stress it predicts diverges as you approach it rather than converging to anything.

Mollification

Both follow from the same source, so both yield to the same fix — but not the fix that first suggests itself. The obvious move is to soften the kernel directly: replace the distance $r$ wherever it appears in the Green’s function with

That does produce a smooth function. It is not a solution of anything, and it is not what Cortez did.

Cortez’s construction regularizes the source, not the kernel. Replace the point force $\delta(\mathbf{x})$ with a smooth, radially symmetric blob $\phi_\varepsilon$ of unit integral, and then solve the elastostatic equations exactly with that forcing. What comes back is a Green’s function that is smooth everywhere and satisfies a partial differential equation — the regularized Cauchy–Navier equation,

For that blob, the regularized Kelvin solution is

The first two terms are the classical Kelvin solution carrying $R_\varepsilon$ in place of $r$ — the naive substitution. The third term is the point. It has no singular counterpart, it vanishes as $\varepsilon \to 0$, and it is what the blob convolution contributes. Drop it and the residual in the regularized Cauchy–Navier equation is of order $10^{-2}$; keep it and the residual is $10^{-16}$ at every Poisson ratio, which the repository checks symbolically in verify_pde_residual. A kernel that merely looks smooth is not the same object as one that solves the equations, and only the second lets you integrate over the element the observation point sits on.

So $\varepsilon$ is not a numerical fudge. It has units of length and it is the width of the fault zone — the physical scale over which slip is distributed — and the stress at a point on the fault is now an ordinary integral with an ordinary value.

What makes this practical rather than merely tidy is that each triangle’s contribution is integrated analytically — in closed form, for a density that is constant, linear or quadratic over the element. There is no quadrature, so there is no quadrature error tying $\varepsilon$ to the element size $h$. The two are independent knobs: you can refine the mesh at fixed fault-zone width, or narrow the fault zone at fixed mesh, and the method does not object. This is the elasticity analogue of a result established for regularized Stokeslets.

Every mollified jump carries an eigenstress

Here is the part that is easy to get wrong, and it is the reason an on-fault stress from this method means something.

Smearing a displacement jump over a width $\varepsilon$ does not only soften the kernel. It introduces a genuine anelastic strain — an eigenstrain $\boldsymbol\varepsilon^{*}$, the distributed version of the slip. The stress you read off a mollified source is therefore the total stress,

and inside the fault zone the second term dominates: it peaks at $\tfrac{3}{4}\mu s/\varepsilon$ and diverges as $\varepsilon \to 0$. Reported as an elastic stress it would be badly wrong, and wrong in a way that gets worse the narrower you make the fault.

So the eigenstress is subtracted. Not approximately — the exact finite-triangle form is known in closed form, element by element, each with its own $\varepsilon$. What remains is the elastic stress, it is finite on the fault, and as $\varepsilon \to 0$ it converges to the finite part of the classical triangular-dislocation solution.

This is not only a fault story. Every mollified double layer carries an eigenstress, including a boundary patch, whose density is a jump between the field inside the body and zero outside. Subtracting it there too removed an error that had previously looked like a mesh limit near boundaries — an error that did not improve under refinement, because it was not a discretization error at all.

What that buys, and what it does not

“Interpretable” is not “exact”, and it is worth being precise about which is which, because the two limits are easy to confuse and they respond to different remedies.

Measured on a manufactured problem whose answer is known in closed form — a uniform strain $\mathbf{u} = \mathbf{A}\mathbf{x}$ imposed on a closed box, whose interior stress is exactly $\mathbf{C}\!:\!\mathbf{A}$ at any distance from the boundary — over a grid of element sizes $h$ and widths $\varepsilon$ given independently:

Read the first two rows: halving the mesh changes nothing, while $\varepsilon$ times 2.4 multiplies the error by 2.4.

What that floor is took four measurements to pin down, and the first three answers were wrong.

It is not the discretization: halving $h$ moves it by under 8 %. It is not the density order: P1 and P2 agree with each other to 0.5 %, and P1 represents $\mathbf{u} = \mathbf{A}\mathbf{x}$ exactly, so if the floor were the piecewise-constant density it would have collapsed. It is not the free term either: the calibrated diagonal and the analytic half-jump agree to 9 %. And it is not simply “the mollified problem is a different problem” — a symmetric unit-integral blob has zero first moment, so a linear field convolves to itself exactly, and a constant field indeed comes back at $3\times10^{-15}$ for every $\varepsilon$ and every order.

What the floor scales with is the domain. Holding $\varepsilon/h$ fixed at 0.145 and growing only the box, the error times $L/\varepsilon$ is 1.060, 1.091, 1.101 at $L$ = 40, 80, 120 km, while the same error times $h/\varepsilon$ moves by a factor of three. So the floor is $\varepsilon/L$, and one per cent asks for $\varepsilon \lesssim 0.01\,L$ no matter how fine the mesh.

The reading that survives: the smeared boundary itself. In the interior the blob is symmetric and preserves linear fields; at a boundary the body lies on one side only, the convolution is truncated, the cancellation fails, and the effective surface sits $O(\varepsilon)$ away from the nominal one — which costs a field with a gradient a relative $\varepsilon/L$. That last step is inference from the scalings rather than derivation. But it is what is left after order and free term were each ruled out by measurement, and it predicts the asymmetry below.

The third row is a different effect again, and here the asymmetry matters. Near a boundary whose density was solved for, the piecewise-constant representation of that density limits the readout and the error rises by a further factor of three to thirteen.

A fault is exempt, and for a stronger reason than it being prescribed data. Because the regularization is applied to the source, the fault’s slip is spread over $\varepsilon$ — and that spreading is the finite-width fault zone the method exists to represent. A point two $\varepsilon$ from the fault is reading the model, not an artefact of it. At a free surface the same spreading has no such warrant: the ground surface is not smeared over $\varepsilon$, so within that distance you are inside a layer the numerics invented. The same quantity, distance in units of $\varepsilon$, is physics on one side and an artefact on the other — which is why this site’s viewer thresholds on distance to boundaries only, and never hides the fault zone.

Neither distance alone governs it — binning the error on $d/h$ leaves a 4.2× spread across $(h, \varepsilon)$ and on $d/\varepsilon$ a 3.8× spread — so both are reported per voxel, and the Examples page lets you threshold on the first and watch the surface shell disappear. One practical consequence, measured: at fixed $\varepsilon$, refining $h$ by four improved the near-boundary error by only 1.27×, while standing off from $d/h = 0.2$ to $2$ improved it thirteenfold. Evaluate deeper; refining the patch alone barely helps.

Solving it at scale

The boundary integral equation is solved by collocation, which gives a dense operator — and at a million elements a dense operator is four million unknowns and hundreds of gigabytes. Two compressed far fields are available, and which one to want depends on what you are doing rather than on which is better:

• a flat H-matrix with adaptive cross approximation, whose rank adapts per block to whatever the geometry needs, and
• a Chebyshev black-box fast multipole method, whose cost is set by an interpolation order that adapts to nothing.

Measured against each other on the model on this site’s front page, the first assembles faster and applies its operator about 55 times more slowly; the second fits in a fraction of the memory. Near the million-element target the ACA operator is quicker but needs some 419 GB, while the FMM fits in 66 and is about 1.5 times slower. So the H-matrix is the choice when the operator will be applied many times — a rate-and-state sequence, an inversion — and the FMM is the choice when the problem would not otherwise fit. Both reproduce the dense LU reference on the four-state showcase to about 2 × 10⁻⁶, which is how we know the question is one of cost rather than of correctness.