This tutorial is pure Python and runs from pip install graintrace — open it in Colab:

Rare-event identification on a 3D field

The 3D analogue of :doc:rei-example-2d: build a small 3D von-Mises stress field with a few high-stress blobs, cluster it, and view a couple of z-slices of the field next to the recovered clusters. Pure Python; runs from pip install graintrace.

See :doc:/algorithms/rare-event-identification and :class:~graintrace.ClusterAnalysisIndicator.

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, BoundaryNorm

from graintrace.cluster_indicator import ClusterAnalysisIndicator
from graintrace.similarity_metric_library import SimilarityMetricLibrary

Build a synthetic 3D field

A low-stress background with a few high-stress spheres. Each voxel becomes a row with a Cauchy-stress tensor (sxx..syz).

def make_vms_volume(path, n=16, n_blobs=6, vm_low=50.0, vm_high=200.0,
                    radius_range=(2, 4), seed=42):
    rng = np.random.default_rng(seed)
    vm = np.full((n, n, n), vm_low)
    zz, yy, xx = np.ogrid[:n, :n, :n]
    for _ in range(n_blobs):
        cx, cy, cz = rng.integers(0, n, size=3)
        r = rng.integers(*radius_range)
        vm[(xx - cx) ** 2 + (yy - cy) ** 2 + (zz - cz) ** 2 <= r * r] = vm_high
    rows = []
    eid = 0
    for k in range(n):
        for j in range(n):
            for i in range(n):
                eid += 1
                t = vm[k, j, i]
                rows.append(dict(id=eid, x=float(i), y=float(j), z=float(k),
                                 sxx=t + rng.normal(0, t * 0.05),
                                 syy=rng.normal(0, t * 0.02), szz=rng.normal(0, t * 0.02),
                                 sxy=rng.normal(0, t * 0.02), sxz=rng.normal(0, t * 0.02),
                                 syz=rng.normal(0, t * 0.02)))
    pd.DataFrame(rows).to_csv(path, index=False)
    return n

csv_path = "rei3d_field.csv"
n = make_vms_volume(csv_path)
print(f"wrote {csv_path}: {n}x{n}x{n} volume")

wrote rei3d_field.csv: 16x16x16 volume

Cluster the volume

indicator = ClusterAnalysisIndicator(csv_path, coord_cols=("x", "y", "z"))
spec = SimilarityMetricLibrary().von_mises_stress()
run = indicator.run(method_type="scipy_hierarchical", spec=spec,
                    threshold=0.01, method="average", criterion="distance")
result = run["points"]
print("n clusters:", result["cluster_label"].nunique())

/opt/hostedtoolcache/Python/3.11.17/x64/lib/python3.11/site-packages/sklearn/manifold/_mds.py:735: FutureWarning: The default value of `init` will change from 'random' to 'classical_mds' in 1.10. To suppress this warning, provide some value of `init`.
  warnings.warn(
/opt/hostedtoolcache/Python/3.11.17/x64/lib/python3.11/site-packages/sklearn/manifold/_mds.py:752: FutureWarning: The `dissimilarity` parameter is deprecated and will be removed in 1.10. Use `metric` instead.
  warnings.warn(

n clusters: 24

View z-slices

def von_mises(df):
    s = {k: df[k].to_numpy() for k in ("sxx", "syy", "szz", "sxy", "sxz", "syz")}
    t1 = ((s["sxx"] - s["syy"]) ** 2 + (s["syy"] - s["szz"]) ** 2
          + (s["szz"] - s["sxx"]) ** 2) / 2.0
    return np.sqrt(t1 + 3.0 * (s["sxy"] ** 2 + s["sxz"] ** 2 + s["syz"] ** 2))

vm = von_mises(result).reshape(n, n, n)
labels = result["cluster_label"].to_numpy()
uniq = np.unique(labels)
idx = {lab: i for i, lab in enumerate(uniq)}
lab = np.vectorize(idx.get)(labels).reshape(n, n, n)
cmap = ListedColormap(plt.get_cmap("tab20")(np.linspace(0, 1, len(uniq))))
norm = BoundaryNorm(np.arange(len(uniq) + 1) - 0.5, len(uniq))

slices = [n // 4, n // 2, 3 * n // 4]
fig, ax = plt.subplots(2, len(slices), figsize=(9, 6))
for c, kz in enumerate(slices):
    ax[0, c].imshow(vm[kz], origin="lower", cmap="viridis")
    ax[0, c].set_title(f"von-Mises, z={kz}")
    ax[1, c].imshow(lab[kz], origin="lower", cmap=cmap, norm=norm)
    ax[1, c].set_title(f"clusters, z={kz}")
plt.tight_layout()
plt.show()

See also

• Algorithm: :doc:/algorithms/rare-event-identification

• Full pipeline with graph clustering and VTK export: :doc:rare-event-identification

• 2D version: :doc:rei-example-2d