From geometry to crack propagation: a square-plate SENT simulation¶

from __future__ import annotations
import copy
import os
from pathlib import Path
import re
import shutil
import subprocess
import sys

from IPython.display import HTML, Image, Markdown, display
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
from matplotlib.animation import FuncAnimation, PillowWriter
import numpy as np
import yaml


def find_repository_root(start: Path) -> Path:
    # Locate a source checkout from the notebook working directory.
    start = start.resolve()
    for candidate in (start, *start.parents):
        if (candidate / "pyproject.toml").is_file() and (candidate / "examples").is_dir():
            return candidate
    raise RuntimeError(
        "Run this notebook from a cloned PhAST repository. "
        "See the installation commands above."
    )

ROOT = find_repository_root(Path.cwd())
os.chdir(ROOT)
sys.path.insert(0, str(ROOT / "src"))
os.environ["PYTHONPATH"] = str(ROOT / "src") + os.pathsep + os.environ.get("PYTHONPATH", "")

import phast
from phast.core.mesh import FEMMesh
from phast.physics.material import Material

plt.rcParams.update({"font.family": "STIXGeneral", "mathtext.fontset": "stix"})

def require_finite(values, name):
    """Reject empty or non-finite data rather than disguising a failed solve."""
    array = np.asarray(values)
    if array.size == 0 or not np.isfinite(array).all():
        raise ValueError(
            f"Empty or non-finite values in {name}; inspect the solver log "
            "before plotting or interpreting this run."
        )
    return array

output_root = Path(os.environ.get("PHAST_NOTEBOOK_OUTPUT_ROOT", ROOT / "runs"))
base = output_root.expanduser().resolve() / "tutorial_square_plate"
RUN_DIR = base
suffix = 2
while RUN_DIR.exists():
    RUN_DIR = Path(f"{base}_{suffix}")
    suffix += 1
RUN_DIR.mkdir(parents=True)

GEO_SOURCE = Path("examples/dynamic/B3_dynamic_sent/mesh.geo")
BASE_CONFIG = Path("examples/dynamic/B3_dynamic_sent/config.yaml")
GEO_PATH = RUN_DIR / "mesh.geo"
MESH_PATH = RUN_DIR / "mesh.msh"
CONFIG_PATH = RUN_DIR / "config.yaml"

print(f"Repository: {ROOT.name}")
print(f"Tutorial output: {RUN_DIR}")
print(f"PhAST import: {Path(phast.__file__).resolve().relative_to(ROOT)}")

shutil.copy2(GEO_SOURCE, GEO_PATH)
geo_text = GEO_PATH.read_text(encoding="utf-8")
print(f"Loaded geometry: {GEO_PATH}")
print("\nGeometry entities and physical groups:\n")
for line in geo_text.splitlines():
    stripped = line.strip()
    if stripped.startswith(("Point(", "Line(", "Physical Curve", "Physical Surface")):
        print(stripped)

gmsh_executable = shutil.which("gmsh")
if gmsh_executable is None:
    raise RuntimeError(
        "Gmsh was not found. Reinstall the PhAST environment or follow "
        "the Gmsh recovery advice in docs/install.md."
    )

completed = subprocess.run(
    [gmsh_executable, str(GEO_PATH), "-2", "-format", "msh2", "-o", str(MESH_PATH)],
    cwd=ROOT,
    text=True,
    capture_output=True,
    timeout=60,
)
if completed.returncode != 0:
    raise RuntimeError(completed.stderr or completed.stdout)

mesh = FEMMesh(str(MESH_PATH), device="cpu")
nodes = require_finite(mesh.nodes.detach().cpu().numpy(), "mesh coordinates")
elements = mesh.elements.detach().cpu().numpy()
triangulation = mtri.Triangulation(nodes[:, 0], nodes[:, 1], elements)

boundary_colours = {
    "top": "#2f6f9f", "bottom": "#2f6f9f",
    "left": "#4f4f4f", "right": "#4f4f4f",
    "notch_upper": "#b24a3b", "notch_lower": "#b24a3b",
}
fig, ax = plt.subplots(figsize=(8.2, 7.0), constrained_layout=True)
ax.triplot(triangulation, color="#c8c8c8", linewidth=0.35)
for name, colour in boundary_colours.items():
    indices = mesh.node_sets[name].detach().cpu().numpy()
    points = nodes[indices]
    ax.scatter(points[:, 0], points[:, 1], s=13, color=colour, label=name, zorder=3)
ax.set(aspect="equal", xlabel="x [mm]", ylabel="y [mm]",
       title="Generated mesh and named boundary node sets")
ax.legend(ncol=3, loc="upper center", bbox_to_anchor=(0.5, -0.09), frameon=False)
plt.show()
print(mesh.summary())

source_config = yaml.safe_load(BASE_CONFIG.read_text(encoding="utf-8"))
m = source_config["material"]
material = Material(
    E=m["E"], nu=m["nu"], Gc=m["Gc"], l0=m["l0"], rho=m["rho"],
    eta_residual=m["eta_residual"], energy_split=m["energy_split"],
    pf_model=m["pf_model"], plane_stress=m.get("plane_stress", False),
)
h_min = float(mesh.h_min)
print(material)
print(f"h_min / l0 = {h_min / float(material.l0):.3f}")

display(Markdown(
    "~~~yaml\n"
    + yaml.safe_dump(
        {"boundary_conditions": source_config["boundary_conditions"],
         "loading": source_config["loading"]},
        sort_keys=False,
    )
    + "~~~"
))

NUM_STEPS = 6101
SNAPSHOT_EVERY = 50
PRINT_EVERY = 100
DT_SAFETY = 0.8

dt = DT_SAFETY * h_min / float(material.c_p)
t_total = NUM_STEPS * dt
assert (NUM_STEPS - 1) % SNAPSHOT_EVERY == 0

config = copy.deepcopy(source_config)
config["problem"]["name"] = "Square-plate SENT: executable teaching calculation"
config["geometry"] = {"mesh_path": "mesh.msh"}
config["loading"].update({"num_steps": NUM_STEPS, "dt": dt, "t_total": t_total})
config["solver"].update({
    "solver_type": "explicit",
    "time_integrator": "central_difference",
    "dt_safety": DT_SAFETY,
    "damage_every": 1,
    "use_multigrid": False,
    "preconditioner": "jacobi",
    "damage_max_iter": 1000,
})
config["output"].update({
    "trajectory": True,
    "h5": True,
    "trajectory_format": "h5",
    "h5_every": SNAPSHOT_EVERY,
    "print_every": PRINT_EVERY,
    "plots": False,
    "gif": False,
    "profile": True,
})
config["device"] = {"device": "cpu", "compile": False}
CONFIG_PATH.write_text(yaml.safe_dump(config, sort_keys=False), encoding="utf-8")

print(f"Configuration: {CONFIG_PATH}")
print("Selected route: dynamic / explicit mechanics / implicit AT2 damage")
print(f"Number of steps: {NUM_STEPS}")
print(f"Stable dt: {dt:.6e} s")
print(f"End time: {t_total * 1e6:.3f} microseconds")
print(f"Stored snapshots: {(NUM_STEPS - 1) // SNAPSHOT_EVERY + 1}")

preflight = subprocess.run(
    [sys.executable, "-m", "phast", "run", str(CONFIG_PATH), "--validate-only"],
    cwd=ROOT, text=True, capture_output=True, timeout=60,
)
print(preflight.stdout, end="")
print(preflight.stderr, end="")
if preflight.returncode != 0:
    raise RuntimeError("Configuration preflight failed.")

if any((RUN_DIR / name).exists()
       for name in ("result.vtkhdf", "training_data.h5", "training_data.zarr")):
    raise RuntimeError(
        f"{RUN_DIR} already contains a result. Restart the kernel to create "
        "a new numbered directory rather than overwriting evidence."
    )

run_command = [
    sys.executable, "-m", "phast", "run", str(CONFIG_PATH),
    "--device", "cpu", "--output_dir", str(RUN_DIR),
]
progress = display(
    HTML(f"<progress value='0' max='{NUM_STEPS}' style='width:70%'></progress> "
         f"0 / {NUM_STEPS} steps"),
    display_id=True,
)
process = subprocess.Popen(
    run_command, cwd=ROOT, stdout=subprocess.PIPE, stderr=subprocess.STDOUT,
    text=True, bufsize=1,
)
assert process.stdout is not None
for line in process.stdout:
    print(line, end="")
    match = re.search(r"Step\s+(\d+)", line)
    if match:
        completed_step = int(match.group(1)) + 1
        progress.update(
            HTML(f"<progress value='{completed_step}' max='{NUM_STEPS}' "
                 f"style='width:70%'></progress> "
                 f"{completed_step:,} / {NUM_STEPS:,} steps")
        )
return_code = process.wait()
if return_code != 0:
    raise RuntimeError(f"PhAST exited with status {return_code}.")
progress.update(
    HTML(f"<progress value='{NUM_STEPS}' max='{NUM_STEPS}' "
         f"style='width:70%'></progress> "
         f"{NUM_STEPS:,} / {NUM_STEPS:,} steps complete")
)

result = phast.load_result(RUN_DIR)
print("Stored fields:", result.field_names())

final_step = NUM_STEPS - 1
final_fields = {}
for name in ("displacement", "strain", "stress", "damage"):
    if not result.has_field(name):
        if name in {"displacement", "damage"}:
            raise RuntimeError(f"Required {name} trajectory is missing.")
        print(f"Optional {name} output was not stored; its plot will be omitted.")
        continue
    field = require_finite(result.field(name, step=-1), name)  # last stored step
    final_fields[name] = field
    print(f"{name:12s} shape={field.shape}, min={field.min():.6g}, max={field.max():.6g}")

displacement = final_fields["displacement"]
damage = final_fields["damage"]
strain = final_fields.get("strain")
stress = final_fields.get("stress")

# Stored steps are numbered 0, 1, 2, ...; the history maps solver steps to them.
stored_steps = np.asarray(result.steps())
stored_times = require_finite(result.times(), "stored times")
history = result.history_columns()
solver_steps = np.asarray(history["increment"])[np.asarray(history["step_index"]) >= 0]
if len(solver_steps) != len(stored_steps) or int(solver_steps[-1]) != final_step:
    raise RuntimeError("The final solver step is not the last stored step.")
print(f"{len(stored_steps)} stored steps; solver steps "
      f"{solver_steps[0]}, {solver_steps[1]}, ..., {solver_steps[-1]}")

mesh_arrays = result.mesh_arrays()
coordinates = require_finite(mesh_arrays["points"], "stored mesh coordinates")
connectivity = np.asarray(mesh_arrays["cells"]["T3"])
if displacement.shape != (len(coordinates), 2) or damage.shape != (len(coordinates),):
    raise ValueError("Stored nodal fields do not match the mesh.")
for name, field in (("strain", strain), ("stress", stress)):
    if field is not None and field.shape != (len(connectivity), 3):
        raise ValueError(f"Unexpected in-plane {name} shape: {field.shape}")

plot_mesh = mtri.Triangulation(coordinates[:, 0], coordinates[:, 1], connectivity)

u_magnitude = require_finite(np.linalg.norm(displacement, axis=1), "displacement magnitude")
epsilon_yy = None if strain is None else strain[:, 1]
in_plane_equivalent_stress = None
if stress is not None:
    sigma_xx, sigma_yy, sigma_xy = stress.T
    in_plane_equivalent_stress = require_finite(
        np.sqrt(np.maximum(
            sigma_xx**2 - sigma_xx * sigma_yy + sigma_yy**2 + 3.0 * sigma_xy**2,
            0.0,
        )),
        "in-plane equivalent stress",
    )

fig, axes = plt.subplots(2, 2, figsize=(12, 9), constrained_layout=True)
specs = [
    ("Displacement magnitude, |u| [mm]", u_magnitude, "cividis", "nodal"),
    ("Vertical strain, epsilon_yy [-]", epsilon_yy, "RdBu_r", "element"),
    ("In-plane equivalent stress [MPa]", in_plane_equivalent_stress, "viridis", "element"),
    ("Phase-field damage, d [-]", damage, "Reds", "nodal"),
]
for ax, (title, field, cmap, location) in zip(axes.flat, specs):
    if field is None:
        ax.text(0.5, 0.5, f"{title}\nOptional field not stored",
                ha="center", va="center", transform=ax.transAxes)
        ax.axis("off")
        continue
    field = require_finite(field, title)
    limits = dict(vmin=0.0, vmax=1.0) if cmap == "Reds" else {}
    artist = (
        ax.tripcolor(plot_mesh, field, shading="gouraud", cmap=cmap, **limits)
        if location == "nodal"
        else ax.tripcolor(plot_mesh, facecolors=field, shading="flat", cmap=cmap, **limits)
    )
    ax.set(aspect="equal", xlabel="x [mm]", ylabel="y [mm]", title=title)
    fig.colorbar(artist, ax=ax, shrink=0.82)
fig.suptitle(
    f"Square-plate SENT: step {final_step} at {stored_times[-1] * 1e6:.2f} microseconds",
    fontsize=14,
)
field_figure = RUN_DIR / "final_fields.png"
fig.savefig(field_figure, dpi=160)
plt.show()
print(f"Saved: {field_figure}")

frame_indices = np.unique(np.linspace(0, len(stored_steps) - 1, 62, dtype=int))
damage_frames = [
    require_finite(result.field("damage", step=int(stored_steps[index])),
                   f"damage at solver step {solver_steps[index]}")
    for index in frame_indices
]
if any(field.shape != (len(coordinates),) for field in damage_frames):
    raise ValueError("A damage frame does not match the mesh nodes.")
fig, ax = plt.subplots(figsize=(7.2, 6.2), constrained_layout=True)
artist = ax.tripcolor(
    plot_mesh, damage_frames[0],
    shading="gouraud", cmap="Reds", vmin=0.0, vmax=1.0,
)
ax.set(aspect="equal", xlabel="x [mm]", ylabel="y [mm]")
fig.colorbar(artist, ax=ax).set_label("damage d")

def update_animation(frame_number: int):
    index = frame_indices[frame_number]
    artist.set_array(damage_frames[frame_number])
    ax.set_title(
        f"Crack evolution: solver step {solver_steps[index]} | "
        f"t = {stored_times[index] * 1e6:.2f} microseconds"
    )
    return (artist,)

animation = FuncAnimation(
    fig, update_animation, frames=len(frame_indices), interval=100, blit=False
)
animation_path = RUN_DIR / "crack_evolution.gif"
animation.save(animation_path, writer=PillowWriter(fps=10), dpi=115)
plt.close(fig)
display(Image(filename=str(animation_path)))
print(f"Saved: {animation_path}")

high_damage_nodes = coordinates[damage >= 0.95]
if len(high_damage_nodes):
    print(
        "Nodes with d >= 0.95 span "
        f"x = {high_damage_nodes[:, 0].min():.3f} to "
        f"{high_damage_nodes[:, 0].max():.3f} mm."
    )
print("\nTop-level run artifacts:")
for artifact in sorted(RUN_DIR.iterdir()):
    if artifact.name not in {"result.vtkhdf", "training_data.h5"}:
        print(f"  {artifact.name}")
print("  result.vtkhdf     (result file: fields, history, mesh; opens in ParaView)")
print("  training_data.h5  (single-file trajectory store, requested above)")