Step-by-step PhAST problem setup¶

import os
import sys
import subprocess
from pathlib import Path

IN_COLAB = "google.colab" in sys.modules
PHAST_REPOSITORY = "https://github.com/CEMS-Lab/PhAST.git"
PHAST_REVISION = "main"  # Mutable public branch, not an immutable release.

if IN_COLAB:
    subprocess.run(["apt-get", "update", "-qq"], check=True)
    subprocess.run(["apt-get", "install", "-y", "-qq", "gmsh"], check=True)
    checkout = Path("/content/PhAST-main")
    if checkout.exists():
        if not (checkout / ".git").is_dir():
            raise RuntimeError(
                f"{checkout} is not a Git checkout. Restart the Colab runtime."
            )
        origin = subprocess.check_output(
            ["git", "-C", str(checkout), "config", "--get", "remote.origin.url"],
            text=True,
        ).strip()
        if origin != PHAST_REPOSITORY:
            raise RuntimeError(f"Unexpected cached checkout origin: {origin!r}")
    else:
        subprocess.run(
            ["git", "clone", "--branch", PHAST_REVISION, "--depth", "1",
             PHAST_REPOSITORY, str(checkout)],
            check=True,
        )
    subprocess.run(
        ["git", "-C", str(checkout), "fetch", "--depth", "1",
         "origin", "refs/heads/main"],
        check=True,
    )
    expected_commit = subprocess.check_output(
        ["git", "-C", str(checkout), "rev-parse", "FETCH_HEAD"], text=True
    ).strip()
    resolved_commit = subprocess.check_output(
        ["git", "-C", str(checkout), "rev-parse", "HEAD"], text=True
    ).strip()
    tracked_changes = subprocess.check_output(
        ["git", "-C", str(checkout), "status", "--porcelain", "--untracked-files=no"],
        text=True,
    ).strip()
    if resolved_commit != expected_commit or tracked_changes:
        raise RuntimeError(
            "The cached checkout differs from public main. Restart the Colab "
            "runtime rather than overwriting local changes."
        )
    print(f"Public source: {PHAST_REVISION} (mutable)")
    print(f"Resolved commit: {resolved_commit}")
    subprocess.run(
        [sys.executable, "-m", "pip", "install", "-q", "-e", str(checkout), "jupyterlab"],
        check=True,
    )
    os.chdir(checkout)
else:
    print("Local/docs execution: using the current checkout; no packages installed.")

import json
import shutil

from PIL import Image as PILImage
import matplotlib.pyplot as plt
import meshio
import numpy as np
import yaml

here = Path.cwd().resolve()
repo_root = next(
    (candidate for candidate in (here, *here.parents)
     if (candidate / "pyproject.toml").is_file()
     and (candidate / "src" / "phast").is_dir()),
    None,
)
if repo_root is None:
    raise RuntimeError("Start this notebook from the current PhAST source checkout.")
source_src = repo_root / "src"
sys.path.insert(0, str(source_src))
os.environ["PYTHONPATH"] = str(source_src) + os.pathsep + os.environ.get("PYTHONPATH", "")

import phast

if not Path(phast.__file__).resolve().is_relative_to(source_src):
    raise RuntimeError("Restart the kernel: PhAST was imported from a different checkout.")

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", repo_root / "runs"))
base = output_root.expanduser().resolve() / "notebook_sent"
run_dir = base
suffix = 2
while run_dir.exists():
    run_dir = base.with_name(f"{base.name}_{suffix}")
    suffix += 1
run_dir.mkdir(parents=True)

geo_path = run_dir / "mesh.geo"
mesh_path = run_dir / "mesh.msh"
config_path = run_dir / "config.yaml"
output_dir = run_dir / "results"

print("PhAST:", getattr(phast, "__version__", "source checkout"))
print("Repository root:", repo_root)
print("Working directory:", run_dir)

L = 1.0          # plate width [mm]
H = 1.0          # plate height [mm]
a0 = 0.50        # initial notch/crack length [mm]
y_crack = 0.50   # crack vertical location [mm]
notch_gap = 0.012 # small geometric opening used for robust meshing [mm]
h_bulk = 0.08    # coarse element size [mm]
h_crack = 0.02   # local element size near the crack [mm]

geo_text = f"""
SetFactory("OpenCASCADE");

L = {L};
H = {H};
a0 = {a0};
yc = {y_crack};
g = {notch_gap};
h_bulk = {h_bulk};
h_crack = {h_crack};

Point(1) = {{0, 0, 0, h_bulk}};
Point(2) = {{L, 0, 0, h_bulk}};
Point(3) = {{L, H, 0, h_bulk}};
Point(4) = {{0, H, 0, h_bulk}};
Point(5) = {{0, yc + 0.5*g, 0, h_crack}};
Point(6) = {{a0, yc + 0.5*g, 0, h_crack}};
Point(7) = {{a0, yc - 0.5*g, 0, h_crack}};
Point(8) = {{0, yc - 0.5*g, 0, h_crack}};

Line(1) = {{1, 2}};
Line(2) = {{2, 3}};
Line(3) = {{3, 4}};
Line(4) = {{4, 5}};
Line(5) = {{5, 6}};
Line(6) = {{6, 7}};
Line(7) = {{7, 8}};
Line(8) = {{8, 1}};

Curve Loop(10) = {{1, 2, 3, 4, 5, 6, 7, 8}};
Plane Surface(20) = {{10}};

Field[1] = Distance;
Field[1].CurvesList = {{5, 6, 7}};
Field[2] = Threshold;
Field[2].InField = 1;
Field[2].SizeMin = h_crack;
Field[2].SizeMax = h_bulk;
Field[2].DistMin = 0.02;
Field[2].DistMax = 0.18;
Background Field = 2;

Physical Surface("body") = {{20}};
Physical Curve("bottom") = {{1}};
Physical Curve("right") = {{2}};
Physical Curve("top") = {{3}};
Physical Curve("left") = {{4, 8}};
Physical Curve("notch") = {{5, 6, 7}};
""".strip()

geo_path.write_text(geo_text + "\n", encoding="utf-8")
print(geo_path)
print(geo_text[:600] + "\n...")

fig, ax = plt.subplots(figsize=(6, 5))
ax.plot([0, L, L, 0, 0], [0, 0, H, H, 0], color="0.15", lw=2)
ax.fill([0, a0, a0, 0], [y_crack - notch_gap/2, y_crack - notch_gap/2, y_crack + notch_gap/2, y_crack + notch_gap/2], color="crimson", alpha=0.22)
ax.plot([0, a0], [y_crack, y_crack], color="crimson", lw=3, label="initial notch / damage seed")
ax.plot([0, L], [0, 0], color="royalblue", lw=4, label="fixed bottom")
ax.plot([0, L], [H, H], color="darkorange", lw=4, label="prescribed top displacement")
ax.set_aspect("equal")
ax.set_xlabel("x [mm]")
ax.set_ylabel("y [mm]")
ax.set_title("Geometry and named boundaries")
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.12), ncol=1)
fig.tight_layout()
fig.savefig(run_dir / "geometry_preview.png", dpi=180)
plt.show()

gmsh = shutil.which("gmsh")
if gmsh is None:
    raise RuntimeError("gmsh executable was not found on PATH")

cmd = [gmsh, "-2", str(geo_path), "-format", "msh2", "-o", str(mesh_path)]
print(" ".join(cmd))
subprocess.run(cmd, check=True, timeout=60)
print(mesh_path, mesh_path.stat().st_size, "bytes")

mesh_summary = phast.inspect_mesh(mesh_path)
print("points:", mesh_summary["n_points"])
print("cell blocks:", mesh_summary["cells"])
print(json.dumps(mesh_summary["named_groups"], indent=2))

required_groups = {"body", "bottom", "top", "notch"}
missing = required_groups - set(mesh_summary["named_groups"])
if missing:
    raise RuntimeError(f"Missing required physical groups: {sorted(missing)}")

mesh = meshio.read(mesh_path)
points = require_finite(mesh.points[:, :2], "mesh coordinates")

triangles = []
lines = []
line_tags = []
for block_index, block in enumerate(mesh.cells):
    if block.type == "triangle":
        triangles.append(block.data)
    if block.type == "line":
        lines.append(block.data)
        tags = mesh.cell_data.get("gmsh:physical", [])[block_index]
        line_tags.append(tags)

triangles = np.vstack(triangles) if triangles else np.empty((0, 3), dtype=int)
lines = np.vstack(lines) if lines else np.empty((0, 2), dtype=int)
line_tags = np.concatenate(line_tags) if line_tags else np.empty((0,), dtype=int)

id_to_name = {int(data[0]): name for name, data in mesh.field_data.items()}
colors = {"bottom": "royalblue", "top": "darkorange", "left": "0.35", "right": "0.35", "notch": "crimson"}

fig, ax = plt.subplots(figsize=(7, 6))
if len(triangles):
    ax.triplot(points[:, 0], points[:, 1], triangles, color="0.82", linewidth=0.5)
for edge, tag in zip(lines, line_tags):
    name = id_to_name.get(int(tag), str(tag))
    xy = points[edge]
    ax.plot(xy[:, 0], xy[:, 1], color=colors.get(name, "0.4"), lw=2.0 if name == "notch" else 1.4)
ax.set_aspect("equal")
ax.set_xlabel("x [mm]")
ax.set_ylabel("y [mm]")
ax.set_title("Mesh with named physical curves")
fig.tight_layout()
fig.savefig(run_dir / "mesh_preview.png", dpi=180)
plt.show()

num_steps = 2          # quick tutorial run; increase for smoother damage evolution
u_top = 0.0002         # prescribed final displacement [mm]
l0 = 0.04             # phase-field length scale [mm]
enable_damage_preseed = False

problem = (
    phast.Problem("Notebook SENT tutorial")
    .mesh(str(mesh_path))
    .region("body", kind="domain", from_mesh="body")
    .region("bottom", from_mesh="bottom")
    .region("top", from_mesh="top")
    .region("notch", from_mesh="notch")
    .material(
        "tutorial_material",
        region="body",
        E=210000.0,
        nu=0.3,
        Gc=2.7,
        l0=l0,
        rho=7.8e-09,
        eta_residual=1.0e-07,
        energy_split="spectral",
        pf_model="AT2",
        plane_stress=False,
    )
    .boundary_condition("fix", region="bottom", dof="x", name="fix_bottom_x")
    .boundary_condition("fix", region="bottom", dof="y", name="fix_bottom_y")
    .boundary_condition("displacement", region="top", dof="y", value=1.0, name="pull_top")
    .analysis_step(
        "load",
        kind="quasi_static",
        controls={"protocol": "simple", "num_steps": num_steps, "dt": u_top / num_steps},
        active_boundary_conditions=["fix_bottom_x", "fix_bottom_y", "pull_top"],
    )
    .solver(
        "quasi_static",
        stagger_tol=1.0e-6,
        max_stagger=50,
        preconditioner="jacobi",
        damage_tol=1.0e-5,
        static_tol=1.0e-7,
        damage_max_iter=500,
        static_max_iter=500,
        fail_on_mechanics_nonconvergence=True,
        fail_on_stagger_nonconvergence=True,
        backend="auto",
        device="cpu",
    )
    .outputs(
        fields=[{"name": "trajectory", "every": 1, "format": "h5"}],
        histories=[{"name": "reaction_force", "region": "bottom", "dof": "y"}],
        output_dir=str(output_dir),
        plots=False,  # Plot only after checking the stored numerical fields.
        profile=True,
        gif=False,
        print_every=1,
    )
)

if enable_damage_preseed:
    problem.initial_condition("damage", region="notch", value=1.0)

problem.validate_setup()

setup_preview = run_dir / "initial_conditions.png"
problem.plot_setup(output=setup_preview)

img = np.asarray(PILImage.open(setup_preview).convert("RGB"))
fig, ax = plt.subplots(figsize=(8, 5))
ax.imshow(img)
ax.axis("off")
ax.set_title("PhAST setup preview")
plt.show()
print(setup_preview)

problem.save(config_path)
saved_config = yaml.safe_load(config_path.read_text(encoding="utf-8"))
if saved_config.get("schema_version") != 2:
    raise RuntimeError("This workflow requires the saved schema-v2 configuration.")
print("Saved configuration used for validation and execution:", config_path)
print(yaml.safe_dump(saved_config, sort_keys=False))

validate_cmd = [sys.executable, "-m", "phast", "run", str(config_path), "--validate-only"]
print(" ".join(validate_cmd))
subprocess.run(validate_cmd, check=True, cwd=repo_root, timeout=60)

run_cmd = [
    sys.executable, "-m", "phast", "run", str(config_path),
    "--output_dir", str(output_dir),
]
print(" ".join(run_cmd))
subprocess.run(run_cmd, check=True, cwd=repo_root, timeout=120)

# Inspect the output contract without assuming every optional writer was enabled.
artifact_groups = {
    "provenance": ["run_manifest.json", "run_metadata.json", "run_lockfile.json"],
    "history": ["history.csv", "results.csv", "reaction_force.csv"],
    "field_or_trajectory": [
        "damage_final.png",
        "displacement_final.png",
        "result.vtkhdf",
        "training_data.h5",
        "training_data.zarr",
    ],
}

print("Artifact contract check")
for category, candidates in artifact_groups.items():
    present = [name for name in candidates if (output_dir / name).exists()]
    print(f"{category:>19}: {present if present else 'not found under standard names'}")

print()
print(
    "Interpretation: availability depends on the selected workflow and output "
    "settings. Inspect the run manifest before treating an absent optional artifact "
    "as a solver failure."
)

for path in sorted(output_dir.iterdir()):
    print(path.name + ("/" if path.is_dir() else ""))

result = phast.load_result(output_dir)
print("\nmetadata keys:", sorted(result.metadata().keys()))
print("histories:", result.history_names())
print("fields:", result.field_names())
print("visuals:", result.visuals())

stored_steps = []
result_file = output_dir / "result.vtkhdf"
if result_file.is_file():
    # The result file numbers its stored steps 0, 1, 2, ...; history/step_index
    # maps each solver step to its stored step (-1 when no fields were stored).
    stored_steps = [int(step) for step in result.steps()]
    history = result.history_columns()
    solver_steps = np.asarray(history["increment"])[np.asarray(history["step_index"]) >= 0]
    if not stored_steps or len(solver_steps) != len(stored_steps) or int(solver_steps[-1]) != num_steps - 1:
        raise RuntimeError("The requested final step is missing from the result file.")
    print("Stored steps:", stored_steps, "| solver steps:", solver_steps.tolist())
else:
    print("No result file is available; field plots and animation will be omitted.")

if "load_displacement" in result.history_names():
    rows = result.history("load_displacement")
    if rows:
        columns = set(rows[0])
        if {"applied_disp", "reaction_force"} <= columns:
            x_key, y_key, force_unit = "applied_disp", "reaction_force", "N"
        elif {"displacement", "reaction_kN"} <= columns:
            x_key, y_key, force_unit = "displacement", "reaction_kN", "kN"
        else:
            raise ValueError(f"Unrecognised load-displacement columns: {sorted(columns)}")
        applied = require_finite(
            np.asarray([row[x_key] for row in rows], dtype=float), x_key
        )
        reaction = require_finite(
            np.asarray([row[y_key] for row in rows], dtype=float), y_key
        )
        fig, ax = plt.subplots(figsize=(6, 4))
        ax.plot(applied, reaction, "o-")
        ax.set(xlabel="Applied displacement [mm]", ylabel=f"Reaction [{force_unit}]",
               title="Two-step response, not a converged fracture curve")
        fig.tight_layout()
        plt.show()
    else:
        print("The load-displacement history is empty; no curve is drawn.")
else:
    print("No load-displacement history was stored; no curve is fabricated.")


def write_scalar_field_png(values, title, path, cmap="viridis", limits=None):
    arr = require_finite(values, title)
    if arr.ndim != 1 or arr.size not in (points.shape[0], triangles.shape[0]):
        raise ValueError(f"{title}: field shape {arr.shape} does not match this mesh")
    colour_limits = {} if limits is None else dict(vmin=limits[0], vmax=limits[1])
    fig, ax = plt.subplots(figsize=(6, 5))
    if arr.size == points.shape[0]:
        im = ax.tripcolor(
            points[:, 0], points[:, 1], triangles, arr,
            shading="gouraud", cmap=cmap, **colour_limits,
        )
    else:
        im = ax.tripcolor(
            points[:, 0], points[:, 1], triangles, facecolors=arr,
            shading="flat", cmap=cmap, **colour_limits,
        )
    ax.set(aspect="equal", title=title, xlabel="x [mm]", ylabel="y [mm]")
    fig.colorbar(im, ax=ax, shrink=0.82)
    fig.tight_layout()
    fig.savefig(path, dpi=180)
    plt.close(fig)
    return path


checked_images = []
for field_name, title, cmap, component, limits in [
    ("damage", "Final phase-field damage, d [-]", "Reds", None, (0.0, 1.0)),
    ("displacement", "Final displacement magnitude [mm]", "cividis", "magnitude", None),
    ("stress", "Final in-plane stress, sigma_yy [MPa]", "RdBu_r", 1, None),
    ("strain", "Final vertical strain, epsilon_yy [-]", "RdBu_r", 1, None),
]:
    if not stored_steps or not result.has_field(field_name):
        print(f"No stored {field_name} field; its plot is omitted.")
        continue
    values = require_finite(result.field(field_name, step=stored_steps[-1]), field_name)
    if component == "magnitude":
        if values.ndim != 2 or values.shape[1] != 2:
            raise ValueError(f"Unexpected displacement shape: {values.shape}")
        values = np.linalg.norm(values, axis=1)
    elif component is not None:
        if values.ndim != 2 or values.shape[1] != 3:
            raise ValueError(f"Unexpected in-plane {field_name} shape: {values.shape}")
        values = values[:, component]
    checked_images.append(write_scalar_field_png(
        values, title, output_dir / f"{field_name}_final.png", cmap=cmap, limits=limits
    ))

if checked_images:
    fig, axes = plt.subplots(
        (len(checked_images) + 1) // 2, 2,
        figsize=(11, 5 * ((len(checked_images) + 1) // 2)),
    )
    for ax in np.asarray(axes).flat:
        ax.axis("off")
    for ax, path in zip(np.asarray(axes).flat, checked_images):
        ax.imshow(np.asarray(PILImage.open(path).convert("RGB")))
        ax.set_title(path.name)
    fig.tight_layout()
    plt.show()
else:
    print("No available finite trajectory fields to plot.")

from IPython.display import Image, display

gif_path = output_dir / "damage_evolution.gif"
if stored_steps and result.has_field("damage"):
    damage_frames = [
        require_finite(result.field("damage", step=step), f"damage at step {step}")
        for step in stored_steps
    ]
    if any(field.shape != (len(points),) for field in damage_frames):
        raise ValueError("A damage snapshot does not match the mesh nodes.")
    frames = []
    for step_index, damage in zip(stored_steps, damage_frames):
        fig, ax = plt.subplots(figsize=(6, 5))
        im = ax.tripcolor(
            points[:, 0], points[:, 1], triangles, damage,
            shading="gouraud", cmap="Reds", vmin=0.0, vmax=1.0,
        )
        ax.set(aspect="equal", title=f"Damage, stored step {step_index}",
               xlabel="x [mm]", ylabel="y [mm]")
        fig.colorbar(im, ax=ax, shrink=0.82)
        fig.tight_layout()
        fig.canvas.draw()
        frames.append(PILImage.fromarray(
            np.asarray(fig.canvas.buffer_rgba())[:, :, :3].copy()
        ))
        plt.close(fig)
    frames[0].save(
        gif_path, save_all=True, append_images=frames[1:], duration=800, loop=0
    )
    display(Image(filename=str(gif_path)))
else:
    print("No stored damage trajectory; no animation is fabricated.")