diff --git a/docs/advanced/figures/stress_glyphs_enechelon.png b/docs/advanced/figures/stress_glyphs_enechelon.png new file mode 100644 index 00000000..06873043 Binary files /dev/null and b/docs/advanced/figures/stress_glyphs_enechelon.png differ diff --git a/docs/advanced/figures/stress_glyphs_listric.png b/docs/advanced/figures/stress_glyphs_listric.png new file mode 100644 index 00000000..9e479496 Binary files /dev/null and b/docs/advanced/figures/stress_glyphs_listric.png differ diff --git a/docs/advanced/figures/stress_glyphs_sinker3d.png b/docs/advanced/figures/stress_glyphs_sinker3d.png new file mode 100644 index 00000000..577c157a Binary files /dev/null and b/docs/advanced/figures/stress_glyphs_sinker3d.png differ diff --git a/docs/advanced/figures/stress_glyphs_thrust.png b/docs/advanced/figures/stress_glyphs_thrust.png new file mode 100644 index 00000000..6ee9196c Binary files /dev/null and b/docs/advanced/figures/stress_glyphs_thrust.png differ diff --git a/docs/advanced/index.md b/docs/advanced/index.md index bc723aae..cb47f8de 100644 --- a/docs/advanced/index.md +++ b/docs/advanced/index.md @@ -53,6 +53,13 @@ What a crack model cannot do at a branch point, measured in a shear box **[→ Branching Faults and Junctions](fault-branching-junctions.md)** +### Visualising the Stress Tensor +Principal-stress crosses sampled at seed points (2-D and 3-D) and +stress trajectories — the tensor equivalents of velocity arrows and +streamlines. + +**[→ Stress Visualisation](stress-visualisation.md)** + ### Custom Meshes Create complex geometries with gmsh for research problems. @@ -127,6 +134,7 @@ crossing-fault-zones gouge-zones fault-branching-junctions fault-mechanics-examples +stress-visualisation custom-meshes curved-boundary-conditions mesh-adaptation diff --git a/docs/advanced/stress-visualisation.md b/docs/advanced/stress-visualisation.md new file mode 100644 index 00000000..0c186620 --- /dev/null +++ b/docs/advanced/stress-visualisation.md @@ -0,0 +1,197 @@ +# Visualising the Stress Tensor + +Scalar fields get colormaps and velocity gets arrows; the stress tensor +needs its own glyph. `underworld3.visualisation` provides +principal-stress glyphs — sampled at seed points, the way velocity +arrows sample the velocity — and stress trajectories, the curvilinear +net traced by the principal directions. Both work from the same +recovered stress fields you already checkpoint. + +## Principal-stress glyphs + +At each seed point we diagonalise the stress and draw one bar per +principal axis: bar length proportional to the principal-value +magnitude, blue for compressive ($\lambda < 0$), red for tensile. In +2-D this is the classical stress cross; in 3-D each seed carries three +orthogonal bars, best drawn on one or two section planes rather than a +filled volume. + +```python +import underworld3 as uw +import underworld3.visualisation as vis +import sympy + +# Recovered stress components (P1 projections of the deviatoric +# stress) plus pressure give the full stress. Project once after the +# solve; do not pass raw solver derivative expressions to a plot. +# A scalar variable's .sym is a 1x1 Matrix — index it before +# assembling the tensor. +stress = sympy.Matrix([[Txx.sym[0] - P.sym[0], Txy.sym[0]], + [Txy.sym[0], Tyy.sym[0] - P.sym[0]]]) + +pl = vis.plot_stress_glyphs(mesh, stress, num_seeds=24, + save_png=True, dir_fname="stress_glyphs.png") +``` + +Seeds default to a regular grid over the mesh bounding box, filtered +to points inside the mesh (an annulus seeds no glyphs in its hole). +Pass `seeds` explicitly to sample section planes in 3-D or to avoid +regions — in fault models, keep seeds out of the weak zones, where the +recovered stress mixes materials across the interface. + +```{figure} figures/stress_glyphs_thrust.png +:alt: Two-panel figure for a thrust-ramp model. Top panel, principal stress crosses on a 26 by 13 seed grid over a 2 by 1 section with two grey parabolic ramps rising from a basal decollement, drawn over a faint peach strain-rate wash in which the quiet wedges riding each ramp glow white and the blind ramp tips smudge darker. Far-field crosses are blue-horizontal (compression); in the wedges the crosses turn red and near-vertical (tension), largest just above the blind tips. Bottom panel, stress trajectories over the same wash: dark sigma-1 lines run horizontally and arch smoothly over each ramp tip, pale sigma-3 lines rise near-vertically between them, crossing at right angles everywhere. +:name: fig-stress-glyphs-thrust + +Principal-stress crosses (top) and stress trajectories (bottom) for a +blind-thrust model, over a faint strain-rate wash (log scale, +background-percentile limits, low opacity — the second invariant of +the strain rate from the same recovered stress). The colour +convention matches the RdBu_r field convention used across the +documentation: blue compressive, red tensile. The wash ties the +stress geometry to the deformation it drives: the near-rigid thrust +wedges glow white while the blind tips concentrate both quantities. +``` + +### The pressure gauge matters for colours, not directions + +For incompressible models the full stress is +$\sigma = \tau - p\,I$ and the pressure datum is a gauge choice. +Shifting that datum shifts every principal value equally, so it can +flip bars between red and blue — but it cannot rotate the principal +directions or reorder the principal values. State the gauge in your +caption (the examples here demean the pressure), or add the +lithostatic reference before plotting if absolute compression +matters. + +The same argument answers a common question about map-view regime +colouring in the style of the World Stress Map (red normal, green +strike-slip, blue thrust): the regime classification is +gauge-invariant, but in a 2-D incompressible plane-strain model it is +also *degenerate* — the out-of-plane deviatoric stress is zero while +the in-plane deviatoric principals are $\pm s$, so the out-of-plane +stress is always the intermediate principal stress and every map-view +point classifies as strike-slip. Regime colouring only carries +information in 3-D models. + +## Stress trajectories (2-D) + +Trajectories integrate the principal *direction* field into curves — +the classical stress-trajectory diagrams of structural geology. A +principal direction is defined only modulo 180°, so ordinary +streamline tools cannot draw this field: the integrator in +`direction_trajectories` sign-aligns each evaluated eigenvector with +the previous heading, and places lines evenly (Jobard–Lehmann +occupancy) so the figure stays legible. The two families cross at +right angles wherever both are drawn — a built-in correctness check. + +The direction callable is yours to build, which keeps the integrator +independent of how the stress is stored. From nodal arrays: + +```python +import numpy as np + +def make_direction(interpolate_stress, family="compressive"): + # interpolate_stress(p) -> (sxx, syy, sxy) at point p, or None + def direction_at(p): + values = interpolate_stress(p) + if values is None: + return None + sxx, syy, sxy = values + mean_dev = 0.5 * (sxx - syy) + radius = np.hypot(mean_dev, sxy) + if radius < 1.0e-12: # isotropic point: direction undefined + return None + angle = 0.5 * np.arctan2(sxy, mean_dev) # most-tensile axis + if family == "compressive": + angle += 0.5 * np.pi + return np.array([np.cos(angle), np.sin(angle)]) + return direction_at + +lines = vis.direction_trajectories( + direction_at, candidate_seeds, inside, + step=0.008, separation=0.04, +) +trajectory_lines = vis.trajectories_to_pv_lines(lines) +``` + +Draw the compressive family dark and the tensile family pale, as in +the figure above. In 3-D the analogue of a trajectory is a surface; +we do not attempt those — draw glyphs on section planes instead. + +## Two more regimes, same recipe + +The same two panels for the extensional and strike-slip companions of +the thrust model — nested listric normal faults above a detachment, +and en-echelon segments in dextral simple shear. Nothing changes in +the code except the checkpoint being loaded. + +```{figure} figures/stress_glyphs_listric.png +:alt: Two-panel figure for three nested listric normal faults soling into a basal detachment, over a faint peach strain-rate wash with white quiet triangles in the footwall beneath each sole. Top panel, principal stress crosses: red-horizontal tension dominates the upper plate and deepens between the fault traces; blue compression concentrates near the detachment, with strongly rotated mixed crosses hugging each curved fault. Below the detachment the crosses return to uniform red-horizontal. Bottom panel, trajectories: dark sigma-1 lines hang near-vertically in the extending upper plate, bending to meet each listric sole, while pale sigma-3 lines run horizontally; beneath the detachment the net is an undisturbed rectangular grid. +:name: fig-stress-glyphs-listric + +Nested listric normal faults (extension): horizontal tension aloft, +compression concentrating under the soles, and an undisturbed +trajectory grid beneath the detachment — the faults decouple the two +plates. +``` + +```{figure} figures/stress_glyphs_enechelon.png +:alt: Two-panel map-view figure for three en-echelon fault segments in dextral simple shear, over a faint strain-rate wash with bright white lobes at the relay steps between segment tips. Top panel, principal stress crosses at 45 degrees far from the faults, shrinking and rotating through the relay zones, with a blue compressive bridge linking overlapping tips. Bottom panel, trajectories: the dark sigma-1 family sweeps diagonally across the box and kinks sharply as it hands across each relay step; the pale sigma-3 family crosses it orthogonally. +:name: fig-stress-glyphs-enechelon + +En-echelon segments in dextral simple shear (map view): the conjugate +cross pattern of the far field collapses through the relay steps, and +the trajectory families kink as stress hands across from segment to +segment. +``` + +## 3-D glyphs + +`principal_stress_glyphs` accepts `(n, 3, 3)` tensors and returns +three bars per seed. Seed one or two planes through the feature of +interest: + +```python +u = np.linspace(0.05, 0.95, 13) +gx, gz = np.meshgrid(u, u) +plane = np.column_stack([gx.ravel(), + np.full(gx.size, 0.5), # y = centre plane + gz.ravel()]) +pl = vis.plot_stress_glyphs(mesh, stress, seeds=plane) +``` + +```{figure} figures/stress_glyphs_sinker3d.png +:alt: A unit cube drawn in outline with a grey sphere of radius 0.16 near the centre, slightly above mid-height. Three-bar principal stress glyphs are drawn on a vertical section plane and a horizontal section plane through the sphere, via the one-call plot_stress_glyphs with the cube edges and sphere added to the returned plotter. Above the sphere the bars are red and near-vertical (tension as material is pulled down behind the sinker); below and beside it they are blue (compression), fanning outward on the horizontal plane beneath the sphere. Bar length decays with distance from the sphere. +:name: fig-stress-glyphs-sinker + +Three-bar principal-stress glyphs on two section planes through a +Stokes sinker: a tensile (red) column above the sinking sphere, a +compressive (blue) fan below and around it. Rendered with +`plot_stress_glyphs(..., show=False)` so the cube outline and sphere +could be added to the returned plotter before the screenshot. +``` + +## Building blocks + +The plot function is a convenience wrapper; every step is available +separately for custom figures: + +| Function | Purpose | +|---|---| +| `tensor_fn_to_pv_points(pv_mesh, uw_fn)` | Evaluate a `dim`×`dim` sympy tensor at points | +| `principal_stress_glyphs(coords, stress, scale)` | Bar segments + `"tensile"` cell array | +| `direction_trajectories(direction_at, seeds, inside, step, separation)` | Evenly spaced mod-180° trajectories | +| `trajectories_to_pv_lines(lines)` | Bundle polylines for `add_mesh` | +| `plot_stress_glyphs(mesh, stress, ...)` | One-call cross plot | + +The figures on this page come from checkpointed models (the three +fault-interaction examples and a Stokes sinker): the solve writes the +mesh, the velocity, the pressure, and the recovered stress components +with `mesh.write_timestep`, and the glyph plots load them back with +`read_timestep` — no re-solving to restyle a figure. The faint +background in the 2-D panels is the strain-rate second invariant from +the same recovered stress ($\dot\varepsilon = \tau/2$ at matrix +viscosity 1), drawn first at low opacity with log scaling and +background-percentile colour limits, so the glyphs carry the figure +and the wash only whispers where deformation concentrates. diff --git a/src/underworld3/visualisation/__init__.py b/src/underworld3/visualisation/__init__.py index 755b96cf..6b586730 100644 --- a/src/underworld3/visualisation/__init__.py +++ b/src/underworld3/visualisation/__init__.py @@ -24,5 +24,14 @@ swarm_to_pv_cloud, ) +# Principal-stress glyphs and stress trajectories (glyphs.py) +from .glyphs import ( + tensor_fn_to_pv_points, + principal_stress_glyphs, + direction_trajectories, + trajectories_to_pv_lines, + plot_stress_glyphs, +) + # Import parallel visualization utilities from . import parallel diff --git a/src/underworld3/visualisation/glyphs.py b/src/underworld3/visualisation/glyphs.py new file mode 100644 index 00000000..5e82292b --- /dev/null +++ b/src/underworld3/visualisation/glyphs.py @@ -0,0 +1,410 @@ +"""Principal-stress glyphs and stress trajectories. + +The stress tensor is sampled at seed points, the way velocity arrows +sample the velocity — not drawn everywhere. Each seed carries one bar +per principal axis: bar length proportional to the principal-value +magnitude, colour by sign (blue compressive, red tensile, matching the +RdBu_r field convention). In 2-D that is a cross; in 3-D, three +orthogonal bars. Stress trajectories integrate the principal +*direction* field into evenly spaced curves (2-D only; the 3-D +analogue is a trajectory surface, which we do not attempt — draw +glyphs on section planes instead). + +Sign conventions +---------------- +Principal values are of the stress tensor as supplied. For the full +stress :math:`\\sigma = \\tau - p I` in an incompressible model, the +pressure datum is a gauge: the compressive/tensile split (and so the +bar colours) is relative to that datum, while principal *directions* +and the ordering of principal values are gauge-invariant. State the +gauge in the caption (e.g. "demeaned pressure"), or add the +lithostatic reference before plotting if absolute signs matter. +""" + + +def tensor_fn_to_pv_points(pv_mesh, uw_fn): + """Evaluate an Underworld tensor function at PyVista mesh points. + + Parameters + ---------- + pv_mesh : pyvista.DataSet + PyVista mesh or point cloud to evaluate at. + uw_fn : sympy.Matrix + Square (``dim x dim``) Underworld tensor function, e.g. a + stress built from recovered component variables. + + Returns + ------- + numpy.ndarray + Tensor values at mesh points, shape ``(n_points, dim, dim)``. + Units are stripped; the units string is stored as + ``pv_mesh._last_tensor_units``. + """ + import numpy as np + import underworld3 as uw + + dim = uw_fn.shape[0] + if uw_fn.shape != (dim, dim): + raise ValueError(f"Expected a square tensor, got shape {uw_fn.shape}") + + if hasattr(pv_mesh, "_coord_array"): + coords = pv_mesh._coord_array[:, 0:dim] + else: + coords = pv_mesh.points[:, 0:dim] + + tensor_values = uw.function.evaluate(uw_fn, coords) + + tensor_units = None + if hasattr(tensor_values, "units") and tensor_values.units is not None: + tensor_units = str(tensor_values.units) + pv_mesh._last_tensor_units = tensor_units + + if hasattr(tensor_values, "magnitude"): + tensor_values = tensor_values.magnitude + + return np.asarray(tensor_values).reshape(-1, dim, dim) + + +def principal_stress_glyphs(coords, stress, scale): + """Principal-axis bar glyphs for symmetric tensors at seed points. + + Each seed point contributes ``dim`` line segments, one per + principal axis, centred on the seed. The half-length of the bar + for principal value :math:`\\lambda_k` along unit eigenvector + :math:`\\hat{e}_k` is :math:`\\mathrm{scale} \\cdot |\\lambda_k|`. + + Parameters + ---------- + coords : numpy.ndarray + Seed coordinates, shape ``(n, 2)`` or ``(n, 3)``. + stress : numpy.ndarray + Symmetric tensors at the seeds, shape ``(n, dim, dim)``. + The input is symmetrised (``eigh`` reads only one triangle, + so a small asymmetry from recovered components would + otherwise be ignored silently). + scale : float + Bar half-length per unit principal value, in mesh + coordinates. A good default is + ``0.45 * seed_spacing / |lambda|_max``. + + Returns + ------- + pyvista.PolyData + One line cell per bar, with cell array ``"tensile"`` + (1.0 where :math:`\\lambda_k \\ge 0`, else 0.0). Split with + ``glyphs.threshold(0.5, scalars="tensile")`` to colour the + two signs separately. + """ + import numpy as np + import pyvista as pv + + coords = np.asarray(coords, dtype=float) + n, dim = coords.shape + stress = 0.5 * (stress + np.transpose(stress, (0, 2, 1))) + + lam, vec = np.linalg.eigh(stress) + + segments = [] + tensile = [] + for k in range(dim): + half = (scale * np.abs(lam[:, k]))[:, None] * vec[:, :, k] + segments.append(np.stack([coords - half, coords + half], axis=1)) + tensile.append(lam[:, k] >= 0) + segments = np.vstack(segments).reshape(-1, dim) + tensile = np.concatenate(tensile).astype(float) + + if dim == 2: + segments = np.column_stack([segments, np.zeros(len(segments))]) + + n_bars = n * dim + lines = np.column_stack( + [ + np.full(n_bars, 2), + np.arange(0, 2 * n_bars, 2), + np.arange(1, 2 * n_bars, 2), + ] + ).ravel() + glyphs = pv.PolyData(segments, lines=lines) + glyphs.cell_data["tensile"] = tensile + return glyphs + + +def direction_trajectories( + direction_at, seeds, inside, step, separation, max_steps=2000 +): + """Evenly spaced trajectories of a direction field (2-D). + + A principal-stress direction is defined only modulo 180 degrees, + so ordinary streamline tools cannot integrate it: the integrator + here carries orientation continuity, sign-aligning each evaluated + eigenvector with the previous heading before the RK2 midpoint + step. Line placement follows Jobard & Lehmann: an occupancy grid + at ``separation`` keeps a new seed at least one spacing from + existing lines, while a running line stops only when it enters a + cell another line has actually traversed — two different tests, + because using the wide corridor for both chops lines into stubs. + + Parameters + ---------- + direction_at : callable + ``direction_at(p)`` with ``p`` shape ``(2,)`` returning a unit + direction vector, or ``None`` at isotropic points and outside + the field. The returned sign is arbitrary (mod-180 field). + seeds : numpy.ndarray + Candidate seed points, shape ``(m, 2)``. Offer more than you + expect to be used; the occupancy grid thins them. + inside : callable + ``inside(p) -> bool`` gating the domain (bounding box, keep + clear of fault zones, ...). + step : float + RK2 step length in mesh coordinates. + separation : float + Target spacing between neighbouring trajectories. + max_steps : int, optional + Cap on integration steps per direction from a seed. + + Returns + ------- + list of numpy.ndarray + Polylines, each of shape ``(n_points, 2)``. Lines shorter + than about four separations are dropped. Render with + :func:`trajectories_to_pv_lines`. + """ + import numpy as np + + occupied = set() # cells a line traversed: stops a converging line + corridor = set() # widened by one ring: blocks new seeds only + + def cell(p): + return ( + int(np.floor(p[0] / separation)), + int(np.floor(p[1] / separation)), + ) + + lines = [] + for seed in seeds: + if cell(seed) in corridor: + continue + u0 = direction_at(np.asarray(seed, dtype=float)) + if u0 is None: + continue + branches = [] + for sense in (+1.0, -1.0): + p = np.array(seed, dtype=float) + u_prev = sense * u0 + path = [p.copy()] + for _ in range(max_steps): + u1 = direction_at(p) + if u1 is None: + break + if np.dot(u1, u_prev) < 0: + u1 = -u1 + p_mid = p + 0.5 * step * u1 + if not inside(p_mid): + break + u2 = direction_at(p_mid) + if u2 is None: + break + if np.dot(u2, u1) < 0: + u2 = -u2 + p = p + step * u2 + if not inside(p) or cell(p) in occupied: + break + u_prev = u2 + path.append(p.copy()) + branches.append(np.array(path)) + line = np.vstack([branches[0][::-1], branches[1][1:]]) + if (len(line) - 1) * step >= 4.0 * separation: + lines.append(line) + # A line claims its cells only AFTER integrating, so it + # never blocks itself; the ring keeps future seeds away. + for q in line: + ci, cj = cell(q) + occupied.add((ci, cj)) + for di in (-1, 0, 1): + for dj in (-1, 0, 1): + corridor.add((ci + di, cj + dj)) + return lines + + +def trajectories_to_pv_lines(lines): + """Bundle polylines from :func:`direction_trajectories` into PolyData. + + Parameters + ---------- + lines : list of numpy.ndarray + Polylines of shape ``(n_points, 2)`` or ``(n_points, 3)``. + + Returns + ------- + pyvista.PolyData or None + One polyline cell per input line (``None`` for an empty list). + """ + import numpy as np + import pyvista as pv + + points, cells, offset = [], [], 0 + for line in lines: + if line.shape[1] == 2: + line = np.column_stack([line, np.zeros(len(line))]) + points.append(line) + cells.extend([len(line), *range(offset, offset + len(line))]) + offset += len(line) + if not points: + return None + return pv.PolyData(np.vstack(points), lines=np.asarray(cells)) + + +def plot_stress_glyphs( + mesh, + stress, + seeds=None, + num_seeds=24, + scale=None, + compressive_colour="#2166ac", + tensile_colour="#b2182b", + line_width=2.5, + show_edges=False, + save_png=False, + dir_fname="", + title="", + cpos="xy", + window_size=(750, 750), + show=True, +): + """Plot principal-stress crosses sampled at seed points. + + Parameters + ---------- + mesh : uw.discretisation.Mesh + The mesh the stress lives on. + stress : sympy.Matrix + Square (``dim x dim``) stress function. Use recovered + (projected) component variables rather than raw solver + derivative expressions — see the visualisation guide. + seeds : numpy.ndarray, optional + Seed coordinates, shape ``(n, dim)``. Default is a regular + grid over the mesh bounding box, filtered to points inside + the mesh (so an annulus hole seeds nothing). In 3-D, prefer + passing seeds on one or two section planes; a filled volume + grid of three-bar glyphs is hard to read. + num_seeds : int, optional + Grid resolution across the longest box edge for the default + seeding. + scale : float, optional + Bar half-length per unit principal value. Default scales the + largest bar to 0.45 of the seed spacing. + compressive_colour, tensile_colour : str, optional + Bar colours for negative / non-negative principal values. + line_width : float, optional + Bar line width in pixels. + show_edges : bool, optional + Draw the mesh edge skeleton faintly beneath the glyphs. + save_png : bool, optional + Save a screenshot to ``dir_fname``. + dir_fname : str, optional + Path for the screenshot when ``save_png`` is True. + title : str, optional + Text placed on the figure. + cpos : str or list, optional + PyVista camera position (``"xy"`` for 2-D sections). + window_size : tuple of int, optional + Render window size in pixels. + show : bool, optional + Call ``show()`` before returning. Pass ``False`` to add + overlays to the returned plotter and screenshot it yourself — + actors added after ``show()`` are ignored by a finalized + scene. + + Returns + ------- + pyvista.Plotter + The plotter. With ``show=False`` it is still open: add + overlays, then ``show()`` or ``screenshot()``. + """ + import numpy as np + import pyvista as pv + + from .visualisation import mesh_to_pv_mesh + + dim = mesh.dim + pvmesh = mesh_to_pv_mesh(mesh) + + if seeds is None: + bounds = np.asarray(pvmesh.bounds).reshape(3, 2) + extents = bounds[:, 1] - bounds[:, 0] + spacing = extents[:dim].max() / num_seeds + axes = [ + np.arange(bounds[k, 0] + 0.5 * spacing, bounds[k, 1], spacing) + for k in range(dim) + ] + grids = np.meshgrid(*axes, indexing="ij") + seeds = np.column_stack([g.ravel() for g in grids]) + # Keep only seeds inside the mesh: an interior point is its own + # closest point in the containing cell, an exterior point is not. + probe = seeds + if dim == 2: + probe = np.column_stack([seeds, np.zeros(len(seeds))]) + _cells, closest = pvmesh.find_closest_cell( + probe, return_closest_point=True + ) + distance = np.linalg.norm(closest - probe, axis=1) + seeds = seeds[distance < 1.0e-6 * extents.max()] + else: + seeds = np.asarray(seeds, dtype=float) + # The auto-scale needs the ACTUAL seed spacing — user seeds + # owe nothing to num_seeds. Mean nearest-neighbour distance, + # subsampled so the pairwise matrix stays small. + sample = seeds + if len(sample) > 2048: + sample = sample[:: len(sample) // 2048 + 1] + offsets = sample[:, None, :] - sample[None, :, :] + distance2 = np.sum(offsets**2, axis=2) + np.fill_diagonal(distance2, np.inf) + spacing = float(np.sqrt(distance2.min(axis=1)).mean()) + + cloud = pv.PolyData( + np.column_stack([seeds, np.zeros(len(seeds))]) + if dim == 2 + else seeds + ) + stress_values = tensor_fn_to_pv_points(cloud, stress) + + if scale is None: + lam = np.linalg.eigvalsh( + 0.5 * (stress_values + np.transpose(stress_values, (0, 2, 1))) + ) + scale = 0.45 * spacing / np.abs(lam).max() + + glyphs = principal_stress_glyphs(seeds, stress_values, scale) + + pl = pv.Plotter(window_size=window_size) + pl.set_background("white") + if show_edges: + pl.add_mesh( + pvmesh.extract_all_edges(), + color="#d9d9d9", + line_width=0.5, + lighting=False, + ) + for threshold_value, invert, colour in ( + (0.5, True, compressive_colour), + (0.5, False, tensile_colour), + ): + part = glyphs.threshold( + threshold_value, scalars="tensile", invert=invert + ) + if part.n_cells: + pl.add_mesh( + part, color=colour, line_width=line_width, lighting=False + ) + if len(title): + pl.add_text(title, font_size=11, color="black") + + pl.camera_position = cpos + if show: + pl.show() + if save_png: + pl.screenshot(dir_fname) + + return pl diff --git a/tests/test_0848_stress_glyphs.py b/tests/test_0848_stress_glyphs.py new file mode 100644 index 00000000..29739528 --- /dev/null +++ b/tests/test_0848_stress_glyphs.py @@ -0,0 +1,145 @@ +"""Principal-stress glyph and trajectory geometry. + +These are figure primitives, so what can be asserted is the geometry +they BUILD, not what it looks like: that a known tensor produces bars +of the right length, direction, and sign class, and that the mod-180 +trajectory integrator follows a direction field whose eigenvector +sign flips underfoot — the case an ordinary streamline integrator +gets wrong by reversing mid-line. +""" +import numpy as np +import pytest + +import underworld3.visualisation as vis + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + +import pyvista + +pyvista.OFF_SCREEN = True + + +@pytest.fixture(autouse=True) +def _close_every_plotter(): + yield + pyvista.close_all() + + +def test_uniaxial_compression_cross(): + # sigma = diag(-2, 0): one compressive bar of half-length 2*scale + # along x, one zero-length tensile bar along y. + coords = np.array([[0.0, 0.0]]) + stress = np.array([[[-2.0, 0.0], [0.0, 0.0]]]) + glyphs = vis.principal_stress_glyphs(coords, stress, scale=0.25) + + assert glyphs.n_cells == 2 + tensile = glyphs.cell_data["tensile"] + assert sorted(tensile) == [0.0, 1.0] + + # Bar k's endpoints are consecutive points 2k, 2k+1. + spans = np.ptp(np.asarray(glyphs.points).reshape(-1, 2, 3), axis=1) + lengths = np.linalg.norm(spans, axis=1) + compressive = int(np.argmin(tensile)) + assert np.isclose(lengths[compressive], 2 * 0.25 * 2.0) + assert np.isclose(lengths[int(np.argmax(tensile))], 0.0) + assert np.allclose(spans[compressive][1:], 0.0) + + +def test_shear_gives_45_degree_cross(): + # Pure shear sigma_xy = 1: principal axes at 45 degrees, one + # tensile and one compressive bar of equal length. + coords = np.array([[0.0, 0.0]]) + stress = np.array([[[0.0, 1.0], [1.0, 0.0]]]) + glyphs = vis.principal_stress_glyphs(coords, stress, scale=1.0) + + spans = np.ptp(np.asarray(glyphs.points).reshape(-1, 2, 3), axis=1) + assert np.allclose(spans[:, 0], spans[:, 1]) + assert sorted(glyphs.cell_data["tensile"]) == [0.0, 1.0] + + +def test_three_bars_in_3d(): + coords = np.array([[0.0, 0.0, 0.0]]) + stress = np.array([np.diag([-3.0, 1.0, 2.0])]) + glyphs = vis.principal_stress_glyphs(coords, stress, scale=1.0) + + assert glyphs.n_cells == 3 + assert sorted(glyphs.cell_data["tensile"]) == [0.0, 1.0, 1.0] + + +def test_trajectory_survives_eigenvector_sign_flip(): + # A uniform horizontal direction field whose reported sign + # alternates with x — legitimate for eigenvectors (mod 180). The + # integrator must keep heading one way and cross the whole box. + def direction_at(p): + sign = 1.0 if np.sin(20.0 * p[0]) >= 0 else -1.0 + return sign * np.array([1.0, 0.0]) + + def inside(p): + return 0.0 <= p[0] <= 1.0 and 0.0 <= p[1] <= 1.0 + + seeds = np.array([[0.5, 0.5]]) + lines = vis.direction_trajectories( + direction_at, seeds, inside, step=0.01, separation=0.05 + ) + assert len(lines) == 1 + line = lines[0] + assert np.ptp(line[:, 0]) > 0.9 # spans the box + assert np.ptp(line[:, 1]) < 1.0e-12 # never turns + + +def test_annulus_default_seeds_avoid_the_hole(): + # The default seed grid spans the bounding box; on an annulus the + # box centre is not in the mesh, and evaluating there would fail. + import sympy + import underworld3 as uw + + mesh = uw.meshing.Annulus(radiusInner=0.5, radiusOuter=1.0, + cellSize=0.2) + x, y = mesh.X + stress = sympy.Matrix([[x, y], [y, -x]]) + + # show=False: the test asserts what the plot BUILDS, and a CI + # runner must never enter a render/interactor path — a worker + # that calls show() can hang the whole xdist session (run + # 32097850915 sat from 90% to the 2 h timeout). + pl = vis.plot_stress_glyphs(mesh, stress, num_seeds=12, show=False) + pl.close() + + # Rebuild the same default seeding to inspect it directly. + import pyvista as pv + + pvmesh = vis.mesh_to_pv_mesh(mesh) + bounds = np.asarray(pvmesh.bounds).reshape(3, 2) + spacing = (bounds[:2, 1] - bounds[:2, 0]).max() / 12 + axes = [ + np.arange(bounds[k, 0] + 0.5 * spacing, bounds[k, 1], spacing) + for k in range(2) + ] + gx, gy = np.meshgrid(*axes, indexing="ij") + seeds = np.column_stack([gx.ravel(), gy.ravel(), np.zeros(gx.size)]) + _cells, closest = pvmesh.find_closest_cell( + seeds, return_closest_point=True + ) + inside = np.linalg.norm(closest - seeds, axis=1) < 1.0e-6 + radii = np.linalg.norm(seeds[inside, :2], axis=1) + assert inside.sum() > 0 + assert radii.min() > 0.45 # nothing seeded in the hole + assert radii.max() < 1.0 + 1.0e-6 + + +def test_trajectories_respect_separation(): + def direction_at(p): + return np.array([1.0, 0.0]) + + def inside(p): + return 0.0 <= p[0] <= 1.0 and 0.0 <= p[1] <= 1.0 + + gx, gy = np.meshgrid(np.linspace(0.1, 0.9, 9), np.linspace(0.1, 0.9, 9)) + seeds = np.column_stack([gx.ravel(), gy.ravel()]) + lines = vis.direction_trajectories( + direction_at, seeds, inside, step=0.01, separation=0.1 + ) + assert len(lines) > 1 + heights = sorted(line[0, 1] for line in lines) + gaps = np.diff(heights) + assert gaps.min() > 0.05 # no two lines share a corridor