{ "cells": [ { "cell_type": "markdown", "id": "r3d-title", "metadata": {}, "source": [ "# Rare-event identification on a 3D field\n", "\n", "The 3D analogue of :doc:`rei-example-2d`: build a small 3D von-Mises stress\n", "field with a few high-stress blobs, cluster it, and view a couple of z-slices of\n", "the field next to the recovered clusters. Pure Python; runs from `pip install\n", "graintrace`.\n", "\n", "See :doc:`/algorithms/rare-event-identification` and\n", ":class:`~graintrace.ClusterAnalysisIndicator`." ] }, { "cell_type": "code", "execution_count": null, "id": "r3d-imports", "metadata": { "execution": { "iopub.execute_input": "2026-10-09T23:30:10.642730Z", "iopub.status.busy": "2026-10-09T23:30:10.642242Z", "iopub.status.idle": "2026-10-09T23:30:11.211714Z", "shell.execute_reply": "2026-10-09T23:30:11.210309Z" } }, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from matplotlib.colors import ListedColormap, BoundaryNorm\n", "\n", "from graintrace.cluster_indicator import ClusterAnalysisIndicator\n", "from graintrace.similarity_metric_library import SimilarityMetricLibrary" ] }, { "cell_type": "markdown", "id": "r3d-gen-md", "metadata": {}, "source": [ "## Build a synthetic 3D field\n", "\n", "A low-stress background with a few high-stress spheres. Each voxel becomes a row\n", "with a Cauchy-stress tensor (`sxx..syz`)." ] }, { "cell_type": "code", "execution_count": null, "id": "r3d-gen", "metadata": { "execution": { "iopub.execute_input": "2026-10-09T23:30:11.216759Z", "iopub.status.busy": "2026-10-09T23:30:11.216531Z", "iopub.status.idle": "2026-10-09T23:30:11.318258Z", "shell.execute_reply": "2026-10-09T23:30:11.316730Z" } }, "outputs": [], "source": [ "def make_vms_volume(path, n=16, n_blobs=6, vm_low=50.0, vm_high=200.0,\n", " radius_range=(2, 4), seed=42):\n", " rng = np.random.default_rng(seed)\n", " vm = np.full((n, n, n), vm_low)\n", " zz, yy, xx = np.ogrid[:n, :n, :n]\n", " for _ in range(n_blobs):\n", " cx, cy, cz = rng.integers(0, n, size=3)\n", " r = rng.integers(*radius_range)\n", " vm[(xx - cx) ** 2 + (yy - cy) ** 2 + (zz - cz) ** 2 <= r * r] = vm_high\n", " rows = []\n", " eid = 0\n", " for k in range(n):\n", " for j in range(n):\n", " for i in range(n):\n", " eid += 1\n", " t = vm[k, j, i]\n", " rows.append(dict(id=eid, x=float(i), y=float(j), z=float(k),\n", " sxx=t + rng.normal(0, t * 0.05),\n", " syy=rng.normal(0, t * 0.02), szz=rng.normal(0, t * 0.02),\n", " sxy=rng.normal(0, t * 0.02), sxz=rng.normal(0, t * 0.02),\n", " syz=rng.normal(0, t * 0.02)))\n", " pd.DataFrame(rows).to_csv(path, index=False)\n", " return n\n", "\n", "csv_path = \"rei3d_field.csv\"\n", "n = make_vms_volume(csv_path)\n", "print(f\"wrote {csv_path}: {n}x{n}x{n} volume\")" ] }, { "cell_type": "markdown", "id": "r3d-cluster-md", "metadata": {}, "source": [ "## Cluster the volume" ] }, { "cell_type": "code", "execution_count": null, "id": "r3d-cluster", "metadata": { "execution": { "iopub.execute_input": "2026-10-09T23:30:11.321607Z", "iopub.status.busy": "2026-10-09T23:30:11.321446Z", "iopub.status.idle": "2026-10-09T23:32:40.259368Z", "shell.execute_reply": "2026-10-09T23:32:40.256948Z" } }, "outputs": [], "source": [ "indicator = ClusterAnalysisIndicator(csv_path, coord_cols=(\"x\", \"y\", \"z\"))\n", "spec = SimilarityMetricLibrary().von_mises_stress()\n", "run = indicator.run(method_type=\"scipy_hierarchical\", spec=spec,\n", " threshold=0.01, method=\"average\", criterion=\"distance\")\n", "result = run[\"points\"]\n", "print(\"n clusters:\", result[\"cluster_label\"].nunique())" ] }, { "cell_type": "markdown", "id": "r3d-viz-md", "metadata": {}, "source": [ "## View z-slices" ] }, { "cell_type": "code", "execution_count": null, "id": "r3d-viz", "metadata": { "execution": { "iopub.execute_input": "2026-10-09T23:32:40.265404Z", "iopub.status.busy": "2026-10-09T23:32:40.264888Z", "iopub.status.idle": "2026-10-09T23:32:40.879486Z", "shell.execute_reply": "2026-10-09T23:32:40.877779Z" } }, "outputs": [], "source": [ "def von_mises(df):\n", " s = {k: df[k].to_numpy() for k in (\"sxx\", \"syy\", \"szz\", \"sxy\", \"sxz\", \"syz\")}\n", " t1 = ((s[\"sxx\"] - s[\"syy\"]) ** 2 + (s[\"syy\"] - s[\"szz\"]) ** 2\n", " + (s[\"szz\"] - s[\"sxx\"]) ** 2) / 2.0\n", " return np.sqrt(t1 + 3.0 * (s[\"sxy\"] ** 2 + s[\"sxz\"] ** 2 + s[\"syz\"] ** 2))\n", "\n", "vm = von_mises(result).reshape(n, n, n)\n", "labels = result[\"cluster_label\"].to_numpy()\n", "uniq = np.unique(labels)\n", "idx = {lab: i for i, lab in enumerate(uniq)}\n", "lab = np.vectorize(idx.get)(labels).reshape(n, n, n)\n", "cmap = ListedColormap(plt.get_cmap(\"tab20\")(np.linspace(0, 1, len(uniq))))\n", "norm = BoundaryNorm(np.arange(len(uniq) + 1) - 0.5, len(uniq))\n", "\n", "slices = [n // 4, n // 2, 3 * n // 4]\n", "fig, ax = plt.subplots(2, len(slices), figsize=(9, 6))\n", "for c, kz in enumerate(slices):\n", " ax[0, c].imshow(vm[kz], origin=\"lower\", cmap=\"viridis\")\n", " ax[0, c].set_title(f\"von-Mises, z={kz}\")\n", " ax[1, c].imshow(lab[kz], origin=\"lower\", cmap=cmap, norm=norm)\n", " ax[1, c].set_title(f\"clusters, z={kz}\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "r3d-seealso", "metadata": {}, "source": [ "**See also**\n", "\n", "- Algorithm: :doc:`/algorithms/rare-event-identification`\n", "- Full pipeline with graph clustering and VTK export: :doc:`rare-event-identification`\n", "- 2D version: :doc:`rei-example-2d`" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }