Visualisation: principal-stress glyphs and stress trajectories

[lmoresi]

What this adds

A way to see the stress tensor, in the same spirit as velocity arrows: sampled at seed points, not drawn everywhere.

• principal_stress_glyphs(coords, stress, scale) — one bar per principal axis at each seed: length proportional to the principal-value magnitude, with a "tensile" cell array so the two signs colour separately (blue compressive, red tensile, matching the RdBu_r field convention). A cross in 2-D; three orthogonal bars in 3-D.
• direction_trajectories(direction_at, seeds, inside, step, separation) — stress trajectories. A principal direction is defined only mod 180°, so ordinary streamline tools cannot integrate it; this integrator sign-aligns each evaluated eigenvector with the previous heading, and places lines evenly (Jobard–Lehmann occupancy with separate seed-blocking and line-stopping tests). 2-D only — in 3-D the analogue is a trajectory surface, which we do not attempt.
• tensor_fn_to_pv_points, trajectories_to_pv_lines — the evaluation and bundling steps, mirroring the existing *_fn_to_pv_points helpers.
• plot_stress_glyphs(mesh, stress, ...) — the one-call wrapper beside plot_vector. Default seeding is a grid over the bounding box filtered to points inside the mesh, so an annulus seeds nothing in its hole.

Docs page docs/advanced/stress-visualisation.md with worked figures (a blind-thrust fault network and a 3-D Stokes sinker, both rendered from checkpoints). It records two facts a user needs: the pressure datum is a gauge that can flip bar colours but never rotate axes, and World-Stress-Map-style regime colouring is degenerate in 2-D incompressible plane strain (the out-of-plane stress is always the intermediate principal stress), so it is deferred until there is a 3-D use case.

Tests

tests/test_0848_stress_glyphs.py (level_1, tier_a): glyph geometry for uniaxial compression, pure shear, and a 3-D diagonal tensor; the eigenvector sign-flip case that a naive integrator fails by reversing mid-line; trajectory separation; and the annulus default-seeding case end-to-end through plot_stress_glyphs. All pass in the worktree environment; the style gate is clean.

Underworld development team with AI support from Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[lmoresi]

Adversarial review

We attacked the diff along the ways a glyph routine can silently draw the wrong thing.

Found and fixed (32b397d): with user-supplied seeds, the auto-scale estimated spacing as extent / num_seeds — a grid parameter the caller's seeds owe nothing to. A 13×13 section plane with the default num_seeds=24 drew bars at roughly half the intended size. Now estimated from the mean nearest-neighbour distance of the seeds themselves (subsampled above 2048 seeds to bound the pairwise matrix).

Checked and held:

• Eigen-decomposition trusts one triangle. numpy.linalg.eigh reads only the lower triangle, so an asymmetric input (recovered components are not exactly symmetric) would be half-ignored silently. The builder symmetrises first; the docstring says why.
• Mod-180° integration. A direction field's eigenvector sign is arbitrary; a naive integrator reverses mid-line and draws a folded stub. The test suite includes a field whose reported sign alternates underfoot (sin(20x) flip) and asserts the trajectory crosses the whole box without turning.
• Occupancy split. Using one occupancy set for both seed-blocking and line-stopping chops trajectories into dashes (we hit exactly this in the prototype). The landed integrator keeps a wide corridor for seeds and a traversed-cell set for stopping, and a line claims its cells only after integrating so it never blocks itself.
• Non-box domains. Default seeding on an annulus: the bounding-box grid is filtered by closest-cell distance, so the hole seeds nothing and evaluate is never called outside the mesh. Asserted end-to-end in test_annulus_default_seeds_avoid_the_hole.
• Units. tensor_fn_to_pv_points strips Pint magnitudes before PyVista sees them and stashes the units string, matching the scalar/vector helpers.
• Style gate clean; tests level_1/tier_a, 6/6 pass in the worktree environment.

Known limits, stated in the docs rather than papered over: closed trajectory orbits are traced once per integration sense (cosmetic overdraw); 3-D trajectories are deliberately out of scope (surfaces, not curves); the compressive/tensile split is relative to the pressure gauge, which can flip colours but never rotate axes.

[lmoresi]

Validation against the fault-interaction examples

We drove the landed API end-to-end over the checkpointed fault-interaction models (nested listrics, en-echelon, blind thrusts, all at w = 0.01) and the 3-D sinker: tensor_fn_to_pv_points evaluates the recovered-stress Matrix at the seeds, principal_stress_glyphs builds the crosses, direction_trajectories + trajectories_to_pv_lines draw the trajectory nets, and the sinker goes through the one-call plot_stress_glyphs with section-plane seeds. The figures reproduce the pre-PR prototypes exactly. No solving anywhere — everything loads from read_timestep.

The exercise caught two more defects, both fixed:

• 8f064dd — the docs recipe built the stress Matrix from Txx.sym, but a scalar variable's .sym is a 1×1 Matrix; the recipe as written nested matrices. Now indexed (Txx.sym[0]) with a note saying why.
• 2479bb0 — plot_stress_glyphs promised a plotter callers could decorate, but ran show() unconditionally first, and a finalized scene ignores later actors: the sinker's box and sphere overlays were dropped silently. A show= parameter now leaves the plotter open (camera set) so overlays and screenshots work as documented.

Tests still 6/6 after both fixes.

[lmoresi]

CI: first run timed out — diagnosed and fixed (d69839d)

The first test run lost a worker at 78% (the intermittent VTK-render crash genre) and then hung from 90% until the 2-hour cap. Development is green on the same configuration, so the cause was in this branch: the annulus test drove plot_stress_glyphs through its default show=True, making it the first test in CI to call show() — which on a headless runner can enter an interactor wait and hang the xdist session.

The test now uses the show=False path, so nothing in the test file ever renders — the assertions are about what the plot builds (seed filtering, glyph geometry), not what it draws. This is the same discipline test_0847 established for plotter lifetime, extended to the render call itself.

[lmoresi]

Adversarial review

Reviewed at the PR head. Three findings, two checks that came back clean, and
one link to another open PR.

1. Degenerate points are the failure mode this integrator has, and nothing
covers them. direction_trajectories sign-aligns each evaluated eigenvector
with the previous heading, and test_trajectory_survives_eigenvector_sign_flip
pins that. A sign flip is not the hard case. Where the two principal values
cross — an isotropic point, generic in a 2-D stress field, and exactly what
occurs between the fault tips in the figures on the docs page — the eigenvectors
are not merely sign-ambiguous but arbitrary in orientation, and continuity with
the previous heading cannot recover a direction that is undefined. A trajectory
passing near one will turn wherever eigh happens to point.

That is inherent to the object rather than a defect in the code, and it is worth
saying so on the docs page next to the figures a reader will draw conclusions
from. A test would be a stress field with a known isotropic point, asserting the
integrator terminates or flags rather than producing a smooth-looking line
through it.

2. np.linalg.eigh ordering is load-bearing and undocumented in the glyph
function. The "most-tensile axis" convention depends on eigh returning
ascending eigenvalues. The tests encode the convention independently — the
oracle is the closed-form angle = 0.5 * arctan2(sxy, 0.5 * (sxx - syy)), which
is the right way to check it — but principal_stress_glyphs itself does not say
which end of lam is which, and a reader adding a third colour or reordering
bars has nothing to work from.

3. separation and step are given defaults (0.04, 0.008) in the docs
example with no statement of what they are relative to. Both are lengths in
model units, so the defaults are implicitly tied to a domain of order 1. On an
annulus of radius 6371e3, or a unit box scaled to metres, they produce either one
line or several million steps. The wrapper could derive them from the mesh extent
rather than carrying absolute numbers.

Checked and cleared.

• Symmetrisation is real. eigh reads one triangle, and the input is
  symmetrised at both entry points — 0.5 * (stress + np.transpose(stress, (0, 2, 1)))
  in the glyph path and the equivalent on stress_values. A non-symmetric input
  (a recovered stress that has not been symmetrised) will not silently produce
  eigenvectors of the wrong tensor.
• The :alt: options are single-line. Each is one long line rather than
  wrapped, which is what MyST requires — a wrapped :alt: truncates the alt text
  and renders the remainder as a stray second caption. Worth stating because the
  alt text here is unusually long and the temptation to wrap it is real.

Link to #602. tensor_fn_to_pv_points samples through
uw.function.evaluate, not global_evaluate, so #602's non-finite recovery does
not reach this code. If the trajectory integrator is ever run in parallel it will
want the global path, and at that point #602's substitution of an RBF
extrapolation for a failed interpolation becomes a substitution inside a figure —
worth remembering rather than acting on now.

Underworld development team with AI support from Claude Code