"""
TorsoCAE Validation -- contact mechanics case 2: rigid-ish flat-ended
cylindrical punch pressed into an elastic half-space, compared to the
classical Boussinesq/Sneddon flat-punch contact solution.

Geometry: a short cylinder (radius PUNCH_A) whose flat bottom face contacts
the flat top face of a large box (elastic half-space approximation).
  punch: tag 3 (bottom, z=GAP_INIT) -> contact_slave
         tag 2 (top, z=GAP_INIT+PUNCH_H) -> displacement BC
  box:   tag 5 (bottom, z=-BOX_Z) -> fixed;  tag 6 (top, z=0) -> contact_master
         tags 1-4 (sides) -> roller BCs, same half-space approximation as case_01

**Displacement-controlled, not force-controlled** -- same reasoning as
case_01. The punch's top face is displacement-driven; the solver-reported
`contact_normal_forces` is the true measured contact force, not a preset
target.

PRIMARY validation metric: mean uz on the punch's BOTTOM (contact) face --
not the top face, which is contaminated by the punch's own bulk axial
compression under load (this punch is not perfectly rigid: E1/E2 = 1000,
stiff but finite) -- versus the classical rigid-flat-punch closed form
evaluated at F_meas:
    delta(F) = F / (2 * E_star * a),   E_star = E2 / (1 - nu2^2)
E_star here is the SINGLE half-space's modulus (the formula's rigid-punch
assumption), not the two-body Hertz combination 1/E* = (1-nu1^2)/E1 +
(1-nu2^2)/E2 used in case_01's Hertz formula (that combination applies to
a genuinely deformable sphere, not a rigid punch).

Modulus ratio: E1/E2 = 1000 approximates a rigid punch. Measured stiffness
is independent of load level across 15-60 kN at fixed step count (see
results.json), and 1:1 / 10:1 ratio configurations are recorded there too,
showing error growing smoothly toward the production value as the punch
approaches the rigid limit the closed form assumes.

Box depth (BOX_XY=0.4, BOX_Z=1.2, i.e. depth h = 120x the punch radius a)
and far-field mesh resolution both matter for how well a finite box
approximates an infinite half-space:
  - Near-contact mesh refinement (REFINE_SIZE = a/6 -> a/7, small box):
    error moved 15.37% -> 15.17%, i.e. converged to within 0.2 percentage
    points -- the local contact-zone mesh is not the limiting factor.
  - Far-field mesh refinement at the production depth (size_max
    0.02 -> 0.018 m): measured stiffness excess over theory moved
    0.65% -> 0.18%, same sign both times (a finite box can only be
    stiffer than the ideal half-space, never softer) -- confirms
    convergence rather than a coincidental cancellation (cross-checked
    against the raw displacement field itself, not just the summary
    error, in convergence.json).
  - Box depth alone (h/a = 20 -> 60, fixed far-field mesh size) reduced
    the error from 12.48% -> 7.24%, confirming depth to the fixed bottom
    face is the dominant boundary effect. Depths beyond h/a=120 at the
    same absolute far-field element size are under-resolved (see
    convergence.json's depth_sweep_runs note) and are not used to infer
    the half-space limit.

Run: python3 validation/structural/contact_mechanics/case_02_boussinesq_flat_punch/run_validation.py
(single point takes 45-100+ min depending on mesh/box size, on 1 core,
iterative linear solver; the production config here has ~109k nodes.)
"""
import json
import math
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

PUNCH_A, PUNCH_H = 0.01, 0.005        # m, punch radius / height
E1, NU1 = 210e12, 0.3                 # punch, near-rigid (1000x steel)
E2, NU2 = 210e9, 0.3                  # box (steel)
BOX_XY, BOX_Z = 0.40, 1.20            # m, half-space approximation box, h/a=120 (see module docstring)
BOX_SIZE_MAX = 0.018                  # m, far-field mesh size -- convergence-tested, see module docstring
GAP_INIT = 2.0e-6                     # finite initial gap, avoids degenerate zero-gap point contact at step 0

# Rigid-flat-punch-on-elastic-half-space (Boussinesq): the classical closed
# form treats the punch as rigid and uses the SINGLE half-space's E* =
# E/(1-nu^2), not the two-body Hertz E* -- see module docstring.
E_STAR = E2 / (1.0 - NU2 ** 2)
PUNCH_AREA = math.pi * PUNCH_A ** 2

F_DESIGN = 30000.0   # N, used only to pick a displacement magnitude in the right ballpark

REFINE_SIZE = PUNCH_A / 6.0
REFINE_RADIUS = 3.0 * PUNCH_A


def boussinesq_delta(F: float) -> float:
    return F / (2.0 * E_STAR * PUNCH_A)


def run_case() -> dict:
    delta_design = boussinesq_delta(F_DESIGN)
    total_closure = GAP_INIT + delta_design

    session = TorsoCAESession()
    session.set_physics("structural", submodel="contact_mechanics", backend="mfem")
    session.csg_make_shape("cylinder", {"r": PUNCH_A, "h": PUNCH_H, "cx": 0, "cy": 0, "cz": GAP_INIT},
                            name="punch", body_id="body_0")
    session.csg_select("body_0")
    session.mesh(
        mesh_id="mesh_0", algo_id="delaunay", size_factor=0.5, size_min=0, size_max=PUNCH_A / 1.5, order=1, dim=3,
        point_refinements=[{"x": 0.0, "y": 0.0, "z": GAP_INIT, "size": REFINE_SIZE, "radius": REFINE_RADIUS}],
    )
    session.surface(3, scope_id="body_0").bc("contact_slave")

    session.csg_make_shape("box", {"dx": BOX_XY, "dy": BOX_XY, "dz": BOX_Z,
                                    "cx": -BOX_XY / 2, "cy": -BOX_XY / 2, "cz": -BOX_Z},
                            name="box", body_id="body_1")
    session.csg_select("body_1")
    session.mesh(
        mesh_id="mesh_1", algo_id="delaunay", size_factor=0.7, size_min=0, size_max=BOX_SIZE_MAX, order=1, dim=3,
        point_refinements=[{"x": 0.0, "y": 0.0, "z": 0.0, "size": REFINE_SIZE, "radius": REFINE_RADIUS}],
    )
    session.surface(5, scope_id="body_1").bc("fixed")
    session.surface(1, scope_id="body_1").bc("fixed", frame="global", components=['x'])
    session.surface(2, scope_id="body_1").bc("fixed", frame="global", components=['x'])
    session.surface(3, scope_id="body_1").bc("fixed", frame="global", components=['y'])
    session.surface(4, scope_id="body_1").bc("fixed", frame="global", components=['y'])
    session.surface(6, scope_id="body_1").bc("contact_master", viscous_damping=1.0e6)

    session.solid("punch").material(E=E1, nu=NU1)
    session.solid("box").material(E=E2, nu=NU2)
    session.surface(2, scope_id="body_0").bc("displacement", values=[0.0, 0.0, -total_closure])

    session.set_contact_options(proximity_ratio=2.0, contact_model="frictionless", print_level=2)
    session.set_solver_options(algo="iterative", num_steps=6, dt=0.1, n_cores=1, device="cpu")
    t0 = time.time()
    result = session.compute(mesh_ids=["mesh_0", "mesh_1"])
    elapsed = time.time() - t0

    npz = np.load(result["npz_path"], allow_pickle=True)
    coords = npz["coordinates"]
    disp = npz["displacement"]
    r = np.sqrt(coords[:, 0] ** 2 + coords[:, 1] ** 2)

    # Contact-face (bottom of punch) uz only -- radius-restricted to the
    # punch footprint (r < a, not clipped inside the true edge, where the
    # Boussinesq pressure singularity lives) so box nodes near z=0 can't
    # pollute the mean.
    bottom_mask = (r < PUNCH_A + 1e-9) & (np.abs(coords[:, 2] - GAP_INIT) < 1e-7)
    delta_contact_face = -float(disp[bottom_mask, 2].mean()) - GAP_INIT

    F_meas = float(npz["contact_normal_forces"][0])
    delta_theory_at_Fmeas = boussinesq_delta(F_meas)
    err_pct = abs(delta_contact_face - delta_theory_at_Fmeas) / delta_theory_at_Fmeas * 100.0

    return {
        "F_meas_N": F_meas,
        "delta_contact_face_m": delta_contact_face,
        "delta_theory_at_Fmeas_m": delta_theory_at_Fmeas,
        "err_pct": err_pct,
        "elapsed_s": elapsed,
        "npz_path": result["npz_path"],
    }


def main():
    row = run_case()
    summary = {
        "case": "contact_mechanics / Boussinesq flat cylindrical punch on elastic half-space",
        "material": {"E1_Pa": E1, "nu1": NU1, "E2_Pa": E2, "nu2": NU2, "E_star_Pa": E_STAR},
        "geometry": {"punch_radius_a_m": PUNCH_A, "box_xy_m": BOX_XY, "box_z_m": BOX_Z,
                     "refine_size_m": REFINE_SIZE, "refine_radius_m": REFINE_RADIUS},
        "method": "Displacement-controlled punch top; F_meas is the solver-reported contact "
                  "force, not a preset target. delta_contact_face (mean uz on the punch's "
                  "bottom/contact face, not the top -- see module docstring) is compared to "
                  "the rigid-flat-punch closed form evaluated at F_meas.",
        "point": row,
        "note": "This is a single-run reproduction at the validated configuration. The "
                "checked-in results.json also carries mesh-refinement and domain-size "
                "sensitivity studies (convergence.json, plots/) -- running this script does "
                "not overwrite that unless you intend to update the record too.",
    }
    print(json.dumps(summary, indent=2), flush=True)


if __name__ == "__main__":
    main()