"""
TorsoCAE Validation -- elastoplastic case 1: uniaxial tension bar loaded
past yield, J2 isotropic hardening (submodel="elastoplastic-j2").

Closed form, derived from the code's ACTUAL hardening law
(server/fem_backends/plasticity_j2.py, class HardeningLaw):

    R(alpha) = H*alpha + Q*(1 - exp(-b*alpha))    # alpha = accum. eq. plastic strain

With Q=0, b=0 (isotropic-hardening-only case used here), R(alpha) = H*alpha,
i.e. H is the PLASTIC modulus dR/dalpha = dsigma_yield/dalpha, NOT the
total-strain tangent modulus. For monotonic uniaxial tension, alpha equals
the axial plastic strain eps_p exactly (J2 associative flow), so:

    eps = eps_e + eps_p = sigma/E + alpha
    sigma = sigma_y0 + H*alpha = sigma_y0 + H*(eps - sigma/E)
    => sigma*(1 + H/E) = sigma_y0 + H*eps
    => sigma = (sigma_y0 + H*eps) * E / (E + H)          [eps >= eps_y]
    sigma = E*eps                                         [eps <  eps_y]
    eps_y = sigma_y0 / E

This is algebraically the same relation as the "sigma_y + H*(eps-sigma_y/E)"
bilinear form (both reduce to a post-yield tangent slope E_t = E*H/(E+H)
and agree exactly at eps=eps_y), so both forms are equivalent -- we use the
E*(sigma_y0+H*eps)/(E+H) form directly since it's the most direct algebraic
rearrangement of the code's R(alpha)=H*alpha law.

BCs mirror the existing working template (journals/elastoplastic_bar_j2_chaboche1.py):
full clamp (all 3 components) at x=0, full displacement BC at x=L. This
constrains Poisson contraction at both end faces (grip-like end effect),
so per that template's own documented caveat we do NOT read stress from
nodal fields at the ends. Instead we read the FIXED-END REACTION FORCE
(saved in the npz as reaction_tags/reaction_forces -- exact integral of
sigma.n over that BC surface), which by simple axial equilibrium of a bar
under no body force equals the internal axial force N at EVERY cross
section, including the mid-span far from end effects. sigma_FEM = N/A is
therefore the correct, end-effect-free average axial stress at mid-span,
without needing to hunt down mid-span nodes/quadrature points explicitly.
Nominal (average) strain = target_disp / L.

Run: HWLOC_COMPONENTS=-gl python3 validation/structural/elastoplastic/case_01_uniaxial_bar_j2/run_validation.py
"""
import json
import sys
import time
from pathlib import Path

import numpy as np

sys.path.insert(0, str(Path(__file__).resolve().parents[4] / "server"))
from pytorsocae import TorsoCAESession

CASE_DIR = Path(__file__).resolve().parent

E_VAL, NU_VAL = 210e9, 0.3
SIGMA_Y0, HARDENING_H = 250e6, 2e9
L, W = 500.0, 50.0   # slender bar, matches journals/elastoplastic_bar_j2_chaboche1.py
AREA = W * W
EPS_Y = SIGMA_Y0 / E_VAL

# strains spanning elastic -> well past yield
TARGET_STRAINS = [0.0008, 0.0013, 0.0018, 0.0028, 0.0042]


def sigma_closed_form(eps: float) -> float:
    if eps < EPS_Y:
        return E_VAL * eps
    return (SIGMA_Y0 + HARDENING_H * eps) * E_VAL / (E_VAL + HARDENING_H)


def run_case(eps_target: float) -> dict:
    target_disp = eps_target * L

    session = TorsoCAESession()
    session.csg_make_shape("box", {"dx": L, "dy": W, "dz": W, "cx": 0, "cy": 0, "cz": 0},
                            name="bar", body_id="body_0")
    session.csg_select("body_0")
    session.mesh(mesh_id="mesh_0", algo_id="hex_transfinite", size_factor=0.6,
                 size_min=0, size_max=0, order=1, dim=3)
    session.solid("bar").material(E=E_VAL, nu=NU_VAL, yield_stress=SIGMA_Y0, hardening_H=HARDENING_H)
    session.surface(1, scope_id="body_0").bc("fixed")
    session.surface(2, scope_id="body_0").bc("displacement", values=[target_disp, 0, 0])
    session.set_physics("structural", submodel="elastoplastic-j2", backend="dolfinx")
    session.set_solver_options(num_steps=20, inner_tol=1e-8, max_inner_iter=30, n_cores=2, device="cpu")
    t0 = time.time()
    result = session.compute(mesh_ids=["mesh_0"])
    elapsed = time.time() - t0

    npz = np.load(result["npz_path"])
    reaction_tags = npz["reaction_tags"]
    reaction_forces = npz["reaction_forces"]
    idx = int(np.where(reaction_tags == 1)[0][0])
    N_reaction = float(reaction_forces[idx, 0])
    sigma_fem = abs(N_reaction) / AREA

    sigma_ref = sigma_closed_form(eps_target)
    return {
        "eps_target": eps_target,
        "sigma_fem_Pa": sigma_fem,
        "sigma_closed_form_Pa": sigma_ref,
        "reaction_force_N": N_reaction,
        "max_von_mises_Pa": result.get("max_von_mises"),
        "error_pct": abs(sigma_fem - sigma_ref) / sigma_ref * 100.0,
        "elapsed_s": elapsed,
    }


def main():
    rows = [run_case(e) for e in TARGET_STRAINS]
    max_err = max(r["error_pct"] for r in rows)
    summary = {
        "case": "elastoplastic / uniaxial bar, J2 isotropic hardening",
        "material": {"E_Pa": E_VAL, "nu": NU_VAL, "yield_stress_Pa": SIGMA_Y0, "hardening_H_Pa": HARDENING_H},
        "geometry": {"L_m": L, "cross_section_m": [W, W]},
        "eps_yield": EPS_Y,
        "method": "closed form derived from code's R(alpha)=H*alpha isotropic hardening law "
                  "(plasticity_j2.py HardeningLaw); stress read from fixed-end reaction force / "
                  "area (equilibrium-exact, end-effect free average axial stress)",
        "points": rows,
        "max_error_pct": max_err,
        "pass": max_err < 5.0,
    }
    (CASE_DIR / "results.json").write_text(json.dumps(summary, indent=2))
    print(json.dumps(summary, indent=2), flush=True)
    plot_curve(rows)


def plot_curve(rows: list[dict]) -> None:
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    eps_curve = np.linspace(0, max(r["eps_target"] for r in rows) * 1.05, 300)
    sigma_curve = [sigma_closed_form(e) / 1e6 for e in eps_curve]

    fig, ax = plt.subplots(figsize=(6.4, 4.4), dpi=150)
    ax.plot(eps_curve * 100, sigma_curve, "-", color="#dc2626", linewidth=2,
             label="Closed form (code's R(alpha)=H*alpha)")
    ax.plot([r["eps_target"] * 100 for r in rows], [r["sigma_fem_Pa"] / 1e6 for r in rows],
             "x", color="#2563eb", markersize=9, markeredgewidth=2,
             label="TorsoCAE FEM (reaction force / area)")
    ax.axvline(EPS_Y * 100, color="gray", linestyle=":", linewidth=1, label=f"yield strain {EPS_Y*100:.3f}%")
    ax.set_xlabel("Nominal axial strain (%)")
    ax.set_ylabel("Axial stress (MPa)")
    ax.set_title("J2 elastoplastic uniaxial bar: FEM vs. closed form")
    ax.legend(loc="upper left")
    ax.grid(alpha=0.3)
    fig.tight_layout()
    fig.savefig(CASE_DIR / "plots" / "stress_strain.png")


if __name__ == "__main__":
    main()