"""
TorsoCAE Validation -- thermal_stress case 2: Timoshenko (1925) bimetallic
strip curvature under uniform heating.

This is the case originally assigned in validationlist.txt. An earlier pass
this session substituted a single-material thermal-bending-gradient case
(case_02_thermal_bending_gradient/) after finding dolfinx_backend.py's
thermoelastic driver (_solve_thermal_with_stress) took E/nu/alpha/k as plain
scalars with no per-region variation, and that csg_boolean("union", ...) of
two differently-materialed bodies collapsed to a single merged volume tag.
That substitution was reconsidered after being challenged (correctly): the
"two volumes" framing was the wrong axis to force -- a bimetallic strip's
material discontinuity is a function of position (z), not of mesh topology,
and the codebase already supports spatial EXPRESSION strings for BC values
(ConditionEvaluator, e.g. inlet_velocity="..."). Testing (not just reading
the code) confirmed material params did NOT yet support this: an
expression-valued alpha crashed at `Constant(mesh, float(alpha))`.

Fix (dolfinx_backend.py): added `_material_coeff(mesh, value, bbox)`, which
returns a Constant for a plain number or an interpolated P1 Function for an
expression string, and `_mesh_bbox(...)` (same pattern already used for
inlet_velocity/traction/pressure BC expressions). Wired into
_solve_thermal_with_stress for k, E, nu, alpha -- Lame constants (lam, mu_c)
are now derived via UFL arithmetic on these coefficients, which works
identically whether they are Constants or spatially-varying Functions. This
is a small, scoped port of a pattern the codebase already used for BC
values, not a new subsystem -- a single mesh volume can now carry a
spatially-discontinuous material, which is exactly what a bonded bimetallic
strip needs (no second mesh volume, no CSG union, no interface at all).

Theory (equal-thickness, equal-E composite beam under uniform deltaT,
derived here from first principles via beam equilibrium -- N=0, M=0 -- and
cross-checked against the well-known n=m=1 Timoshenko special case):

    kappa = 1.5 * (alpha_bottom - alpha_top) * deltaT / h

with the tip curling toward the LOWER-expansion layer (same sense as any
bimetallic strip). E is kept equal between the two regions specifically so
the exact, unambiguous n=1 special case of Timoshenko's formula applies
(the general unequal-modulus formula was not used here to avoid the risk of
misremembering its more complex coefficients from memory without an
independent derivation to check them against).

Small-deflection cantilever: delta_tip_z = +kappa * L^2 / 2.

Run: HWLOC_COMPONENTS=-gl python3 validation/thermal/thermal_stress/case_02_timoshenko_bimetal_strip/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 = 0.5, 0.02, 0.02      # m -- cantilever length, width, total (both-layer) thickness
E_VAL, NU_VAL = 150e9, 0.3         # equal for both layers -- isolates the n=1 Timoshenko special case
ALPHA_BOTTOM = 19e-6                # 1/K, z < 0  (brass-like -- higher expansion)
ALPHA_TOP    = 12e-6                # 1/K, z > 0  (steel-like -- lower expansion)
K_VAL = 100.0                       # W/m-K, uniform (irrelevant -- temperature forced uniform via BC)
T_REF = 293.15                      # K
DELTA_T = 100.0                     # K, uniform temperature rise
T_FINAL = T_REF + DELTA_T

ALPHA_EXPR = f"{ALPHA_TOP} + {ALPHA_BOTTOM - ALPHA_TOP}*(z<0)"   # vectorized step, not if/else

KAPPA = 1.5 * (ALPHA_BOTTOM - ALPHA_TOP) * DELTA_T / LZ   # 1/m
TIP_DZ_EXACT = KAPPA * LX ** 2 / 2.0                       # m, curls toward the lower-expansion (top) layer

# Plain (non-boolean) csg box: 1=x=0, 2=x=dx, 3=y=cy, 4=y=cy+dy, 5=z=cz, 6=z=cz+dz
TAG_ROOT = 1


def run_case(size_factor: float) -> 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="Strip", body_id="body_0")
    session.csg_select("body_0")
    # order=2: same bending-locking finding as case_02_thermal_bending_gradient applies
    # identically here (thin, pure-bending-dominated beam) -- start with quadratic elements.
    session.mesh(mesh_id="mesh_0", algo_id="hex_transfinite", size_factor=size_factor,
                 size_min=0, size_max=0, order=2, dim=3)
    session.solid("Strip").material(E=E_VAL, nu=NU_VAL, k=K_VAL, alpha=ALPHA_EXPR, T_ref=T_REF)
    for tag in range(1, 7):
        session.surface(tag, scope_id="body_0").bc("fixed_temp", values=[T_FINAL])
    session.surface(TAG_ROOT, scope_id="body_0").bc("fixed", frame="global", components=['x', 'y', 'z'])
    session.set_physics("thermal", submodel="thermal_stress", backend="dolfinx")
    session.set_solver_options(algo="direct", precond="none", 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"])
    coords = npz["coordinates"]
    disp = npz["displacement"]
    temp = npz["temperature"]

    tip_mask = np.abs(coords[:, 0] - LX) < 1e-9
    tip_dz_fem = float(disp[tip_mask, 2].mean())
    tip_error_pct = abs(tip_dz_fem - TIP_DZ_EXACT) / abs(TIP_DZ_EXACT) * 100.0
    temp_uniformity_max_deviation_K = float(np.abs(temp - T_FINAL).max())

    return {
        "size_factor": size_factor,
        "n_dofs": int(coords.shape[0]),
        "tip_dz_exact_m": TIP_DZ_EXACT,
        "tip_dz_fem_m": tip_dz_fem,
        "tip_error_pct": tip_error_pct,
        "temp_uniformity_max_deviation_K": temp_uniformity_max_deviation_K,
        "max_von_mises_pa": result.get("max_von_mises"),
        "elapsed_s": elapsed,
    }


def main():
    runs = [run_case(sf) for sf in (1.0, 0.5)]
    finest = runs[-1]

    summary = {
        "case": "thermal_stress / Timoshenko (1925) bimetallic strip curvature under uniform heating",
        "citation": "Timoshenko, S. (1925) 'Analysis of Bi-Metal Thermostats', J. Opt. Soc. Am. 11(3); "
                     "equal-thickness equal-modulus (n=m=1) special case kappa=1.5*(alpha1-alpha2)*deltaT/h, "
                     "derived here independently via composite-beam equilibrium (N=0, M=0) and cross-checked "
                     "against that known special case.",
        "backend_fix": "Added expression-valued material coefficient support to dolfinx_backend.py's "
                        "_solve_thermal_with_stress (_material_coeff/_mesh_bbox helpers), mirroring the "
                        "existing BC-value expression pattern (ConditionEvaluator) -- lets alpha (and E, nu, "
                        "k) vary spatially within a single mesh volume, which is what a bonded bimetallic "
                        "strip needs. See module docstring for the full story (this replaced an earlier "
                        "single-material substitute case after being directly challenged and re-tested).",
        "geometry": {"L_m": LX, "width_m": LY, "total_thickness_m": LZ},
        "material": {"E_Pa": E_VAL, "nu": NU_VAL, "alpha_bottom_1_K": ALPHA_BOTTOM,
                     "alpha_top_1_K": ALPHA_TOP, "k_W_mK": K_VAL, "T_ref_K": T_REF},
        "delta_T_K": DELTA_T,
        "kappa_1_per_m": KAPPA,
        "runs": runs,
        "tip_error_pct": finest["tip_error_pct"],
        "temp_uniformity_max_deviation_K": finest["temp_uniformity_max_deviation_K"],
        "pass": bool(finest["tip_error_pct"] < 5.0),
    }
    (CASE_DIR / "results.json").write_text(json.dumps(summary, indent=2))
    print(json.dumps(summary, indent=2), flush=True)


if __name__ == "__main__":
    main()