"""
TorsoCAE Validation -- thermal_only case 1: 1D slab conduction.

Classic textbook steady-state heat conduction benchmark: a slab with fixed
temperatures on its two opposing x-faces and insulated (adiabatic, the
submodel's default) side faces. The exact solution is the linear profile
T(x) = T_hot + (T_cold - T_hot) * x / L, with uniform heat flux
q = k * (T_hot - T_cold) / L throughout.

Because the exact solution is LINEAR in x and TorsoCAE's thermal_only
submodel uses linear (P1) Lagrange elements by default, a correctly
implemented FEM solve should reproduce this profile to machine precision
(any P1 field can represent an exact linear function exactly) -- this is
therefore primarily a solver-correctness check, in the same spirit as the
homogeneous-state hyperelastic/elastoplastic cases elsewhere in this suite.

Run: HWLOC_COMPONENTS=-gl python3 validation/thermal/thermal_only/case_01_1d_slab_conduction/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

LX, LY, LZ = 1.0, 0.1, 0.1     # m, slab along x
K = 50.0                        # W/m-K
T_HOT, T_COLD = 400.0, 300.0    # K

Q_EXACT = K * (T_HOT - T_COLD) / LX   # W/m^2, uniform flux


def t_exact(x: np.ndarray) -> np.ndarray:
    return T_HOT + (T_COLD - T_HOT) * x / LX


def run_case() -> dict:
    session = TorsoCAESession()
    session.csg_make_shape("box", {"dx": LX, "dy": LY, "dz": LZ, "cx": 0.0, "cy": -0.5 * LY, "cz": -0.5 * LZ},
                            name="Slab", body_id="body_0")
    session.csg_select("body_0")
    session.mesh(mesh_id="mesh_0", algo_id="hex_transfinite", size_factor=0.5,
                 size_min=0, size_max=0, order=1, dim=3)
    session.solid("Slab").material(k=K)
    session.surface(1, scope_id="body_0").bc("fixed_temp", values=[T_HOT])
    session.surface(2, scope_id="body_0").bc("fixed_temp", values=[T_COLD])
    session.set_physics("thermal", submodel="thermal_only", backend="dolfinx")
    session.set_solver_options(algo="cg", precond="gamg", tol=1e-12, max_iter=1000, 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"])
    coords = npz["coordinates"]
    temp = npz["temperature"]
    t_ref = t_exact(coords[:, 0])
    temp_err_pct = np.abs(temp - t_ref) / (T_HOT - T_COLD) * 100.0
    max_temp_err_pct = float(temp_err_pct.max())

    flux = npz["heat_flux"]
    flux_err_pct = float(np.abs(flux - Q_EXACT).max() / Q_EXACT * 100.0)

    return {
        "case": "thermal_only / 1D slab conduction",
        "citation": "Standard 1D steady-state Fourier conduction (textbook exact solution)",
        "geometry": {"Lx_m": LX, "Ly_m": LY, "Lz_m": LZ},
        "material": {"k_W_mK": K},
        "bc": {"T_hot_K": T_HOT, "T_cold_K": T_COLD},
        "method": "FEM temperature field at every mesh node compared to the exact linear "
                  "profile T(x) = T_hot + (T_cold-T_hot)*x/L; FEM heat flux magnitude "
                  "(saved per-cell in the npz) compared to the exact uniform flux "
                  "q = k*(T_hot-T_cold)/L.",
        "q_exact_W_m2": Q_EXACT,
        "max_temp_error_pct": max_temp_err_pct,
        "max_flux_error_pct": flux_err_pct,
        "n_nodes": int(coords.shape[0]),
        "pass": bool(max_temp_err_pct < 1.0 and flux_err_pct < 1.0),
        "elapsed_s": elapsed,
        "profile": {
            "x_m": coords[:, 0].tolist(),
            "t_fem_K": temp.tolist(),
        },
    }


def plot_profile(summary: dict) -> None:
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    prof = summary["profile"]
    x = np.array(prof["x_m"])
    t_fem = np.array(prof["t_fem_K"])
    x_curve = np.linspace(0, LX, 200)

    fig, ax = plt.subplots(figsize=(6.4, 4.4), dpi=150)
    ax.plot(x_curve, t_exact(x_curve), "-", color="#dc2626", linewidth=2, label="Exact linear profile")
    ax.plot(x, t_fem, "o", color="#2563eb", markersize=3, alpha=0.5, label="TorsoCAE FEM (all nodes)")
    ax.set_xlabel("x (m)")
    ax.set_ylabel("Temperature (K)")
    ax.set_title(f"1D slab conduction: FEM vs. exact (max err {summary['max_temp_error_pct']:.2e}%)")
    ax.legend(loc="best")
    ax.grid(alpha=0.3)
    fig.tight_layout()
    fig.savefig(CASE_DIR / "plots" / "temperature_profile.png")


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


if __name__ == "__main__":
    main()