#!/usr/bin/env python3
"""Tier-2 item 10 -- is the h/p convergence study measuring one problem?

Why this diagnostic matters
---------------------------
The chromatin input is a 50 x 50 x 50 voxel volume mapped onto the core, which
is about 0.30 x 0.22 x 0.14 um per voxel. The finest mesh has about 7531
vertices, i.e. about 0.43 um nodal spacing -- undersampled by roughly an order
of magnitude -- and the ten grayscale groups alias badly at that spacing.

The assignment is per tetrahedron. The solver records it as
`"assignment": "reference_mesh_tetrahedron_centroid_trilinear_sampling"`
(nonlinear_solver.py, `core_heterogeneity_metadata`) and `_build_material_regions`
builds one material region per grayscale group out of `group_ids`, one id per
tetrahedron. A centroid sample is a point sample, not a cell average, so
refining the mesh RESAMPLES the material field: h-refinement is not converging
to a fixed continuous problem, it is walking across a family of problems whose
material distributions differ.

Direct evidence in the submitted numbers: the P1 family Richardson-extrapolates
to a Control endpoint of 0.066378 mN/m while intermediate-mesh P2 already gives
0.067416 mN/m -- a 1.54 % gap between two discretisations of what should be the
same continuous problem.

The published summary CSV contains the smoking gun on its own, because it
reports the realised phase volume fractions per discretisation:

    coarse_h1p25        0.507108 / 0.354038 / 0.138852
    intermediate_h0p75  0.569446 / 0.310014 / 0.120538
    main_h0p45          0.569430 / 0.310033 / 0.120536

P1 and P2 on the SAME mesh are bit-identical -- p-refinement cannot move the
material field, because the groups are element sets. Different meshes are not.
The pipeline's own acceptance rule (`_evaluate_acceptance`) requires
`max |phase - [0.57, 0.31, 0.12]| <= 0.02` and applies it only to the primary
mesh; the coarse mesh, which is the base point of the Richardson extrapolation,
misses that band by about 3x.

What this script does
---------------------
Three layers, each usable on its own:

  A. published summary (CSV only, no mesh, no solve)
     phase-fraction drift across meshes, the p-refinement null control, the
     project's own phase-fraction acceptance rule applied to EVERY row, and the
     Richardson-vs-P2 gap recomputed from the same table.

  B. material resampling (needs the real tree, meshio and the chromatin .npz)
     calls the project's own `load_reynolds_cellwise_heterogeneity` on each real
     mesh and reports per-group volume fractions, Voigt and Reuss effective
     shear moduli, and their drift across meshes. If the material field were
     mesh-independent these would be identical; the drift IS the contamination.

  C. quadrature sweep (needs solves; --quadrature-sweep)
     re-runs the two-group solve at FIXED mesh and FIXED element order with
     `quadrature_order` raised 2 -> 4 -> 6 via `_common.patch_options`.
     Read it as a null control, and read it carefully: the material field is
     piecewise constant per tetrahedron, so raising the quadrature order cannot
     resample it. Endpoint movement here is therefore integration error at
     frozen material and frozen geometry, and it bounds how much of the reported
     h-drift may be plain quadrature error rather than material resampling.
     Note the coupling: the outer-surface post-processing uses
     `max(4, quadrature_order + 1)`, so part of any movement is post-processing,
     not equilibrium.

Usage:
    python3 10_material_projection_convergence.py --project-root /path/to/tree \
        [--summary-csv RUN_DIR_OR_CSV] [--json-out out.json] [--markdown-out out.md]

    python3 10_material_projection_convergence.py --project-root /path/to/tree \
        --quadrature-sweep --out-root /path/to/out/quadrature --dry-run
"""

from __future__ import annotations

import argparse
import importlib
import json
import math
import sys
from pathlib import Path

import numpy as np

sys.path.insert(0, str(Path(__file__).resolve().parent.parent / "tier1_scripts"))
from _common import (  # noqa: E402
    SUMMARY_NAME,
    discretizations,
    import_runner,
    patch_options,
    read_endpoints,
)

DEFAULT_MESHES = ("coarse_h1p25", "intermediate_h0p75", "main_h0p45")

# Representative element size in um, used for the refinement ratio only.
MESH_SIZE_UM = {"coarse_h1p25": 1.25, "intermediate_h0p75": 0.75, "main_h0p45": 0.45}

# _evaluate_acceptance in the runner: max |phase - target| <= 0.02, checked for
# the primary mesh only. This script applies the same rule to every row.
PHASE_TARGET = (0.57, 0.31, 0.12)
PHASE_TOLERANCE = 0.02
PHASE_COLUMNS = ("phase_low_GSG1_4", "phase_intermediate_GSG5_6", "phase_high_GSG7_10")
PHASE_GROUP_SLICES = {"low_GSG1_4": (0, 4), "intermediate_GSG5_6": (4, 6), "high_GSG7_10": (6, 10)}

ENDPOINT_COLUMN = "mean_shell_resultant_mn_per_m"


def _lazy_import(name: str, purpose: str):
    """Import a heavy optional dependency, or explain what is missing and why."""
    try:
        return importlib.import_module(name)
    except ImportError as exc:
        raise SystemExit(
            f"{name} is required for {purpose}, and it is not importable here.\n"
            f"Run this script in the environment that runs the FEM project "
            f"(the one with meshio, scipy and scikit-fem), or drop the option "
            f"that needs it: layer A works from the summary CSV alone."
        ) from exc


def _import_project(project_root: Path):
    """Return (load_reynolds_cellwise_heterogeneity, tags) from the real tree."""
    source_dir = project_root / "05_Source_Data"
    if not (source_dir / "nonlinear_solver.py").exists():
        source_dir = project_root
    for path in (str(source_dir), str(project_root / "src")):
        if path not in sys.path:
            sys.path.insert(0, path)
    try:
        chromatin = importlib.import_module("nuclear_envelope_fem.chromatin_heterogeneity")
        tags = importlib.import_module("nuclear_envelope_fem.tags")
    except ImportError as exc:
        raise SystemExit(
            "Could not import nuclear_envelope_fem from "
            f"{project_root}.\nThe submission zip does not ship "
            "src/nuclear_envelope_fem/, the gmsh meshes or the chromatin .npz; "
            "point --project-root at the real working tree."
        ) from exc
    return chromatin.load_reynolds_cellwise_heterogeneity, tags


def _runner(project_root: Path):
    """Import the upstream two-group runner, or say what the tree is missing."""
    try:
        return import_runner(project_root)
    except ImportError as exc:
        raise SystemExit(
            f"Importing the two-group runner from {project_root} failed: {exc}.\n"
            "It needs src/nuclear_envelope_fem, run_extreme_finite_deformation.py, "
            "run_vaziri_finan_two_group.py and validate_literature_reproduction.py, "
            "none of which ship in the submission zip. Point --project-root at the "
            "real working tree."
        ) from exc


def _find_mesh(project_root: Path, mesh_name: str) -> Path | None:
    candidates = sorted(project_root.rglob("*.msh"))
    exact = [path for path in candidates if mesh_name in str(path)]
    if exact:
        return exact[0]
    return None


def _resolve_meshes(project_root: Path, names, explicit) -> dict[str, Path]:
    """Map mesh label -> .msh path, preferring the runner's own MESHES table."""
    resolved: dict[str, Path] = {}
    for entry in explicit or ():
        label, _, path = str(entry).partition("=")
        if not path:
            path, label = label, Path(label).parent.name
        resolved[label] = Path(path)
    if resolved:
        return resolved
    try:
        runner = import_runner(project_root)
        table = getattr(runner, "MESHES", {})
    except (SystemExit, ImportError, OSError):
        table = {}
    for name in names:
        path = table.get(name)
        if path is not None and Path(path).exists():
            resolved[name] = Path(path)
            continue
        found = _find_mesh(project_root, name)
        if found is not None:
            resolved[name] = found
    return resolved


def _grid_shape(volume_path: Path) -> tuple[int, ...] | None:
    """Read the chromatin grid shape with numpy alone; None if unreadable."""
    try:
        with np.load(volume_path, allow_pickle=False) as handle:
            for key in handle.files:
                array = handle[key]
                if array.ndim == 3:
                    return tuple(int(size) for size in array.shape)
    except (OSError, ValueError, EOFError):
        return None
    return None


def _tetra_volumes(points: np.ndarray, cells: np.ndarray) -> np.ndarray:
    tetra = points[cells]
    return (
        np.abs(
            np.einsum(
                "ij,ij->i",
                tetra[:, 1] - tetra[:, 0],
                np.cross(tetra[:, 2] - tetra[:, 0], tetra[:, 3] - tetra[:, 0]),
            )
        )
        / 6.0
    )


def _group_shear_moduli(heterogeneity) -> np.ndarray:
    moduli = []
    for material in heterogeneity.materials_by_group:
        modulus = getattr(material, "shear_modulus_pa", None)
        if modulus is None:
            young = getattr(material, "young_modulus_pa", None)
            poisson = getattr(material, "poisson_ratio", None)
            if young is None or poisson is None:
                raise SystemExit(
                    "Cannot read a shear modulus off the material objects; "
                    "inspect materials_by_group and adjust this script."
                )
            modulus = young / (2.0 * (1.0 + poisson))
        moduli.append(float(modulus))
    return np.asarray(moduli, dtype=float)


def _richardson(series: list[tuple[float, float]]) -> dict | None:
    """Richardson estimate from three (h, f) pairs; refuses non-monotone data."""
    if len(series) < 3:
        return None
    series = sorted(series, key=lambda item: -item[0])
    (h1, f1), (h2, f2), (h3, f3) = series[-3:]
    first, second = f2 - f1, f3 - f2
    if second == 0.0 or first / second <= 0.0:
        return {"valid": False, "reason": "non-monotone increments; GCI undefined"}
    ratio = h1 / h2
    if not math.isclose(ratio, h2 / h3, rel_tol=0.05):
        return {"valid": False, "reason": "non-uniform refinement ratio"}
    order = math.log(first / second) / math.log(ratio)
    scale = ratio**order
    return {
        "valid": True,
        "observed_order": order,
        "extrapolated": f3 + second / (scale - 1.0),
        "finest": f3,
    }


def _summary_directory(candidate: Path) -> Path:
    if candidate.is_file():
        if candidate.name != SUMMARY_NAME:
            raise SystemExit(
                f"--summary-csv must point at {SUMMARY_NAME} or at the run "
                f"directory that contains it, got {candidate}"
            )
        return candidate.parent
    return candidate


def _discover_summary(project_root: Path) -> Path | None:
    matches = sorted(
        project_root.rglob(SUMMARY_NAME), key=lambda path: path.stat().st_mtime, reverse=True
    )
    return matches[0].parent if matches else None


def _mesh_of(discretization: str) -> str:
    return discretization[:-3] if discretization[-3:] in ("_P1", "_P2") else discretization


def _order_of(discretization: str) -> str:
    return discretization[-2:] if discretization[-3:] in ("_P1", "_P2") else "?"


def main() -> int:  # noqa: C901 - one linear report, split would hide the flow
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument(
        "--project-root",
        type=Path,
        required=True,
        help="Working tree containing src/nuclear_envelope_fem, meshes and results.",
    )
    parser.add_argument("--mesh-name", nargs="+", default=list(DEFAULT_MESHES))
    parser.add_argument(
        "--mesh-path",
        nargs="+",
        default=None,
        help="Explicit meshes as LABEL=PATH (or bare paths; label taken from the parent).",
    )
    parser.add_argument(
        "--template-volume",
        type=Path,
        default=None,
        help="Chromatin .npz (default: search the project tree).",
    )
    parser.add_argument(
        "--summary-csv",
        type=Path,
        default=None,
        help=f"Run directory or {SUMMARY_NAME} (default: newest under --project-root).",
    )
    parser.add_argument("--outer-axes-um", type=float, nargs=3, default=(7.5, 5.5, 3.5))
    parser.add_argument("--shell-thickness-um", type=float, default=0.12)
    parser.add_argument("--group-count", type=int, default=10)
    parser.add_argument("--minimum-shear-pa", type=float, default=170.0)
    parser.add_argument("--maximum-shear-pa", type=float, default=4000.0)
    parser.add_argument("--sigmoid-midpoint", type=float, default=0.5)
    parser.add_argument("--sigmoid-width", type=float, default=0.083)
    parser.add_argument("--bulk-to-shear-ratio", type=float, default=30.0)
    parser.add_argument(
        "--quadrature-sweep",
        action="store_true",
        help="Layer C: re-solve at fixed mesh with quadrature_order raised.",
    )
    parser.add_argument("--quadrature-orders", type=int, nargs="+", default=(2, 4, 6))
    parser.add_argument(
        "--out-root", type=Path, default=None, help="Required with --quadrature-sweep."
    )
    parser.add_argument("--primary-mesh", default="intermediate_h0p75")
    parser.add_argument("--element-order", choices=("p1", "p2"), default="p1")
    parser.add_argument("--sparse-backend", choices=("cpu", "cuda", "auto"), default="cpu")
    parser.add_argument(
        "--rigid-constraint-linear-solver",
        choices=("augmented_direct", "projected_minres"),
        default="augmented_direct",
    )
    parser.add_argument("--resume", action="store_true")
    parser.add_argument("--json-out", type=Path, default=None)
    parser.add_argument("--markdown-out", type=Path, default=None)
    parser.add_argument(
        "--dry-run",
        action="store_true",
        help="Resolve inputs, run the no-solve layers, print the planned solves and stop.",
    )
    args = parser.parse_args()
    if args.quadrature_sweep and args.out_root is None and not args.dry_run:
        raise SystemExit("--quadrature-sweep needs --out-root")

    lines: list[str] = []

    def emit(text: str = "") -> None:
        print(text)
        lines.append(text)

    report: dict = {"project_root": str(args.project_root)}

    meshes = _resolve_meshes(args.project_root, args.mesh_name, args.mesh_path)
    volume_path = args.template_volume
    if volume_path is None:
        matches = sorted(args.project_root.rglob("*chromatin_50cube.npz"))
        volume_path = matches[0] if matches else None
    summary_dir = (
        _summary_directory(args.summary_csv)
        if args.summary_csv is not None
        else _discover_summary(args.project_root)
    )

    emit("# Material projection and the h/p convergence study")
    emit()
    emit("| input | resolved |")
    emit("|---|---|")
    emit(f"| project root | `{args.project_root}` |")
    emit(f"| summary CSV | `{summary_dir}` |" if summary_dir else "| summary CSV | NOT FOUND |")
    emit(
        f"| chromatin volume | `{volume_path}` |"
        if volume_path
        else "| chromatin volume | NOT FOUND |"
    )
    for name in args.mesh_name:
        emit(f"| mesh {name} | " + (f"`{meshes[name]}` |" if name in meshes else "NOT FOUND |"))
    report["inputs"] = {
        "summary_dir": str(summary_dir) if summary_dir else None,
        "template_volume": str(volume_path) if volume_path else None,
        "meshes": {name: str(path) for name, path in meshes.items()},
    }

    inner_axes = np.asarray(args.outer_axes_um, dtype=float) - float(args.shell_thickness_um)
    grid = _grid_shape(volume_path) if volume_path else None
    if grid is None:
        grid = (50, 50, 50)
        grid_source = "assumed 50^3 (the volume was not readable here)"
    else:
        grid_source = f"read from {Path(volume_path).name}"
    # The loader is handed inner_axes and the voxel grid, so the pitch is taken as
    # the core bounding box over the grid. chromatin_heterogeneity.py does not ship
    # in the zip; this mapping is what reproduces the stated 0.30 x 0.22 x 0.14 um.
    pitch = 2.0 * inner_axes / np.asarray(grid, dtype=float)
    pitch_geometric_mean = float(np.prod(pitch) ** (1.0 / 3.0))
    voxel_volume = float(np.prod(pitch))
    body_volume = float(4.0 / 3.0 * np.pi * np.prod(np.asarray(args.outer_axes_um, dtype=float)))

    emit()
    emit("## 1. Voxel pitch against element size")
    emit()
    emit(f"Chromatin grid {grid[0]}x{grid[1]}x{grid[2]} ({grid_source}) over the core "
         f"bounding box 2 x ({inner_axes[0]:.2f}, {inner_axes[1]:.2f}, {inner_axes[2]:.2f}) um.")
    emit(f"Voxel pitch {pitch[0]:.4f} x {pitch[1]:.4f} x {pitch[2]:.4f} um "
         f"(geometric mean {pitch_geometric_mean:.4f} um, volume {voxel_volume:.6f} um^3).")
    emit()
    emit("| mesh | vertices | nodal spacing (um) | spacing / voxel pitch | voxels per nodal cell "
         "| mean tet volume (um^3) | voxels per tetrahedron |")
    emit("|---|---:|---:|---:|---:|---:|---:|")
    report["voxel_geometry"] = {
        "grid": list(grid),
        "grid_source": grid_source,
        "voxel_pitch_um": pitch.tolist(),
        "voxel_volume_um3": voxel_volume,
        "body_volume_um3": body_volume,
        "meshes": {},
    }

    summary_rows = None
    if summary_dir is not None:
        try:
            summary_rows = read_endpoints(summary_dir)
        except FileNotFoundError as exc:
            emit(f"warning: {exc}")

    vertex_counts: dict[str, int] = {}
    if summary_rows:
        for discretization, _ in summary_rows:
            if _order_of(discretization) != "P1":
                continue
            row = summary_rows[(discretization, "Control")]
            vertex_counts[_mesh_of(discretization)] = int(row["dofs"]) // 3

    mesh_metrics: dict[str, dict] = {}
    for name in args.mesh_name:
        vertices = vertex_counts.get(name)
        mean_tet = None
        if vertices is None and name not in meshes:
            continue
        if vertices is None:
            vertices = 0
        spacing = (body_volume / vertices) ** (1.0 / 3.0) if vertices else float("nan")
        ratio = spacing / pitch_geometric_mean
        metrics = {
            "vertices": vertices,
            "vertices_source": "summary CSV P1 dofs / 3" if vertices else "unknown",
            "nodal_spacing_um": spacing,
            "spacing_over_voxel_pitch": ratio,
            "voxels_per_nodal_cell": ratio**3,
            "mean_tetrahedron_volume_um3": mean_tet,
            "voxels_per_tetrahedron": None,
        }
        mesh_metrics[name] = metrics
        emit(
            f"| {name} | {vertices or '-'} | {spacing:.4f} | "
            f"{metrics['spacing_over_voxel_pitch']:.2f} | "
            f"{metrics['voxels_per_nodal_cell']:.1f} | - | - |"
        )
    report["voxel_geometry"]["meshes"] = mesh_metrics
    if mesh_metrics:
        emit()
        emit("Layer B fills the last two columns from the real meshes. A nodal cell wider "
             "than the voxel pitch cannot represent the field it samples: one grayscale "
             "sample is taken per tetrahedron centroid, and it stands for the whole cell.")

    emit()
    emit("## 2. Published summary: the material field is not mesh-independent")
    emit()
    if not summary_rows:
        emit("No summary CSV found, so layer A is unavailable. Pass --summary-csv.")
    else:
        emit("| discretisation | dofs | low GSG1-4 | intermediate GSG5-6 | high GSG7-10 "
             "| max abs(phase - target) | passes the project's own 0.02 rule |")
        emit("|---|---:|---:|---:|---:|---:|:--:|")
        phases: dict[str, tuple[float, float, float]] = {}
        acceptance: dict[str, dict] = {}
        for discretization in discretizations(summary_rows):
            row = summary_rows[(discretization, "Control")]
            values = tuple(float(row[column]) for column in PHASE_COLUMNS)
            error = float(np.max(np.abs(np.asarray(values) - np.asarray(PHASE_TARGET))))
            passed = error <= PHASE_TOLERANCE
            phases[discretization] = values
            acceptance[discretization] = {"phase_fractions": list(values),
                                          "max_absolute_error": error,
                                          "passes": passed}
            emit(
                f"| {discretization} | {row['dofs']} | {values[0]:.6f} | {values[1]:.6f} | "
                f"{values[2]:.6f} | {error:.4f} | {'yes' if passed else 'NO'} |"
            )
        report["published_phase_fractions"] = acceptance

        by_mesh: dict[str, list[str]] = {}
        for discretization in phases:
            by_mesh.setdefault(_mesh_of(discretization), []).append(discretization)
        p_control_identical = all(
            len({phases[key] for key in group}) == 1 for group in by_mesh.values()
        )
        mesh_stacked = np.asarray([phases[group[0]] for group in by_mesh.values()], dtype=float)
        absolute_drift = float(np.max(np.ptp(mesh_stacked, axis=0)))
        relative_drift = float(
            100.0 * np.max(np.ptp(mesh_stacked, axis=0) / np.mean(mesh_stacked, axis=0))
        )
        emit()
        emit(f"- p-refinement null control: phase fractions identical at fixed mesh across "
             f"P1/P2 -- **{p_control_identical}**. They must be: the groups are element sets, "
             "so raising the element order cannot resample the material field.")
        emit(f"- h-refinement: largest spread across meshes **{absolute_drift:.4f}** in absolute "
             f"volume fraction, **{relative_drift:.2f} %** relative.")
        failed = [key for key, entry in acceptance.items() if not entry["passes"]]
        if failed:
            emit(f"- the project's own material-fidelity rule (max abs(phase - "
                 f"{list(PHASE_TARGET)}) <= {PHASE_TOLERANCE}) is applied upstream to the "
                 f"primary mesh only. Applied to every row it FAILS for: "
                 f"{', '.join(sorted(failed))}. Those discretisations carry a different "
                 "material distribution from the one the paper accepts.")
        report["phase_drift"] = {
            "p_refinement_identical": bool(p_control_identical),
            "absolute_drift": absolute_drift,
            "relative_drift_percent": relative_drift,
            "acceptance_failures": sorted(failed),
            "meshes_compared": int(mesh_stacked.shape[0]),
        }

        emit()
        emit("### Richardson limit against the same-mesh P2 point")
        emit()
        gaps = {}
        for condition in ("Control", "Hyperosmotic"):
            series = [
                (
                    MESH_SIZE_UM[_mesh_of(key)],
                    float(summary_rows[(key, condition)][ENDPOINT_COLUMN]),
                )
                for key in phases
                if _order_of(key) == "P1" and _mesh_of(key) in MESH_SIZE_UM
            ]
            estimate = _richardson(series)
            p2_key = f"{args.primary_mesh}_P2"
            p2_value = (
                float(summary_rows[(p2_key, condition)][ENDPOINT_COLUMN])
                if (p2_key, condition) in summary_rows
                else None
            )
            if estimate and estimate.get("valid") and p2_value is not None:
                gap = 100.0 * (p2_value - estimate["extrapolated"]) / p2_value
                gaps[condition] = gap
                emit(
                    f"- **{condition}** -- P1 observed order {estimate['observed_order']:.2f}, "
                    f"Richardson limit {estimate['extrapolated']:.6f} mN/m; {p2_key} gives "
                    f"{p2_value:.6f} mN/m; gap **{gap:+.2f} %** of the P2 value."
                )
            elif estimate and not estimate.get("valid"):
                emit(f"- **{condition}** -- no Richardson limit ({estimate['reason']}).")
            else:
                emit(f"- **{condition}** -- too few P1 meshes for an h-study.")
        report["richardson_gap_percent"] = gaps

    emit()
    emit("## 3. Material resampling measured on the real meshes")
    emit()
    resampling: dict[str, dict] = {}
    if args.dry_run:
        emit("--dry-run: layer B skipped. It would call the project's own "
             "`load_reynolds_cellwise_heterogeneity` on "
             f"{len(meshes)} mesh(es) with group_count={args.group_count}, "
             f"shear {args.minimum_shear_pa:g}-{args.maximum_shear_pa:g} Pa, "
             f"sigmoid ({args.sigmoid_midpoint:g}, {args.sigmoid_width:g}), "
             f"K/G={args.bulk_to_shear_ratio:g}, and report per-group volume fractions, "
             "Voigt/Reuss effective shear moduli and their drift across meshes. "
             "It needs meshio and src/nuclear_envelope_fem.")
    elif not meshes or volume_path is None:
        emit("Meshes or the chromatin volume were not found, so layer B cannot run. "
             "Point --project-root at the real working tree, or pass --mesh-path and "
             "--template-volume.")
    else:
        meshio = None
        try:
            meshio = _lazy_import("meshio", "reading the gmsh meshes (layer B)")
            load_heterogeneity, tags = _import_project(args.project_root)
        except SystemExit as exc:
            # Layer A has already produced numbers; report the gap instead of dying on it.
            emit(str(exc))
            meshio = None
            report["layer_b_unavailable"] = str(exc)
        targets = meshes if meshio is not None else {}
        for name, mesh_path in targets.items():
            mesh = meshio.read(mesh_path)
            points = np.asarray(mesh.points, dtype=float)[:, :3]
            cells = None
            cell_tags = None
            for block, data in zip(
                mesh.cells, mesh.cell_data.get("gmsh:physical", []), strict=False
            ):
                if block.type == "tetra":
                    cells = np.asarray(block.data, dtype=np.int64)
                    cell_tags = np.asarray(data, dtype=np.int32)
            if cells is None:
                emit(f"- {name}: no tetrahedral block in `{mesh_path}`, skipped.")
                continue
            heterogeneity = load_heterogeneity(
                str(volume_path),
                points,
                cells,
                cell_tags,
                inner_axes,
                group_count=args.group_count,
                minimum_shear_pa=args.minimum_shear_pa,
                maximum_shear_pa=args.maximum_shear_pa,
                sigmoid_midpoint=args.sigmoid_midpoint,
                sigmoid_width=args.sigmoid_width,
                bulk_to_shear_ratio=args.bulk_to_shear_ratio,
            )
            core = cell_tags == tags.CORE_VOLUME
            volumes = _tetra_volumes(points, cells[core])
            groups = np.asarray(heterogeneity.group_ids, dtype=np.int64)[core]
            shear = _group_shear_moduli(heterogeneity)
            fractions = np.asarray(
                [volumes[groups == index].sum() / volumes.sum() for index in range(len(shear))]
            )
            present = fractions > 0.0
            voigt = float(np.sum(fractions[present] * shear[present]))
            reuss = float(1.0 / np.sum(fractions[present] / shear[present]))
            phase_fractions = {
                label: float(fractions[start:stop].sum())
                for label, (start, stop) in PHASE_GROUP_SLICES.items()
            }
            resampling[name] = {
                "mesh": str(mesh_path),
                "core_tetrahedra": int(core.sum()),
                "mean_tetrahedron_volume_um3": float(volumes.mean()),
                "voxels_per_tetrahedron": float(volumes.mean() / voxel_volume),
                "group_volume_fractions": fractions.tolist(),
                "group_shear_moduli_pa": shear.tolist(),
                "phase_volume_fractions": phase_fractions,
                "voigt_shear_pa": voigt,
                "reuss_shear_pa": reuss,
            }
            if name in mesh_metrics:
                mesh_metrics[name]["mean_tetrahedron_volume_um3"] = float(volumes.mean())
                mesh_metrics[name]["voxels_per_tetrahedron"] = float(volumes.mean() / voxel_volume)
        if resampling:
            emit("| mesh | core tets | mean tet volume (um^3) | voxels/tet | Voigt G (Pa) "
                 "| Reuss G (Pa) |")
            emit("|---|---:|---:|---:|---:|---:|")
            for name, entry in resampling.items():
                emit(
                    f"| {name} | {entry['core_tetrahedra']} | "
                    f"{entry['mean_tetrahedron_volume_um3']:.6f} | "
                    f"{entry['voxels_per_tetrahedron']:.2f} | {entry['voigt_shear_pa']:.1f} | "
                    f"{entry['reuss_shear_pa']:.1f} |"
                )
            emit()
            emit("| group | " + " | ".join(resampling) + " |")
            emit("|---|" + "---:|" * len(resampling))
            group_count = len(next(iter(resampling.values()))["group_volume_fractions"])
            for index in range(group_count):
                cells_text = " | ".join(
                    f"{entry['group_volume_fractions'][index]:.6f}" for entry in resampling.values()
                )
                emit(f"| {index + 1} | {cells_text} |")
            voigts = np.asarray([entry["voigt_shear_pa"] for entry in resampling.values()])
            reusses = np.asarray([entry["reuss_shear_pa"] for entry in resampling.values()])
            voigt_drift = float(100.0 * np.ptp(voigts) / np.mean(voigts))
            reuss_drift = float(100.0 * np.ptp(reusses) / np.mean(reusses))
            emit()
            emit(f"- Voigt effective core shear modulus drifts **{voigt_drift:.2f} %** across "
                 f"the meshes ({voigts.min():.1f} to {voigts.max():.1f} Pa).")
            emit(f"- Reuss effective core shear modulus drifts **{reuss_drift:.2f} %** "
                 f"({reusses.min():.1f} to {reusses.max():.1f} Pa).")
            emit("- Identical meshes would give identical numbers. These are the same "
                 "chromatin volume and the same loader; only the tetrahedra changed.")
            report["resampling_drift_percent"] = {"voigt": voigt_drift, "reuss": reuss_drift}
    report["resampling"] = resampling

    emit()
    emit("## 4. Quadrature sweep at fixed mesh")
    emit()
    sweep: dict[str, dict] = {}
    if not args.quadrature_sweep:
        emit("Not requested. Add --quadrature-sweep --out-root DIR to run it.")
    else:
        discretization = f"{args.primary_mesh}_{args.element_order.upper()}"
        emit(f"Fixed mesh `{args.primary_mesh}`, fixed element order "
             f"{args.element_order.upper()}, quadrature_order in "
             f"{list(args.quadrature_orders)} via _common.patch_options.")
        emit("Each order is a full Control/Hyperosmotic pair into "
             f"`{args.out_root or '<--out-root>'}/quadrature_<order>`.")
        if args.dry_run:
            emit()
            emit("--dry-run: nothing solved.")
        else:
            runner = _runner(args.project_root)
            args.out_root.mkdir(parents=True, exist_ok=True)
            for order in args.quadrature_orders:
                out_dir = args.out_root / f"quadrature_{order}"
                print(f"\n=== quadrature_order={order} -> {out_dir} ===", flush=True)
                original = patch_options(runner, quadrature_order=int(order))
                try:
                    runner.run_analysis(
                        out_dir,
                        primary_mesh=args.primary_mesh,
                        run_coarse_p2=False,
                        run_intermediate_p2=(args.element_order == "p2"),
                        resume=args.resume,
                        sparse_backend=args.sparse_backend,
                        rigid_constraint_linear_solver=args.rigid_constraint_linear_solver,
                    )
                finally:
                    runner._options = original
                try:
                    rows = read_endpoints(out_dir)
                    control = float(rows[(discretization, "Control")][ENDPOINT_COLUMN])
                    hyper = float(rows[(discretization, "Hyperosmotic")][ENDPOINT_COLUMN])
                except (FileNotFoundError, KeyError) as exc:
                    emit(f"- quadrature_order {order}: could not read the endpoint ({exc}).")
                    continue
                sweep[str(order)] = {
                    "out_dir": str(out_dir),
                    "control_mn_per_m": control,
                    "hyperosmotic_mn_per_m": hyper,
                    "difference_mn_per_m": hyper - control,
                }
            if sweep:
                baseline = sweep[str(args.quadrature_orders[0])]
                emit()
                emit("| quadrature_order | Control (mN/m) | Hyperosmotic (mN/m) | difference "
                     "(mN/m) | Control movement |")
                emit("|---:|---:|---:|---:|---:|")
                for order, entry in sweep.items():
                    movement = 100.0 * (
                        entry["control_mn_per_m"] - baseline["control_mn_per_m"]
                    ) / baseline["control_mn_per_m"]
                    entry["control_movement_percent"] = movement
                    emit(
                        f"| {order} | {entry['control_mn_per_m']:.6f} | "
                        f"{entry['hyperosmotic_mn_per_m']:.6f} | "
                        f"{entry['difference_mn_per_m']:+.6f} | {movement:+.3f} % |"
                    )
                emit()
                emit("Interpretation: the material field is piecewise constant per tetrahedron, "
                     "so quadrature refinement cannot resample it. Movement here is integration "
                     "error at frozen material and frozen geometry -- partly in the outer-surface "
                     "post-processing, which uses max(4, quadrature_order + 1). Subtract it from "
                     "the h-family drift before attributing that drift to discretisation.")
    report["quadrature_sweep"] = sweep

    emit()
    emit("## Verdict")
    emit()
    drift_percent = None
    if "resampling_drift_percent" in report:
        drift_percent = max(report["resampling_drift_percent"].values())
        drift_source = "Voigt/Reuss effective core shear modulus across the meshes"
    elif "phase_drift" in report:
        drift_percent = report["phase_drift"]["relative_drift_percent"]
        drift_source = "published phase volume fractions across the meshes"
    gap_percent = None
    if report.get("richardson_gap_percent"):
        gap_percent = abs(report["richardson_gap_percent"].get("Control", float("nan")))

    if drift_percent is None:
        emit("Not enough inputs were available to reach a verdict. Layer A needs only the "
             f"{SUMMARY_NAME}; run it first.")
    elif gap_percent is not None and np.isfinite(gap_percent):
        verdict = (
            f"The chromatin material field changes by {drift_percent:.2f} % "
            f"({drift_source}) between the meshes used in the h-refinement study, against a "
            f"{gap_percent:.2f} % discrepancy between the P1 Richardson limit and the "
            f"same-mesh P2 endpoint. "
            + (
                "The material change is of the same order as, or larger than, the "
                "discrepancy attributed to discretisation error, so the h-family is not a "
                "sequence of discretisations of a single continuous problem and its "
                "Richardson extrapolation and GCI cannot be read as discretisation error."
                if drift_percent >= gap_percent
                else "The material change is smaller than the discrepancy, so resampling "
                "explains only part of it and the remainder is genuine discretisation error; "
                "report both rather than the GCI alone."
            )
        )
        emit(verdict)
        report["verdict"] = verdict
        emit()
        emit("Paste-ready limitations sentence:")
        emit()
        emit("> Chromatin stiffness is assigned per tetrahedron by trilinear sampling of a "
             f"{grid[0]}^3 voxel volume at element centroids, so mesh refinement resamples the "
             "material field rather than holding it fixed; across the meshes used in the "
             f"convergence study the realised material distribution changes by "
             f"{drift_percent:.2f} %, against the {gap_percent:.2f} % spread between "
             "discretisations that the study attributes to numerical error, and the reported "
             "mesh study should therefore be read as a sensitivity bound rather than as "
             "convergence to a fixed continuous problem.")
    else:
        emit(f"The material field changes by {drift_percent:.2f} % ({drift_source}) across the "
             "meshes used in the convergence study. No Richardson/P2 comparison was available "
             "here, so quote the drift alongside the endpoint spread.")

    if args.json_out:
        args.json_out.write_text(json.dumps(report, indent=2) + "\n", encoding="utf-8")
        print(f"\nWrote {args.json_out}", file=sys.stderr)
    if args.markdown_out:
        args.markdown_out.write_text("\n".join(lines) + "\n", encoding="utf-8")
        print(f"Wrote {args.markdown_out}", file=sys.stderr)
    return 0


if __name__ == "__main__":
    raise SystemExit(main())