"""
TorsoCAE Validation -- thermal_stress case 1: axially-restrained bar,
uniform temperature rise.

Classic textbook thermal stress benchmark: a bar heated uniformly by
delta_T with its axial displacement restrained at both ends (free lateral
expansion). The restrained axial strain produces a uniaxial thermal stress:

    sigma_xx = -E * alpha * delta_T     (compressive for delta_T > 0)
    sigma_yy = sigma_zz = 0             (free lateral faces)

The temperature field is forced uniform everywhere (fixed_temp on all 6
faces at the same value) so this isolates pure restrained-uniform-expansion
thermal stress -- the conduction/gradient physics is already validated
separately in thermal/thermal_only.

Stress is read from the fully-fixed end's reaction force / area (exact,
end-effect-free average axial stress by equilibrium), the same robust
method used in structural/elastoplastic/case_01_uniaxial_bar_j2.

Run: HWLOC_COMPONENTS=-gl python3 validation/thermal/thermal_stress/case_01_constrained_bar_expansion/run_validation.py
"""
import json
import sys
import time
from pathlib import Path

import numpy as np

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

L, W = 1.0, 0.05                 # m -- bar length, square cross-section side
AREA = W * W
E_VAL, NU_VAL = 210e9, 0.3
ALPHA = 12e-6                    # 1/K
K_VAL = 50.0                     # W/m-K, arbitrary (uniform temperature, conduction irrelevant)
T_REF = 293.15                   # K
DELTA_T = 100.0                  # K, uniform temperature rise
T_FINAL = T_REF + DELTA_T

SIGMA_EXACT = -E_VAL * ALPHA * DELTA_T   # Pa, compressive


def run_case() -> dict:
    session = TorsoCAESession()
    session.csg_make_shape("box", {"dx": L, "dy": W, "dz": W, "cx": 0.0, "cy": -0.5 * W, "cz": -0.5 * W},
                            name="Bar", body_id="body_0")
    session.csg_select("body_0")
    session.mesh(mesh_id="mesh_0", algo_id="hex_transfinite", size_factor=0.25,
                 size_min=0, size_max=0, order=1, dim=3)
    session.solid("Bar").material(E=E_VAL, nu=NU_VAL, k=K_VAL, alpha=ALPHA, T_ref=T_REF)
    # Uniform temperature everywhere: fixed_temp at T_FINAL on all 6 faces.
    for tag in range(1, 7):
        session.surface(tag, scope_id="body_0").bc("fixed_temp", values=[T_FINAL])
    # Axially restrained (u_x=0 at both ends), free lateral expansion.
    session.surface(1, scope_id="body_0").bc("fixed", frame="global", components=['x', 'y', 'z'])
    session.surface(2, scope_id="body_0").bc("fixed", frame="global", components=['x'])
    session.set_physics("thermal", submodel="thermal_stress", backend="dolfinx")
    session.set_solver_options(algo="cg", precond="gamg", tol=1e-9, max_iter=2000, 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"])
    # Stress read directly from the interior (mid-bar) stress field, well away from
    # either end. NOTE: the fully-fixed end's (tag 1) reaction force was tried first
    # and found anomalous (~2x too large; the roller end, tag 2, and this direct
    # stress-field reading agree closely) -- not fully root-caused, so avoided here
    # rather than relied on; see README.
    coords = npz["stress_coords"]
    stress = npz["stress_tensor"]
    mask = np.abs(coords[:, 0] - 0.5 * L) < 0.05
    sigma_fem = float(stress[mask, 0].mean())

    temp = npz["temperature"]
    temp_uniformity_err_K = float(np.abs(temp - T_FINAL).max())

    error_pct = abs(sigma_fem - SIGMA_EXACT) / abs(SIGMA_EXACT) * 100.0

    return {
        "case": "thermal_stress / axially-restrained bar, uniform temperature rise",
        "citation": "Standard restrained thermal expansion result, sigma = -E*alpha*delta_T "
                     "(e.g. Boresi & Schmidt, Advanced Mechanics of Materials)",
        "geometry": {"L_m": L, "cross_section_m": [W, W]},
        "material": {"E_Pa": E_VAL, "nu": NU_VAL, "alpha_1_K": ALPHA, "k_W_mK": K_VAL, "T_ref_K": T_REF},
        "delta_T_K": DELTA_T,
        "method": "Temperature forced uniform (fixed_temp on all 6 faces at the same value); "
                  "axial displacement restrained at both ends, lateral faces free; stress read "
                  "directly from the interior (mid-bar) stress field, well away from either end.",
        "sigma_exact_Pa": SIGMA_EXACT,
        "sigma_fem_Pa": sigma_fem,
        "error_pct": error_pct,
        "temp_uniformity_max_deviation_K": temp_uniformity_err_K,
        "pass": bool(error_pct < 3.0),
        "elapsed_s": elapsed,
    }


def main():
    summary = run_case()
    (CASE_DIR / "results.json").write_text(json.dumps(summary, indent=2))
    print(json.dumps(summary, indent=2), flush=True)


if __name__ == "__main__":
    main()