diff --git a/docs/index.md b/docs/index.md index a013a8a6e..8d5081f78 100644 --- a/docs/index.md +++ b/docs/index.md @@ -143,6 +143,7 @@ tutorial/index tutorial/01_phase_field_primer tutorial/02_visual_glossary tutorial/notebook_square_plate_fracture +tutorial/notebook_layered_dcb tutorial/notebook_dynamic_sent tutorial/notebook_setup tutorial/notebook_mesh_resolution diff --git a/docs/tutorial/index.md b/docs/tutorial/index.md index 003667519..2ed0d1c5b 100644 --- a/docs/tutorial/index.md +++ b/docs/tutorial/index.md @@ -33,6 +33,10 @@ silently migrate them or add 3D fracture support. 5. [Retained Miehe SENT results](notebook_retained_results.ipynb) examines checked-in load-displacement and damage evidence and states the current post-processing boundary. +6. [Layered DCB: understand, edit, and interpret a fracture model](notebook_layered_dcb.ipynb) + explains material regions, the initial crack, boundary conditions, the + implicit staggered solution, and result interpretation. It displays retained + evidence by default and makes a new simulation an explicit choice. Asymmetric three-point bending and L-shaped panel notebooks are not presented as public benchmarks because the current repository does not retain the @@ -55,6 +59,7 @@ checkout, and prints the resolved commit before installation. | Tutorial | Time | What you learn | |---|---:|---| | [Standard simulation workflow](07_standard_simulation_workflow.md) | Self-paced; full DCB runtime is machine-dependent | Use one schema-2 layout across the small SENT exercises and full DCB, edit inputs, and distinguish numerical and physical evidence. | +| [Layered DCB student notebook](notebook_layered_dcb.ipynb) | Self-paced; full solve optional | Connect the physical model and equations to editable YAML, view a geometry schematic, inspect retained fields and animation, and design a material-contrast comparison. | | [Getting started](../getting-started.md) | 10-20 min on a fresh machine | Create an environment, install PyTorch and PhAST, run `doctor`, validate a shipped example, and inspect a result. A prepared teaching environment is faster. | | [Square-plate fracture: complete workflow](notebook_square_plate_fracture.ipynb) | 25-40 min including installation; about 15-60 s for the solve on a prepared laptop | Generate the B3 square-plate mesh, inspect named boundaries, assemble a CFL-safe teaching configuration, run explicit mechanics with an implicit damage update, and create fresh field plots and a crack animation. Runtime is machine-dependent. | | [Dynamic SENT example](notebook_dynamic_sent.ipynb) | 20-30 min | Inspect the existing B3 mesh, named regions, loading, explicit solver route, retained histories, and crack-growth animation without presenting retained evidence as a new run. | diff --git a/docs/tutorial/notebook_layered_dcb.ipynb b/docs/tutorial/notebook_layered_dcb.ipynb new file mode 100644 index 000000000..edb9821b9 --- /dev/null +++ b/docs/tutorial/notebook_layered_dcb.ipynb @@ -0,0 +1,926 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "dcb-00", + "metadata": {}, + "source": [ + "# Understanding a layered DCB simulation with PhAST\n", + "\n", + "**Physical question:** how does a crack opening along a weak bulk layer interact with a circular region that is harder to fracture?\n", + "\n", + "A double-cantilever-beam (DCB) specimen opens in two arms above and below an initial crack. This **DCB-style teaching model** uses that loading idea. It is not an ASTM test reproduction or an experimentally calibrated specimen.\n", + "\n", + "By completing this notebook, you will be able to:\n", + "\n", + "- identify geometry, mesh, material regions, the initial crack, and boundary conditions;\n", + "- connect displacement, strain, stress, and phase-field damage to their physical meanings;\n", + "- explain the implicit staggered solution used in this quasi-static example;\n", + "- edit a copy of the standard YAML and submit it through the normal PhAST commands;\n", + "- distinguish prescribed initial damage from computed crack growth;\n", + "- interpret images, an animation, response histories, and a controlled material comparison.\n", + "\n", + "**Default mode:** inspect the checked-in results without rerunning the fracture simulation. The optional full solve is disabled by `RUN_NEW_SIMULATION = False`. The retained fine calculation took about 968 seconds on its recorded CPU run; this is not a runtime estimate for your computer.\n", + "\n", + "**Scope:** two-dimensional, small-strain, plane-stress, quasi-static Amor AT2 fracture on a structured triangular mesh. The three regions are bulk materials. A weak layer is not an independent cohesive or interphase-debonding law. The retained results do not establish penetration, completed bypass, mesh independence, or physical validation.\n", + "\n", + "This notebook uses the existing PhAST configuration, CLI, result reader, and comparison script. It does not implement another FEM solver.\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-01", + "metadata": {}, + "source": [ + "## 1. Install and open the notebook\n", + "\n", + "Use **Python 3.10 or newer**. Python 3.11 is a suitable choice for a teaching environment. This multi-material route uses CPU calculations in double precision. A GPU, a separate PhAST C++ compilation, and optional external solver packages are not required for this example.\n", + "\n", + "Clone the complete repository, not only this notebook: the example configuration and retained figures are required inputs. Start Jupyter from the repository root and select the **Python (PhAST)** kernel.\n", + "\n", + "### macOS or Linux\n", + "\n", + "~~~bash\n", + "# Clone the source, including example inputs and retained results.\n", + "git clone https://github.com/CEMS-Lab/PhAST.git\n", + "cd PhAST\n", + "\n", + "# Use an installed Python >=3.10 to create an isolated environment.\n", + "python3 -m venv .venv\n", + "source .venv/bin/activate\n", + "\n", + "# Install PhAST and Jupyter into that same environment.\n", + "python -m pip install --upgrade pip\n", + "python -m pip install -e .\n", + "python -m pip install jupyterlab ipykernel\n", + "python -m ipykernel install --user --name phast --display-name \"Python (PhAST)\"\n", + "\n", + "# Open the explainer.\n", + "python -m jupyter lab docs/tutorial/notebook_layered_dcb.ipynb\n", + "~~~\n", + "\n", + "### Windows PowerShell\n", + "\n", + "~~~powershell\n", + "# Clone the same source.\n", + "git clone https://github.com/CEMS-Lab/PhAST.git\n", + "cd PhAST\n", + "\n", + "# The selected Python 3 installation must be version 3.10 or newer.\n", + "py -3 -m venv .venv\n", + "\n", + "# Calling the environment's Python directly avoids activation-policy issues.\n", + ".\\.venv\\Scripts\\python.exe -m pip install --upgrade pip\n", + ".\\.venv\\Scripts\\python.exe -m pip install -e .\n", + ".\\.venv\\Scripts\\python.exe -m pip install jupyterlab ipykernel\n", + ".\\.venv\\Scripts\\python.exe -m ipykernel install --user --name phast --display-name \"Python (PhAST)\"\n", + ".\\.venv\\Scripts\\python.exe -m jupyter lab docs/tutorial/notebook_layered_dcb.ipynb\n", + "~~~\n", + "\n", + "The base installation supplies the solver and plotting dependencies. The optional `python -m pip install -r requirements.txt` route follows the repository dependency list. Consult [Installation](../install.md) for platform-specific and optional backend requirements.\n", + "\n", + "**Already installed?** Do not clone again. Open this notebook from the checkout you installed. Run cells from top to bottom. After changing the parameter cell, restart the kernel and run all cells so that the schematic, configuration, and output paths remain consistent.\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-02", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "# Configuration, plotting, and calls to PhAST; no replacement FEM implementation.\n", + "from __future__ import annotations\n", + "\n", + "import copy\n", + "import csv\n", + "import json\n", + "import math\n", + "import os\n", + "from pathlib import Path\n", + "import platform\n", + "import subprocess\n", + "import sys\n", + "import uuid\n", + "\n", + "if sys.version_info < (3, 10):\n", + " raise RuntimeError(\"Select a notebook kernel with Python 3.10 or newer.\")\n", + "\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.patches import Circle, Rectangle\n", + "import numpy as np\n", + "import yaml\n", + "from IPython.display import Image, Markdown, display\n", + "import phast\n", + "\n", + "# Set this only if automatic discovery cannot locate your cloned repository.\n", + "REPOSITORY_PATH: str | None = None\n", + "\n", + "\n", + "def find_repository(override: str | None) -> Path:\n", + " if override is not None:\n", + " candidates = [Path(override).expanduser().resolve()]\n", + " else:\n", + " current = Path.cwd().resolve()\n", + " package = Path(phast.__file__).resolve()\n", + " candidates = [current, *current.parents, *package.parents]\n", + " for candidate in candidates:\n", + " if (candidate / \"examples/two_material_dcb_beta/config.yaml\").is_file():\n", + " return candidate\n", + " raise FileNotFoundError(\n", + " \"Clone the complete PhAST repository and, if necessary, \"\n", + " \"set REPOSITORY_PATH to that checkout.\"\n", + " )\n", + "\n", + "\n", + "ROOT = find_repository(REPOSITORY_PATH)\n", + "if not Path(phast.__file__).resolve().is_relative_to((ROOT / \"src/phast\").resolve()):\n", + " raise RuntimeError(\n", + " \"The selected kernel imports a different PhAST installation. \"\n", + " \"Install this checkout with python -m pip install -e . and restart the kernel.\"\n", + " )\n", + "\n", + "EXAMPLE = ROOT / \"examples/two_material_dcb_beta\"\n", + "RETAINED = EXAMPLE / \"results\"\n", + "REFERENCE_CONFIG = EXAMPLE / \"config.yaml\"\n", + "reference_config = yaml.safe_load(REFERENCE_CONFIG.read_text(encoding=\"utf-8\"))\n", + "config = copy.deepcopy(reference_config) # Keep the checked-in reference unchanged.\n", + "\n", + "plt.style.use(\"default\")\n", + "plt.rcParams.update({\"figure.facecolor\": \"white\", \"axes.facecolor\": \"white\"})\n", + "print(\"Python:\", platform.python_version())\n", + "print(\"PhAST:\", getattr(phast, \"__version__\", \"version not exposed\"))\n", + "print(\"Repository:\", ROOT)\n", + "print(\"Reference input:\", REFERENCE_CONFIG.relative_to(ROOT))\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-03", + "metadata": {}, + "source": [ + "## 2. State the physical problem\n", + "\n", + "![Reference input diagram showing three coloured material regions, the starter crack, the right clamp, and opposite opening displacements.](../../examples/two_material_dcb_beta/results/material_regions_and_loading.png)\n", + "\n", + "Blue identifies the outer Material 1, grey the weak bulk layer, and ochre the circular Material 2. This is a reference input diagram, not a new simulation result. The diagram and the editable schematic below use the same colours; red remains reserved for damage. The final prescribed displacements are $+0.15$ mm and $-0.15$ mm, giving a relative opening of $0.30$ mm.\n", + "\n", + "The reference rectangle is 6.0 mm long and 1.2 mm high. A weak bulk layer occupies its central 0.6 mm. A circular region of radius 0.20 mm is centred at $(2.2,0)$ mm inside that layer. An initial damaged line extends from the left edge to $x=1.5$ mm.\n", + "\n", + "The right edge is fixed. Opposite vertical displacements open the upper and lower portions of the left edge. The final values are $+0.15$ mm above and $-0.15$ mm below: **0.30 mm total opening**.\n", + "\n", + "| Material key | Meaning | Reference $E$ [MPa] | Reference $G_c$ [N/mm] |\n", + "|---|---|---:|---:|\n", + "| `material_1` | Outer fracture-resistant material | 5,000 | 1.00 |\n", + "| `weak_midline` | Finite-width weak bulk layer | 5,000 | 0.04 |\n", + "| `material_2` | Tougher circular region inside the layer | 5,000 | 0.12 |\n", + "\n", + "The elastic stiffness is initially identical in all three regions. The circle is therefore **tougher, not stiffer**, than the weak layer. This is a three-region model despite the example directory's historical name.\n", + "\n", + "**Before viewing the result:** sketch two possible crack paths. Would a larger $G_c$ necessarily force a crack around the circle? Keep your prediction separate from what the computed fields show.\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-04", + "metadata": {}, + "source": [ + "## 3. Edit the study parameters in one place\n", + "\n", + "The defaults below reproduce the declared fine-reference input values. They do not guarantee bit-identical reproduction of a historical run: its exact generating source revision was not recorded.\n", + "\n", + "Students normally edit this cell. The following cells map these variables into the existing schema-2 configuration, rather than requiring a new simulation program.\n", + "\n", + "**Units:** lengths and displacements are in mm, $E$ is in MPa ($\\mathrm{N/mm^2}$), and $G_c$ is in N/mm. Poisson's ratio $\\nu$ and damage $d$ are dimensionless. The regularisation length $\\ell_0$ is a model parameter, not the mesh spacing. Density does not enter this quasi-static equilibrium calculation.\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-05", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "# Enable only when you intend to run the full fracture calculation.\n", + "RUN_NEW_SIMULATION = False\n", + "\n", + "# Geometry and mesh. An even ny places nodes on the midline y = 0.\n", + "L = 6.0 # Length [mm].\n", + "H = 1.2 # Height [mm].\n", + "layer_height = 0.60 # Weak bulk-layer thickness [mm].\n", + "circle_x = 2.20 # Material 2 centre [mm].\n", + "circle_y = 0.0\n", + "circle_radius = 0.20 # Material 2 radius [mm].\n", + "initial_crack_length = 1.50\n", + "crack_selection_width = 0.012 # Selection distance [mm], not crack opening.\n", + "nx = 160 # Cells along x; each cell is split into two T3 triangles.\n", + "ny = 48\n", + "\n", + "# Material parameters. YAML uses Gc for G_c and l0 for ell_0.\n", + "E = 5000.0 # Outer-material and weak-layer Young's modulus [MPa].\n", + "E_material_2 = 5000.0 # Keep equal to E to isolate the toughness contrast.\n", + "nu = 0.30\n", + "G_c_outer = 1.00 # Outer material fracture resistance [N/mm].\n", + "G_c_layer = 0.04 # Weak-layer fracture resistance [N/mm].\n", + "G_c_material_2 = 0.12 # Circular-region fracture resistance [N/mm].\n", + "ell_0 = 0.10 # Shared regularisation length [mm].\n", + "eta = 1.0e-6 # Shared residual stiffness parameter, dimensionless.\n", + "\n", + "# This is a load sequence, not a duration in seconds.\n", + "final_opening = 0.30 # Total separation of the loaded arms [mm].\n", + "number_of_load_levels = 121 # Includes zero: 120 nominal nonzero increments.\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-06", + "metadata": {}, + "source": [ + "## 4. Geometry, mesh, and named regions\n", + "\n", + "**Geometry** specifies the physical dimensions. **Meshing** divides the domain into finite elements. A T3 element is a three-node triangle with linear interpolation of the nodal fields.\n", + "\n", + "PhAST constructs this mesh using `geometry.type: structured_rectangle`. This example does not require a separate geometry file. Named regions select material elements, the starter crack, and boundary nodes.\n", + "\n", + "The next cell keeps the circle and its complementary selections consistent. It updates the layer rectangles and load bounds when dimensions change. The right clamp uses `side: right` and therefore follows the right edge.\n", + "\n", + "Material assignment uses **triangle centroids**. The circle is approximated by selected triangles; its boundary is not a curved, conforming mesh.\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-07", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "# Geometry and mesh are specified through PhAST's standard configuration.\n", + "config[\"geometry\"][\"parameters\"].update(\n", + " length=L, height=H, origin=[0.0, -H / 2.0], nx=nx, ny=ny\n", + ")\n", + "regions = config[\"regions\"]\n", + "\n", + "# Keep the disk and both complementary outside-circle selectors identical.\n", + "for selector in (\n", + " regions[\"material_2_region\"],\n", + " regions[\"material_1_region\"][\"selectors\"][0],\n", + " regions[\"weak_midline_region\"][\"selectors\"][1],\n", + "):\n", + " selector.update(center=[circle_x, circle_y], radius=circle_radius)\n", + "\n", + "# Update the weak-layer rectangle and the rectangle used by its complement.\n", + "for rectangle in (\n", + " regions[\"weak_midline_region\"][\"selectors\"][0],\n", + " regions[\"material_1_region\"][\"selectors\"][1][\"selector\"],\n", + "):\n", + " rectangle.update(origin=[0.0, -layer_height / 2.0], size=[L, layer_height])\n", + "\n", + "# Starter crack and exterior boundary selections.\n", + "regions[\"initial_crack\"].update({\n", + " \"from\": [0.0, 0.0], \"to\": [initial_crack_length, 0.0],\n", + " \"thickness\": crack_selection_width,\n", + "})\n", + "regions[\"right_clamp\"] = {\n", + " \"kind\": \"node_region\", \"type\": \"boundary\", \"side\": \"right\", \"tolerance\": 1.0e-8,\n", + "}\n", + "regions[\"upper_load\"].update(\n", + " axis=\"x\", value=0.0, bounds={\"y\": [layer_height / 2.0, H / 2.0]}\n", + ")\n", + "regions[\"lower_load\"].update(\n", + " axis=\"x\", value=0.0, bounds={\"y\": [-H / 2.0, -layer_height / 2.0]}\n", + ")\n", + "\n", + "# Here h is the maximum triangle-edge length, not an equivalent-area diameter.\n", + "h = math.hypot(L / nx, H / ny)\n", + "print(f\"Expected mesh: {(nx + 1) * (ny + 1):,} nodes, {2 * nx * ny:,} triangles\")\n", + "print(f\"Maximum triangle edge h = {h:.5f} mm; h / ell_0 = {h / ell_0:.3f}\")\n", + "print(\"These sizes follow from the input; they are not completed simulation results.\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-08", + "metadata": {}, + "source": [ + "For the reference mesh, $h\\approx0.0451$ mm and $h/\\ell_0\\approx0.451$. The often-used criterion $h\\leq\\ell_0/2$ concerns resolution of the regularised crack. It does **not** establish mesh-independent forces, crack trajectories, or physical accuracy.\n", + "\n", + "Refining `nx` and `ny` at fixed $\\ell_0$ examines discretisation effects. Changing $\\ell_0$ changes the model as well as its resolution requirements. These are different experiments.\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-09", + "metadata": {}, + "source": [ + "## 5. Materials and the phase-field variable\n", + "\n", + "Displacement $\\mathbf u$ describes movement. Strain $\\boldsymbol\\varepsilon$ describes local deformation. Stress $\\boldsymbol\\sigma$ describes internal force per area. The scalar phase field $d$ describes damage: $d=0$ is intact and $d=1$ is fully damaged in this convention.\n", + "\n", + "For small strain,\n", + "\n", + "$$\n", + "\\boldsymbol\\varepsilon(\\mathbf u)=\\tfrac12\\left(\\nabla\\mathbf u+(\\nabla\\mathbf u)^{\\mathsf T}\\right).\n", + "$$\n", + "\n", + "The reference model degrades the selected crack-driving elastic energy using\n", + "\n", + "$$\n", + "g(d)=(1-d)^2+\\eta.\n", + "$$\n", + "\n", + "The small positive $\\eta$ avoids exactly zero residual stiffness. It is not a fourth material phase.\n", + "\n", + "A useful energy description, per unit out-of-plane thickness, is\n", + "\n", + "$$\n", + "\\Pi(\\mathbf u,d)=\\int_\\Omega\\left[\n", + "g(d)\\psi_0^+(\\boldsymbol\\varepsilon)+\\psi_0^-(\\boldsymbol\\varepsilon)\n", + "+G_c(\\mathbf x)\\left(\\frac{d^2}{2\\ell_0}+\\frac{\\ell_0}{2}|\\nabla d|^2\\right)\n", + "\\right]\\,\\mathrm d\\Omega-W_{\\mathrm{ext}}.\n", + "$$\n", + "\n", + "Here $\\Omega$ is the specimen domain, $\\psi_0^+$ is the elastic energy selected to drive damage, $\\psi_0^-$ is the retained part, and $W_{\\mathrm{ext}}$ is external work. The AT2 crack term spreads damage over a length controlled by $\\ell_0$.\n", + "\n", + "This input uses the **Amor volumetric-deviatoric energy split**. It is not a Miehe spectral split. The multi-material adapter does not support arbitrary replacement of the split or fracture law.\n", + "\n", + "The implementation uses a history field to retain the maximum previous crack-driving energy. That update is distinguished from unconstrained minimisation of the displayed instantaneous energy.\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-10", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "# Preserve the existing material definitions and assignment list.\n", + "for name, youngs_modulus, fracture_resistance in (\n", + " (\"material_1\", E, G_c_outer),\n", + " (\"weak_midline\", E, G_c_layer),\n", + " (\"material_2\", E_material_2, G_c_material_2),\n", + "):\n", + " parameters = config[\"materials\"][name][\"parameters\"]\n", + " parameters.update(\n", + " E=youngs_modulus, nu=nu, Gc=fracture_resistance, l0=ell_0,\n", + " eta_residual=eta, energy_split=\"amor\", pf_model=\"AT2\", plane_stress=True,\n", + " )\n", + " print(f\"{name:15s}: E = {youngs_modulus:g} MPa, Gc = {fracture_resistance:g} N/mm\")\n", + "\n", + "print(\"\\nAssignments are unchanged:\")\n", + "print(yaml.safe_dump(config[\"assignments\"], sort_keys=False))\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-11", + "metadata": {}, + "source": [ + "## 6. Initial crack, restraints, and opening displacements\n", + "\n", + "The initial crack is a selected row of nodes with maintained damage $d=1$. This defines the starting defect. It does not prescribe the later path around or through Material 2.\n", + "\n", + "| Selection | Applied condition | Reason |\n", + "|---|---|---|\n", + "| `right_clamp` | $u_x=u_y=0$ | Restrain the right edge and remove rigid motion. |\n", + "| `upper_load` | $u_y=+s\\Delta/2$ | Open the upper left arm. |\n", + "| `lower_load` | $u_y=-s\\Delta/2$ | Open the lower left arm. |\n", + "| `initial_crack` | Maintained $d=1$ | Keep the starter crack broken. |\n", + "\n", + "The load factor $s$ increases from zero to one. The final prescribed opening is $\\Delta$. Other exterior boundary segments have the natural zero-traction condition where no displacement is prescribed. The vertical loading conditions do not fix horizontal motion on the left segments.\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-12", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "# Keep the condition types, names, and active-condition list from the reference.\n", + "config[\"initial_conditions\"][0][\"value\"] = 1.0\n", + "for boundary in config[\"boundary_conditions\"]:\n", + " if boundary[\"name\"] == \"open_upper\":\n", + " boundary[\"value\"] = final_opening / 2.0\n", + " elif boundary[\"name\"] == \"open_lower\":\n", + " boundary[\"value\"] = -final_opening / 2.0\n", + "\n", + "print(yaml.safe_dump(config[\"boundary_conditions\"], sort_keys=False))\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-13", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "# Input schematic only. PhAST, not this drawing, constructs the actual FEM mesh.\n", + "figure, axis = plt.subplots(figsize=(11, 3.5), constrained_layout=True)\n", + "axis.add_patch(Rectangle((0.0, -H / 2.0), L, H, facecolor=\"#8AB0C4\", edgecolor=\"0.25\", label=\"Material 1: outer regions\"))\n", + "axis.add_patch(Rectangle(\n", + " (0.0, -layer_height / 2.0), L, layer_height,\n", + " facecolor=\"#D8D8D1\", edgecolor=\"0.45\", label=\"Weak bulk layer\",\n", + "))\n", + "axis.add_patch(Circle(\n", + " (circle_x, circle_y), circle_radius,\n", + " facecolor=\"#D6A24B\", edgecolor=\"0.15\", label=\"Material 2\",\n", + "))\n", + "axis.plot([0.0, initial_crack_length], [0.0, 0.0],\n", + " color=\"black\", linewidth=3, label=\"Prescribed starter crack\")\n", + "axis.plot([L, L], [-H / 2.0, H / 2.0],\n", + " color=\"black\", linewidth=4, label=\"Fixed right edge\")\n", + "load_y = (H + layer_height) / 4.0\n", + "for sign, label in ((1.0, \"+ opening\"), (-1.0, \"- opening\")):\n", + " axis.annotate(\n", + " \"\", xy=(0.0, sign * (load_y + 0.12 * H)), xytext=(0.0, sign * load_y),\n", + " arrowprops={\"arrowstyle\": \"->\", \"color\": \"0.15\", \"lw\": 1.8},\n", + " )\n", + " axis.text(-0.08 * L, sign * load_y, label, ha=\"right\", va=\"center\", fontsize=9)\n", + "axis.set(\n", + " xlim=(-0.30 * L, 1.03 * L), ylim=(-0.8 * H, 0.8 * H),\n", + " xlabel=\"x [mm]\", ylabel=\"y [mm]\",\n", + " title=\"Edited input: geometry and boundary-condition schematic\",\n", + ")\n", + "axis.set_aspect(\"equal\")\n", + "axis.legend(loc=\"upper center\", bbox_to_anchor=(0.5, -0.28), ncol=2, frameon=False)\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-14", + "metadata": {}, + "source": [ + "## 7. What does the staggered solver do?\n", + "\n", + "**Quasi-static** means inertia is omitted: each accepted load level seeks mechanical equilibrium. It does not mean explicit time stepping with a small time step. The load factor is not elapsed physical time.\n", + "\n", + "At one load level, PhAST alternates between two subproblems:\n", + "\n", + "1. Hold damage fixed and solve mechanical equilibrium for displacement. With no body force, the strong balance is $\\nabla\\!\\cdot\\boldsymbol\\sigma=\\mathbf0$.\n", + "2. Evaluate strain and the crack-driving elastic energy. Update the history $\\mathcal H_n=\\max(\\mathcal H_{n-1},\\psi_0^+(\\boldsymbol\\varepsilon_n))$.\n", + "3. Hold that driving field fixed and solve the spatial damage problem subject to bounds and maintained damage conditions.\n", + "4. Recompute the coupled state until the configured subproblem and staggered criteria are met. Otherwise, the run may reduce the load increment or stop.\n", + "\n", + "An **implicit** subproblem determines the unknown current field by satisfying equations or constrained optimality conditions. It does not add a prescribed damage increment at every node. Both mechanics and damage are implicit here. Explicit mechanics in a dynamic example is a different route, not a supported option of this multi-material DCB adapter.\n", + "\n", + "### Why is AT2 damage often called linear?\n", + "\n", + "For fixed $\\mathcal H$ and material coefficients, the unconstrained weak damage equation is\n", + "\n", + "$$\n", + "\\int_\\Omega\\left[\n", + "G_c\\ell_0\\nabla d\\cdot\\nabla q\n", + "+\\left(\\frac{G_c}{\\ell_0}+2\\mathcal H\\right)dq\n", + "\\right]\\,\\mathrm d\\Omega\n", + "=\\int_\\Omega2\\mathcal Hq\\,\\mathrm d\\Omega.\n", + "$$\n", + "\n", + "The test function $q$ represents an admissible variation of damage. With finite-element shape functions $N_i$, interpolation gives $A\\mathbf d=\\mathbf b$, where\n", + "\n", + "$$\n", + "A_{ij}=\\int_\\Omega\\left[\n", + "G_c\\ell_0\\nabla N_i\\cdot\\nabla N_j\n", + "+\\left(\\frac{G_c}{\\ell_0}+2\\mathcal H\\right)N_iN_j\n", + "\\right]\\,\\mathrm d\\Omega,\n", + "\\qquad\n", + "b_i=\\int_\\Omega2\\mathcal HN_i\\,\\mathrm d\\Omega.\n", + "$$\n", + "\n", + "Quadratic degradation gives a **linear unconstrained subproblem when the driving field is fixed**. It does not make the complete fracture calculation linear. Bounds $d_{\\mathrm{previous}}\\leq d\\leq1$, prescribed damage, spatial coupling, iterative linear algebra, and repeated staggered updates still require work.\n", + "\n", + "The projected residual measures the constrained damage problem, rather than requiring the unconstrained residual to vanish at every active bound. The matrix notation above explains the problem; this notebook delegates its implementation to PhAST instead of assembling a second FEM system.\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-15", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "# Load levels and solver iterations are different quantities.\n", + "step = config[\"analysis_steps\"][0]\n", + "step[\"controls\"][\"number_of_steps\"] = number_of_load_levels\n", + "\n", + "# Keep the supported route and the reference convergence controls.\n", + "# Changing this string alone does not enable dynamic multi-material fracture.\n", + "config[\"solver\"][\"type\"] = \"quasi_static\"\n", + "print(\"Analysis:\", step[\"type\"])\n", + "print(\"Nominal nonzero increments:\", number_of_load_levels - 1)\n", + "print(\"\\nNumerical controls inherited from the reference:\")\n", + "print(yaml.safe_dump(config[\"solver\"], sort_keys=False))\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-16", + "metadata": {}, + "source": [ + "The fine reference accepts 120 nonzero increments. With cutbacks, a run can accept more increments than the nominal count. Solver iteration limits count attempts to solve a subproblem, not additional physical load levels.\n", + "\n", + "Automatic cutback is not evidence of success by itself. Inspect accepted-step residuals and final status. Increasing an iteration limit or loosening a tolerance does not establish physical accuracy. This route reports its selected mechanics backend during execution; the requested value `auto` alone does not identify that backend.\n", + "\n", + "## 8. Save one YAML and use the standard commands\n", + "\n", + "The following cells write a new input under `runs/` and leave the checked-in reference unchanged. Each notebook execution receives a new directory so results from different parameter choices are not mixed.\n", + "\n", + "The command sequence is:\n", + "\n", + "~~~text\n", + "python -m phast explain-config CONFIG\n", + "python -m phast run CONFIG --validate-only\n", + "python -m phast run CONFIG --output_dir RUN\n", + "~~~\n", + "\n", + "Here `CONFIG` and `RUN` are placeholders. The code below supplies actual paths. Explanation describes the input without creating a mesh. Multi-material validation constructs the geometric mesh to check selections and assignments, but does not solve the FEM problem. Neither step establishes numerical convergence or physical validity.\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-17", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "# All new files are outside the curated example directory.\n", + "output_root = Path(os.environ.get(\"PHAST_NOTEBOOK_OUTPUT_ROOT\", str(ROOT / \"runs\")))\n", + "STUDY_DIR = output_root / f\"dcb_student_{uuid.uuid4().hex[:10]}\"\n", + "STUDY_DIR.mkdir(parents=True, exist_ok=False)\n", + "CONFIG_PATH = STUDY_DIR / \"config.yaml\"\n", + "RUN_DIR = STUDY_DIR / \"solution\"\n", + "\n", + "config[\"name\"] = \"Layered DCB student study\"\n", + "config[\"outputs\"][\"directory\"] = str(RUN_DIR.resolve())\n", + "for field in config[\"outputs\"][\"fields\"]:\n", + " if field[\"name\"] == \"trajectory\":\n", + " field.update(format=\"h5\", every=1)\n", + "CONFIG_PATH.write_text(yaml.safe_dump(config, sort_keys=False), encoding=\"utf-8\")\n", + "print(\"Student input:\", CONFIG_PATH)\n", + "print(\"New result directory, if solving is enabled:\", RUN_DIR)\n", + "print(\"The checked-in reference and retained results have not been modified.\")\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-18", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "# Use the notebook kernel's Python rather than a different shell environment.\n", + "def run_cli(arguments: list[str]) -> None:\n", + " completed = subprocess.run(\n", + " [sys.executable, \"-m\", \"phast\", *arguments],\n", + " cwd=ROOT, capture_output=True, text=True, check=False,\n", + " )\n", + " print(completed.stdout)\n", + " if completed.stderr:\n", + " print(completed.stderr)\n", + " if completed.returncode != 0:\n", + " raise RuntimeError(\"PhAST rejected the command. Read the first error above before proceeding.\")\n", + "\n", + "\n", + "run_cli([\"explain-config\", str(CONFIG_PATH)])\n", + "run_cli([\"run\", str(CONFIG_PATH), \"--validate-only\"])\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-19", + "metadata": {}, + "source": [ + "## 9. Optional: execute your edited problem\n", + "\n", + "With `RUN_NEW_SIMULATION = False`, the next cell does not start a simulation. It selects the separately retained reference results for the interpretation exercises.\n", + "\n", + "To compute new results, set that parameter to `True`, restart the kernel, and run all cells. This full reference calculation is not a two-second installation check. PhAST prints accepted increments, opening, reaction, crack-front position, and a residual as the solve progresses. A separate log is also saved.\n", + "\n", + "A failed solve raises an error here. Partial output must not be presented as a completed reference run. Use a new study directory when retrying; the trajectory writer does not silently overwrite an existing HDF5 file.\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-20", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "run_completed = False\n", + "if not RUN_NEW_SIMULATION:\n", + " print(\"No new FEM simulation requested. Later result cells use the retained fine reference.\")\n", + "else:\n", + " if RUN_DIR.exists():\n", + " raise FileExistsError(\"The result directory already exists. Restart and run all cells for a new study.\")\n", + " command = [\n", + " sys.executable, \"-m\", \"phast\", \"run\", str(CONFIG_PATH), \"--output_dir\", str(RUN_DIR),\n", + " ]\n", + " log_path = STUDY_DIR / \"solver.log\"\n", + " # Stream the existing runner's progress, without another solver loop.\n", + " with log_path.open(\"w\", encoding=\"utf-8\") as log:\n", + " with subprocess.Popen(\n", + " command, cwd=ROOT, stdout=subprocess.PIPE,\n", + " stderr=subprocess.STDOUT, text=True, bufsize=1,\n", + " ) as process:\n", + " if process.stdout is None:\n", + " raise RuntimeError(\"The solver output stream was not available.\")\n", + " for line in process.stdout:\n", + " print(line, end=\"\", flush=True)\n", + " log.write(line)\n", + " return_code = process.wait()\n", + " if return_code != 0:\n", + " raise RuntimeError(f\"The solve failed. Inspect {log_path}; partial output is not a converged run.\")\n", + " run_completed = True\n", + " print(\"\\nThe PhAST command completed. Inspect its status and residuals below.\")\n", + "\n", + "RESULT_DIR = RUN_DIR if run_completed else RETAINED\n", + "RESULT_LABEL = \"New student calculation\" if run_completed else \"Retained fine-reference calculation\"\n", + "print(\"Selected result:\", RESULT_LABEL)\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-21", + "metadata": {}, + "source": [ + "## 10. Read a result before drawing a conclusion\n", + "\n", + "A successful command and a visually plausible crack are different kinds of evidence. Inspect the recorded status, finite values, convergence measures, and the meaning of each observable.\n", + "\n", + "The response table is `history.csv`. Its `opening` column is relative arm opening in mm. The two-dimensional reaction is interpreted **per unit out-of-plane thickness**, in N/mm for this unit system. Comparison with an experimental force in N requires a specimen thickness and an appropriate two-dimensional idealisation.\n", + "\n", + "The initial crack already has $d=1$. Therefore, `maximum_damage = 1` alone says nothing about new crack growth.\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-22", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "def read_history(path: Path) -> dict[str, np.ndarray]:\n", + " required = (\"opening\", \"reaction_force\", \"connected_crack_front_x\", \"projected_damage_residual\")\n", + " with path.open(\"r\", encoding=\"utf-8\", newline=\"\") as stream:\n", + " rows = list(csv.DictReader(stream))\n", + " if not rows:\n", + " raise ValueError(f\"No accepted steps were found in {path.name}.\")\n", + " history = {\n", + " key: np.asarray([float(row[key]) for row in rows], dtype=float)\n", + " for key in required\n", + " }\n", + " if not all(np.isfinite(values).all() for values in history.values()):\n", + " raise ValueError(\"A selected history contains non-finite values. Investigate rather than hiding them.\")\n", + " return history\n", + "\n", + "\n", + "history = read_history(RESULT_DIR / \"history.csv\")\n", + "summary = json.loads((RESULT_DIR / \"summary.json\").read_text(encoding=\"utf-8\"))\n", + "if not summary.get(\"all_steps_converged\", False):\n", + " raise RuntimeError(\"The selected run does not report convergence of all accepted steps.\")\n", + "\n", + "for key in (\n", + " \"n_nodes\", \"n_elements\", \"accepted_steps\", \"cutbacks\", \"final_opening\",\n", + " \"final_reaction_force\", \"final_connected_crack_front_x\", \"maximum_projected_damage_residual\",\n", + "):\n", + " print(f\"{key}: {summary[key]}\")\n", + "print(\"These are recorded numerical results, not experimental validation.\")\n", + "\n", + "figure, axes = plt.subplots(1, 3, figsize=(13, 3.5), constrained_layout=True)\n", + "axes[0].plot(history[\"opening\"], history[\"reaction_force\"], color=\"0.20\")\n", + "axes[0].set_ylabel(\"Opening reaction per unit thickness [N/mm]\")\n", + "axes[1].step(history[\"opening\"], history[\"connected_crack_front_x\"], where=\"post\", color=\"0.20\")\n", + "axes[1].set_ylabel(\"Maximum seed-connected x at d >= 0.80 [mm]\")\n", + "axes[2].plot(history[\"opening\"], history[\"projected_damage_residual\"], color=\"0.20\")\n", + "result_input = config if run_completed else reference_config\n", + "axes[2].axhline(\n", + " result_input[\"solver\"][\"projected_residual_tolerance\"],\n", + " color=\"0.45\", linestyle=\"--\", label=\"Configured projected-residual tolerance\",\n", + ")\n", + "axes[2].set_ylabel(\"Projected damage residual [solver convention]\")\n", + "axes[2].legend(fontsize=8, frameon=False)\n", + "for axis in axes:\n", + " axis.set_xlabel(\"Opening [mm]\")\n", + " axis.grid(alpha=0.2)\n", + "figure.suptitle(RESULT_LABEL)\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-23", + "metadata": {}, + "source": [ + "### Reading the retained fine-reference result\n", + "\n", + "The retained fine calculation reports 120 accepted increments, no cutbacks, and 0.30 mm final opening. Its maximum recorded projected damage residual is about $6.91\\times10^{-6}$, below the configured $10^{-5}$ criterion. These are recorded results, not a new solve by this notebook.\n", + "\n", + "The seed-connected front reaches $x=2.0625$ mm at threshold $d\\geq0.80$. This is the largest $x$ coordinate connected to the starter crack through damaged mesh edges, not the crack's arc length. Passing the circle's leftmost $x$ coordinate does not establish penetration: a branch can advance above or below it without entering the disk.\n", + "\n", + "The retained field review found no seed-connected $d\\geq0.80$ node inside the tough disk. The images show branching near its flanks, not completed bypass. Both geometry and material contrast influence the result.\n", + "\n", + "## 11. Displacement, strain, stress, and damage images\n", + "\n", + "**This section always shows the retained fine reference.** The figures are not updated by editing the parameter cell. A separate cell below displays new images only after a fresh solve completes.\n", + "\n", + "Displacement is motion in mm. Strain is dimensionless deformation. Stress components are in MPa. Damage is dimensionless. Read component names and any deformation scaling on each figure. Colour intensity is not comparable between different physical quantities.\n", + "\n", + "![Retained fine-reference material regions and property fields.](../../examples/two_material_dcb_beta/results/material_fields.png)\n", + "\n", + "![Retained fine-reference displacement fields at the final opening.](../../examples/two_material_dcb_beta/results/displacement_final.png)\n", + "\n", + "![Retained fine-reference strain fields at the final opening.](../../examples/two_material_dcb_beta/results/strain_final.png)\n", + "\n", + "![Retained fine-reference stress fields at the final opening.](../../examples/two_material_dcb_beta/results/stress_final.png)\n", + "\n", + "![Retained fine-reference damage at 0.30 mm opening; branching does not establish completed bypass.](../../examples/two_material_dcb_beta/results/damage_final.png)\n", + "\n", + "### Follow the crack evolution\n", + "\n", + "Identify the prescribed starter crack first, then newly damaged material. Grey outlines locate material boundaries; they are not additional crack constraints. This is a quasi-static load sequence, so animation speed is not physical crack speed.\n", + "\n", + "![Retained fine-reference damage evolution under increasing opening, without demonstrated penetration or completed bypass.](../../examples/two_material_dcb_beta/results/damage_evolution.gif)\n", + "\n", + "The final damage image is the static fallback if the animation does not play.\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-24", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "# New figures are generated by PhAST's existing example runner.\n", + "if run_completed:\n", + " display(Markdown(\"### New student calculation: generated field images\"))\n", + " for filename in (\n", + " \"initial_conditions.png\", \"material_fields.png\", \"displacement_final.png\",\n", + " \"strain_final.png\", \"stress_final.png\", \"damage_final.png\", \"damage_evolution.gif\",\n", + " ):\n", + " display(Markdown(f\"**{filename}**: new calculation, not the retained reference.\"))\n", + " display(Image(filename=str(RUN_DIR / filename)))\n", + "else:\n", + " print(\"No fresh-run figures exist for the edited input. The images above are the retained reference.\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-25", + "metadata": {}, + "source": [ + "## 12. Inspect numerical fields through PhAST\n", + "\n", + "A new run requests `training_data.h5`. It contains stored fields and mesh data and remains local, rather than being added to the public repository. The compact retained directory intentionally omits that full trajectory.\n", + "\n", + "Use `phast.load_result` for a completed local run. Its field reader returns stored quantities; it does not infer raw fields from PNG images. A derived scalar such as von Mises stress requires the documented post-processing convention, especially for a two-dimensional model.\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-26", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "if run_completed:\n", + " result = phast.load_result(RUN_DIR)\n", + " print(\"Available stored fields:\", result.field_names())\n", + " for field_name in (\"displacement\", \"strain\", \"stress\", \"damage\"):\n", + " values = np.asarray(result.field(field_name)) # Final stored step by default.\n", + " if not np.isfinite(values).all():\n", + " raise ValueError(f\"Non-finite values in {field_name} require investigation.\")\n", + " print(field_name, \"shape:\", values.shape, \"dtype:\", values.dtype)\n", + " damage = np.asarray(result.field(\"damage\"))\n", + " if damage.min() < -1.0e-8 or damage.max() > 1.0 + 1.0e-8:\n", + " raise ValueError(\"The final damage field lies outside its expected numerical bounds.\")\n", + "else:\n", + " print(\"Raw fine-reference HDF5 fields are not part of the compact public results.\")\n", + " print(\"Run a new calculation before using this field-loading section.\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-27", + "metadata": {}, + "source": [ + "## 13. What changes without the circle's toughness contrast?\n", + "\n", + "A useful comparison changes one physical input while holding the other choices fixed. Two complete existing inputs are available:\n", + "\n", + "- `examples/two_material_dcb_beta/comparisons/tough_region.yaml`: disk $G_c=0.12$ N/mm;\n", + "- `examples/two_material_dcb_beta/comparisons/uniform_layer.yaml`: disk $G_c=0.04$ N/mm, equal to the weak layer.\n", + "\n", + "Both use a **separate coarse $120\\times48$ mesh**, not the fine grid above. Outer material remains fracture-resistant: the control is a uniform weak layer, not a homogeneous specimen. Elastic properties and prescribed loading are held fixed.\n", + "\n", + "At 0.30 mm opening, the retained tough-region case has about 0.60938 N/mm reaction and a connected front at 2.10 mm. The control has about 0.45847 N/mm and a front at 2.40 mm. Both cases branch. The circle cannot therefore be identified as the sole cause of branching.\n", + "\n", + "The control required three cutbacks and has different accepted openings. Compare each history against its own opening coordinate, not by subtracting corresponding CSV rows. The next command only replots retained data. It does **not** run either solver or use your edited input.\n" + ] + }, + { + "cell_type": "code", + "id": "dcb-28", + "metadata": {}, + "execution_count": null, + "outputs": [], + "source": [ + "# Reuse the checked-in comparison script and its provenance/checksum checks.\n", + "comparison_directory = STUDY_DIR / \"retained_comparison\"\n", + "comparison = subprocess.run(\n", + " [\n", + " sys.executable, str(EXAMPLE / \"compare_results.py\"),\n", + " \"--output-dir\", str(comparison_directory),\n", + " ],\n", + " cwd=ROOT, capture_output=True, text=True, check=False,\n", + ")\n", + "print(comparison.stdout)\n", + "if comparison.returncode != 0:\n", + " raise RuntimeError(comparison.stderr or \"The retained comparison could not be replotted.\")\n", + "for filename in (\"matched_control_comparison.png\", \"matched_control_damage.png\"):\n", + " display(Markdown(\"**Retained coarse comparison. Not a new calculation of the edited problem.**\"))\n", + " display(Image(filename=str(comparison_directory / filename)))\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-29", + "metadata": {}, + "source": [ + "The comparison retains the original triangle connectivity. The coarse tough-region run did not write HDF5; the uniform-layer control did, and its stored mesh supplied that connectivity. Both new comparison inputs request HDF5 for local reruns. Exact historical generating-source identity remains unattested, as recorded in `results/comparison_provenance.json`.\n", + "\n", + "To compute new coarse comparisons, submit each existing YAML through the **same** explain/check/run sequence, with distinct output directories. Do not compare a newly edited fine run directly with a retained coarse result and attribute the difference solely to material properties.\n", + "\n", + "## 14. Guided student experiments\n", + "\n", + "Predict an outcome first. Afterwards, record the input, source version, solver status, mesh measure, and result directory.\n", + "\n", + "| Exercise | Change | Keep fixed | Question |\n", + "|---|---|---|---|\n", + "| Remove toughness contrast | Set `G_c_material_2 = G_c_layer`. | Geometry, elastic properties, mesh, loading, and solver controls. | Does the crack advance farther at the same opening? Do the fields support your explanation? |\n", + "| Separate stiffness from toughness | Vary `E_material_2` at fixed `G_c_material_2`. | Other properties and loading. | How does elastic redistribution change the crack-driving field? |\n", + "| Examine discretisation | Increase `nx` and `ny`, retaining an even `ny`. | Geometry, $\\ell_0$, material properties, and loading. | How much do reaction and connected-front position change? One pair is not a completed convergence study. |\n", + "| Examine the model length | Change `ell_0` and choose adequate mesh resolution. | State all other choices. | Which changes may reflect regularisation rather than discretisation error? |\n", + "| Move the obstacle | Change `circle_x` or `circle_radius`. | Decide whether the disk must remain wholly inside the layer. | Does the starter crack end before the disk? Are the selections meaningful? |\n", + "\n", + "For a short laboratory session, interpret the retained results first. Plan long solves separately instead of making participants wait for a full fracture calculation before understanding the model.\n", + "\n", + "### Questions for a short written explanation\n", + "\n", + "1. Why does the reference isolate a toughness contrast rather than a stiffness contrast inside the layer?\n", + "2. Why does $\\max d=1$ not establish crack growth?\n", + "3. What is fixed when the damage subproblem is described as linear?\n", + "4. Why can 121 nominal load levels lead to more than 120 accepted increments?\n", + "5. Why is a weak finite-width bulk layer not a cohesive interface law?\n", + "6. Which independent data would be needed before calling this an experimentally validated specimen model?\n" + ] + }, + { + "cell_type": "markdown", + "id": "dcb-30", + "metadata": {}, + "source": [ + "## 15. Troubleshooting and limits\n", + "\n", + "| Observation | First action |\n", + "|---|---|\n", + "| `ModuleNotFoundError` or an unexpected version | Select the environment used for `pip install -e .`, confirm the checkout printed by the first cell, and restart the kernel. |\n", + "| Reference configuration not found | Use the complete repository, not a single downloaded notebook. Set `REPOSITORY_PATH` if necessary. |\n", + "| Empty crack, load, or material selection | Inspect dimensions and selectors. An odd `ny` may omit midline nodes. A small disk may contain no triangle centroids. |\n", + "| Geometry changes alter the material arrangement | Check circle/layer overlap, the starter tip, and load segments. Accepted input does not establish that the model answers your intended question. |\n", + "| Repeated cutbacks or solver failure | Preserve the log. Inspect units, regions, tolerances, and residuals. Do not hide non-finite values or present partial fields as a completed run. |\n", + "| Damage is already one before visible extension | The starter crack is prescribed. Inspect the connected front and newly damaged material instead. |\n", + "| Retained HDF5 file missing | This is intentional. A fresh local run creates the requested trajectory. |\n", + "| Result directory already exists | Restart and run all cells for a new study identifier. Do not overwrite another parameter choice. |\n", + "| OneDrive reports conflicts | Finish writing before synchronisation elsewhere. HDF5 reduces file count but does not make concurrent writers safe. |\n", + "\n", + "This example does not establish a general interface, plasticity, three-dimensional composite, dynamic multi-material, or learned-damage workflow. A neural predictor cannot be inserted into this particular adapter merely by changing its solver string. Check the [capability matrix](../user_guide/capability_matrix.md) before selecting another route.\n", + "\n", + "The small-strain assumption, diffuse crack length relative to the circular region, elementwise boundary approximation, and confinement by the outer material affect interpretation. A convincing animation alone does not settle those modelling questions.\n", + "\n", + "## Ask for help and contribute\n", + "\n", + "If an installation, configuration, execution, or interpretation step is unclear, [open a GitHub issue](https://github.com/CEMS-Lab/PhAST/issues/new/choose). Include the YAML, command, operating system, Python and PhAST versions, source revision if known, and complete first error or solver log. State whether you inspected retained data or ran a new calculation.\n", + "\n", + "Students, developers, and researchers are welcome to report unclear wording, propose a controlled comparison, review a result, or contribute a documented improvement. No familiarity with solver internals is required to ask a question.\n", + "\n", + "See [One configuration workflow for student simulations](07_standard_simulation_workflow.md) for the common input convention. The [DCB example on GitHub](https://github.com/CEMS-Lab/PhAST/tree/main/examples/two_material_dcb_beta) contains the full configuration, result README, and evidence limitations.\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "file_extension": ".py", + "mimetype": "text/x-python", + "pygments_lexer": "ipython3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/two_material_dcb_beta/README.md b/examples/two_material_dcb_beta/README.md index 3a88454ec..bfac290d2 100644 --- a/examples/two_material_dcb_beta/README.md +++ b/examples/two_material_dcb_beta/README.md @@ -108,8 +108,51 @@ The initial crack is maintained by a damage constraint. No damaged path is prescribed through or around Material 2. Penetration, deflection, branching, or arrest must be identified from the computed fields rather than assumed. +## Material Regions And Input Conditions + +![Input diagram showing blue outer material, grey weak layer, ochre circular Material 2, the starter crack, right clamp, and prescribed opening displacements.](results/material_regions_and_loading.png) + +The example contains **three bulk material regions**, despite its historical +folder name. Blue denotes the outer Material 1, grey denotes the weak layer, +and ochre denotes the circular Material 2. All three have the same Young's +modulus, $E = 5{,}000$ MPa. Their fracture toughnesses are respectively +$G_c = 1.00$, $0.04$, and $0.12$ N/mm. Material 2 is therefore tougher than +the weak layer, not stiffer. The colours describe the input assignments; +they are not the red damage scale used in the result figures. + +The plate measures $6.0 \times 1.2$ mm. The weak layer is $0.6$ mm high, +the circular region has radius $0.20$ mm, and the initial crack is $1.5$ mm +long. The right edge is fixed in both directions. The upper and lower +left-edge loading segments receive vertical displacements of $+0.15$ mm +and $-0.15$ mm at the final load factor, giving a relative opening of +$0.30$ mm. Horizontal displacement is not prescribed on these loading +segments. Other tractions are zero. The material boundaries do not have +an independent cohesive or interface law. + +After installing PhAST, run this command from the repository root to draw +the setup before solving: + +```bash +# Read the YAML, build its mesh and selections, and draw the inputs only. +python examples/two_material_dcb_beta/plot_setup.py \ + --config examples/two_material_dcb_beta/config.yaml \ + --output-dir runs/dcb_setup +``` + +The output is `runs/dcb_setup/material_regions_and_loading.png`. Use +`--config` with your edited YAML to inspect a modified setup. This command +uses the same PhAST mesh, material assignments, and node selections as the +solver; it neither solves the fracture problem nor establishes convergence. +New full calculations also write `material_regions_and_loading.png`. + ## Standard Results +For a step-by-step explanation, use the +[layered DCB student notebook](../../docs/tutorial/notebook_layered_dcb.ipynb). +It connects the equations to editable geometry, materials, loads, solver +settings, and result interpretation. Retained images and animation are the +default; running a new simulation is an explicit choice. + Retained lightweight evidence is in [results/README.md](results/README.md). It is not a new run of the documentation commands: diff --git a/examples/two_material_dcb_beta/plot_setup.py b/examples/two_material_dcb_beta/plot_setup.py new file mode 100644 index 000000000..56e02c066 --- /dev/null +++ b/examples/two_material_dcb_beta/plot_setup.py @@ -0,0 +1,36 @@ +"""Draw the DCB input configuration without running a fracture calculation.""" +from __future__ import annotations + +import argparse +from pathlib import Path + +from phast.workflow import problem_spec_from_yaml +from phast.workflow.multimaterial_fracture import ( + _assign_materials, + _build_mesh, + _resolve_bcs, + _resolve_regions, +) +from phast.workflow.setup_plot import save_material_setup + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--config", type=Path, default=Path(__file__).with_name("config.yaml")) + parser.add_argument("--output-dir", type=Path, default=Path("runs/dcb_setup")) + args = parser.parse_args() + spec = problem_spec_from_yaml(args.config) + # Reuse the solver's mesh and selections, rather than a second model setup. + nodes, elements = _build_mesh(spec) + node_regions, element_regions = _resolve_regions(spec, nodes, elements) + _, _, material_ids, material_names = _assign_materials(spec, element_regions, len(elements)) + path = save_material_setup( + spec, nodes, elements, material_ids, material_names, node_regions, + _resolve_bcs(spec, node_regions), + args.output_dir / "material_regions_and_loading.png", + ) + print(f"Input diagram saved to {path}; no fracture simulation was run.") + + +if __name__ == "__main__": + main() diff --git a/examples/two_material_dcb_beta/results/README.md b/examples/two_material_dcb_beta/results/README.md index 8936a3afd..837fc713e 100644 --- a/examples/two_material_dcb_beta/results/README.md +++ b/examples/two_material_dcb_beta/results/README.md @@ -12,6 +12,27 @@ public documentation dependency. The generating source fingerprint was not recorded at runtime and remains unknown. A later public-availability fingerprint does not establish which source generated historical outputs. +## Material Regions And Input Conditions + +![Reference input diagram with three coloured bulk regions, a black starter crack, a fixed right edge, and opposite vertical displacements at the left edge.](material_regions_and_loading.png) + +This new diagram represents the reference YAML inputs, not a new fracture +calculation. Blue is the outer Material 1, grey is the weak bulk layer, and +ochre is the circular Material 2. All regions have $E = 5{,}000$ MPa; their +fracture toughnesses are $1.00$, $0.04$, and $0.12$ N/mm, respectively. +Material 2 is tougher than the weak layer but has the same elastic modulus. + +The right edge has $u_x = u_y = 0$. The upper and lower loading segments +on the left edge reach $u_y = +0.15$ mm and $u_y = -0.15$ mm, respectively. +Their horizontal displacement remains free. The final relative opening is +$0.30$ mm. Black marks on the centreline identify the prescribed initial +crack, not newly propagated damage. Other tractions are zero. + +Colours show bulk assignments, not a cohesive interface or a computed damage +field. The [example instructions](../README.md#material-regions-and-input-conditions) +explain how to regenerate the diagram after editing the geometry, materials, +or boundary conditions. Retained simulation figures below remain unchanged. + ## What The Calculation Demonstrates At a total opening of 0.30 mm, the seeded crack bifurcates within the weak diff --git a/examples/two_material_dcb_beta/results/material_regions_and_loading.png b/examples/two_material_dcb_beta/results/material_regions_and_loading.png new file mode 100644 index 000000000..53af5638a Binary files /dev/null and b/examples/two_material_dcb_beta/results/material_regions_and_loading.png differ diff --git a/examples/two_material_dcb_beta/results/visual_manifest.json b/examples/two_material_dcb_beta/results/visual_manifest.json index 27da8ff1c..8b9642f47 100644 --- a/examples/two_material_dcb_beta/results/visual_manifest.json +++ b/examples/two_material_dcb_beta/results/visual_manifest.json @@ -67,5 +67,14 @@ { "path": "matched_control_damage.png", "artifact_type": "image" + }, + { + "path": "material_regions_and_loading.png", + "artifact_type": "image", + "source": "../config.yaml", + "width_px": 1800, + "height_px": 1200, + "generation_command": "python examples/two_material_dcb_beta/plot_setup.py --config examples/two_material_dcb_beta/config.yaml --output-dir examples/two_material_dcb_beta/results", + "claim_role": "Input representation only; no new fracture calculation." } ] diff --git a/src/phast/workflow/multimaterial_fracture.py b/src/phast/workflow/multimaterial_fracture.py index 097290fca..6a679f7f0 100644 --- a/src/phast/workflow/multimaterial_fracture.py +++ b/src/phast/workflow/multimaterial_fracture.py @@ -28,6 +28,7 @@ from ..mesh import FEMMesh from ..io_utils import init_h5, write_h5_snapshot from .specs import ProblemSpec +from .setup_plot import save_material_setup class MultimaterialWorkflowError(RuntimeError): @@ -870,7 +871,11 @@ def run_multimaterial_fracture_spec(spec: ProblemSpec, *, output_dir=None) -> in "config_sha256": hashlib.sha256(config_text.encode()).hexdigest(), "source_config": spec.source_path, "solver_parameters": dict(spec.solver.parameters)}, indent=2) + "\n", encoding="utf-8") - visuals = ["initial_conditions.png", "material_fields.png", "damage_final.png", + save_material_setup( + spec, nodes_np, elements_np, material_ids, material_names, node_regions, + resolved_bcs, output / "material_regions_and_loading.png", + ) + visuals = ["initial_conditions.png", "material_regions_and_loading.png", "material_fields.png", "damage_final.png", "displacement_final.png", "strain_final.png", "stress_final.png", "load_displacement.png", "convergence_and_crack_front.png", "damage_evolution.gif"] diff --git a/src/phast/workflow/setup_plot.py b/src/phast/workflow/setup_plot.py new file mode 100644 index 000000000..fade3415f --- /dev/null +++ b/src/phast/workflow/setup_plot.py @@ -0,0 +1,227 @@ +"""Input-only material and boundary-condition diagrams for 2D teaching cases. + +The caller supplies the mesh, material assignments, and boundary selections +resolved by the calculation. No displacement or damage solution is computed. +""" +from __future__ import annotations + +from collections.abc import Sequence +from pathlib import Path +from typing import Any + +import matplotlib.pyplot as plt +from matplotlib.colors import BoundaryNorm, ListedColormap +import numpy as np + +from .specs import BoundaryConditionSpec, ProblemSpec + + +MATERIAL_COLOURS = { + "material_1": "#8AB0C4", + "weak_midline": "#D8D8D1", + "material_2": "#D6A24B", +} +MATERIAL_LABELS = { + "material_1": "Material 1: outer regions", + "weak_midline": "Weak bulk layer", + "material_2": "Material 2: circular region", +} + + +def save_material_setup( + spec: ProblemSpec, + nodes: np.ndarray, + elements: np.ndarray, + material_ids: np.ndarray, + material_names: Sequence[str], + node_regions: dict[str, np.ndarray], + boundary_conditions: Sequence[tuple[BoundaryConditionSpec, Any]], + output_path: str | Path, +) -> Path: + """Draw the assigned T3 materials and selected nodes in reference coordinates. + + This diagram is intended for rectangular 2D teaching geometries. Material + units are labelled for the N-mm convention when geometry units are mm; + otherwise the table explicitly refers to the configuration's input units. + Arrow lengths indicate direction only, not the displacement magnitude. + """ + output_path = Path(output_path) + output_path.parent.mkdir(parents=True, exist_ok=True) + xmin, ymin = nodes.min(axis=0)[:2] + xmax, ymax = nodes.max(axis=0)[:2] + length, height = xmax - xmin, ymax - ymin + units = spec.geometry.units if spec.geometry is not None else "input units" + colours = [ + MATERIAL_COLOURS.get(name, f"C{index % 10}") + for index, name in enumerate(material_names) + ] + material_by_name = {material.name: material for material in spec.materials} + colour_map = ListedColormap(colours) + colour_norm = BoundaryNorm(np.arange(len(colours) + 1) - 0.5, len(colours)) + + # Apply the public figure style locally, without changing later field plots. + with plt.rc_context(): + plt.style.use("default") + plt.rcParams.update({ + "font.family": "serif", + "font.serif": ["STIXGeneral", "DejaVu Serif"], + "mathtext.fontset": "stix", + "text.usetex": False, + "axes.unicode_minus": False, + "font.size": 11, + "figure.facecolor": "white", + "axes.facecolor": "white", + "text.color": "black", + }) + figure = plt.figure(figsize=(12, 8)) + figure.text(0.075, 0.95, "Layered DCB: materials and input conditions", + fontsize=21, weight="bold") + figure.text(0.075, 0.912, + "Reference geometry and prescribed data, not a computed deformation or crack path.", + fontsize=12) + axis = figure.add_axes((0.075, 0.49, 0.86, 0.34)) + axis.tripcolor(nodes[:, 0], nodes[:, 1], elements, + facecolors=material_ids, shading="flat", + cmap=colour_map, norm=colour_norm, rasterized=True) + axis.plot([xmin, xmax, xmax, xmin, xmin], + [ymin, ymin, ymax, ymax, ymin], color="black", linewidth=0.8) + + # The black seed is drawn from the same node selection used by PhAST. + for initial in spec.initial_conditions: + if initial.field != "damage" or initial.region not in node_regions: + continue + indices = np.flatnonzero(node_regions[initial.region]) + if not indices.size: + continue + points = nodes[indices, :2] + axis.scatter(points[:, 0], points[:, 1], s=8, color="black", zorder=4) + span = np.ptp(points[:, 0]) + axis.annotate( + f"Initial crack: $d={float(initial.value):g}$\n" + f"Selected length = {span:g} {units}", + xy=points.mean(axis=0), + xytext=(xmin + 0.22 * length, ymin - 0.40 * height), + ha="center", va="center", fontsize=10, + arrowprops={"arrowstyle": "-", "color": "black", "lw": 0.8}, + ) + + # Merge x/y constraints on the same edge into one readable annotation. + fixed_regions: dict[str, set[int]] = {} + vertical_values: list[float] = [] + for boundary, indices in boundary_conditions: + components = (0, 1) if boundary.component is None else (boundary.component,) + if boundary.kind == "fix": + fixed_regions.setdefault(boundary.region, set()).update(components) + continue + if boundary.kind != "prescribe": + continue + points = nodes[np.asarray(indices, dtype=int), :2] + centre = points.mean(axis=0) + axis.scatter(points[:, 0], points[:, 1], s=9, color="black", zorder=4) + value = float(boundary.value) + for component in components: + direction = 1.0 if value >= 0.0 else -1.0 + start = centre.copy() + start[0] -= 0.035 * length + end = start.copy() + end[component] += direction * 0.22 * height + axis.annotate("", xy=end, xytext=start, + arrowprops={"arrowstyle": "-|>", "lw": 1.8, "color": "black"}) + symbol = "x" if component == 0 else "y" + axis.text( + xmin - 0.075 * length, centre[1] + direction * 0.14 * height, + f"{boundary.name or boundary.region}\n$u_{symbol}={value:+g}$ {units}", + ha="right", va="center", fontsize=11, + ) + if component == 1: + vertical_values.append(value) + for region, components in fixed_regions.items(): + points = nodes[np.flatnonzero(node_regions[region]), :2] + centre = points.mean(axis=0) + axis.scatter(points[:, 0], points[:, 1], marker="+", s=25, + linewidths=1.2, color="black", zorder=5) + label = ", ".join(f"$u_{'x' if component == 0 else 'y'}=0$" + for component in sorted(components)) + axis.annotate( + f"{region.replace('_', ' ')}\n{label}", xy=centre, + xytext=(xmax + 0.055 * length, centre[1]), ha="left", va="center", + fontsize=11, arrowprops={"arrowstyle": "-", "lw": 0.8, "color": "black"}, + ) + + # Label the declared circular geometry; its material assignment remains + # the element-centre selection shown by the coloured triangles. + for region in spec.regions: + selector = region.selector + if selector.get("type") != "circle": + continue + centre = np.asarray(selector["center"], dtype=float) + radius = float(selector["radius"]) + axis.annotate( + f"Circular region: $R={radius:g}$ {units}", xy=centre, + xytext=(centre[0] + 0.16 * length, ymax + 0.15 * height), + ha="left", va="center", fontsize=10, + arrowprops={"arrowstyle": "-", "lw": 0.8, "color": "black"}, + ) + dimension_y = ymax + 0.48 * height + axis.annotate("", xy=(xmin, dimension_y), xytext=(xmax, dimension_y), + arrowprops={"arrowstyle": "<->", "lw": 0.8, "color": "black"}) + axis.text((xmin + xmax) / 2.0, dimension_y + 0.035 * height, + f"Length = {length:g} {units}", ha="center", va="bottom", fontsize=10) + axis.set(xlim=(xmin - 0.32 * length, xmax + 0.26 * length), + ylim=(ymin - 0.60 * height, ymax + 0.75 * height), + xlabel=f"x [{units}]", ylabel=f"y [{units}]", aspect="equal") + axis.set_yticks([ymin, (ymin + ymax) / 2.0, ymax]) + axis.spines[["top", "right"]].set_visible(False) + + table_axis = figure.add_axes((0.075, 0.275, 0.86, 0.155)) + table_axis.set_axis_off() + modulus_label = "$E$ [MPa]" if units == "mm" else "$E$ [input units]" + toughness_label = "$G_c$ [N/mm]" if units == "mm" else "$G_c$ [input units]" + rows = [] + for name in material_names: + parameters = material_by_name[name].parameters + rows.append([ + MATERIAL_LABELS.get(name, name.replace("_", " ")), + f"{float(parameters['E']):g}", + f"{float(parameters['Gc']):g}", + ]) + table = table_axis.table(cellText=rows, + colLabels=["Material region", modulus_label, toughness_label], + colWidths=[0.56, 0.22, 0.22], cellLoc="left", + colLoc="left", bbox=(0, 0, 1, 1)) + table.auto_set_font_size(False) + table.set_fontsize(12) + for (row, column), cell in table.get_celld().items(): + cell.set_edgecolor("white") + cell.set_linewidth(2) + if row == 0: + cell.set_facecolor("#F2F2F0") + cell.set_text_props(weight="bold") + elif column == 0: + cell.set_facecolor(colours[row - 1]) + + first_material = spec.materials[0].parameters + step = spec.analysis_steps[0] + increments = max(int(step.controls.get("number_of_steps", 2)) - 1, 1) + opening = "" + if vertical_values and min(vertical_values) < 0 < max(vertical_values): + opening = (f"Final relative opening: {max(vertical_values) - min(vertical_values):g} {units}. " + "Prescribed values above are at load factor 1.\n") + figure.text( + 0.075, 0.229, + f"Geometry: {length:g} x {height:g} {units}. " + f"Mesh: {len(nodes):,} nodes; {len(elements):,} linear triangles.\n" + f"Shared inputs: $\\nu={float(first_material['nu']):g}$; " + f"$\\ell_0={float(first_material['l0']):g}$ {units}. " + f"Analysis: {step.kind.replace('_', ' ')}; {increments} nominal load increments.", + fontsize=11, va="top", linespacing=1.6, + ) + figure.text(0.075, 0.136, opening + + "Material colours identify bulk regions, not damage. Unprescribed tractions are zero.", + fontsize=11, va="top", linespacing=1.6) + figure.text(0.075, 0.045, + "Source: YAML selections resolved on the PhAST mesh. Arrow lengths are schematic. " + "No simulation was run to produce this diagram.", fontsize=9) + figure.savefig(output_path, dpi=150, facecolor="white") + plt.close(figure) + return output_path