r"""
Cantilever tip deflection — Euler-Bernoulli closed form
=======================================================

Slender clamped cantilever under a transverse tip load has
tip deflection

.. math::

    \delta = \frac{P L^{3}}{3 E I}, \qquad I = \frac{w h^{3}}{12}.

Reference: Timoshenko, S. P. *Strength of Materials, Part I*,
3rd ed., Van Nostrand, 1955, §5.4.

This gallery example drives a multi-refinement
:class:`~femorph_solver.validation.ConvergenceStudy` and
renders the deformed shape at the finest mesh.

Vendor cross-references
-----------------------

======================================  =====================  ============================================
Source                                   Reported δ [m]         Problem ID / location
======================================  =====================  ============================================
Closed form (Euler-Bernoulli)            3.81 × 10⁻⁵            Timoshenko SoM Part I §5.4
NAFEMS Background to Benchmarks §3.1     3.81 × 10⁻⁵            NAFEMS BtB-3.1 (cantilever)
MAPDL Verification Manual                3.81 × 10⁻⁵            VM-2 (beam stresses and deflections)
Abaqus Verification Manual               3.81 × 10⁻⁵            AVM 1.4.3 (cantilever with end shear)
CalculiX ccx/test/beamp.inp              3.80 × 10⁻⁵            CalculiX 2.21 statics test
======================================  =====================  ============================================
"""

from __future__ import annotations

import numpy as np
import pyvista as pv

from femorph_solver.validation import ConvergenceStudy
from femorph_solver.validation.problems import CantileverEulerBernoulli

# %%
# Problem + convergence ladder
# ----------------------------
#
# Three in-plane refinements (slenderness L/h = 20); three
# transverse layers keeps the shear-strain field resolved at
# the HEX8 EAS formulation used under the hood.

problem = CantileverEulerBernoulli()
study = ConvergenceStudy(
    problem=problem,
    refinements=[
        {"nx": 20, "ny": 3, "nz": 3},
        {"nx": 40, "ny": 3, "nz": 3},
        {"nx": 80, "ny": 3, "nz": 3},
    ],
)
records = study.run()

# %%
# Print the comparison table

rec = records[0]
pub = rec.results[0].published
print(f"{problem.name}: {problem.description}")
print(f"\n  source : {pub.source}")
print(f"  formula: {pub.formula}")
print(f"  published: {pub.value:.6e} {pub.unit}\n")
print(f"  {'mesh':<20s} {'n_dof':>8s} {'computed':>14s} {'rel err':>10s}  pass")
for r in rec.results:
    mesh_str = "  ".join(f"{k}={v}" for k, v in r.mesh_params.items())
    pass_str = "✓" if r.passed else "✗"
    print(
        f"  {mesh_str:<20s} {r.n_dof:>8d} {r.computed:+14.6e} "
        f"{r.rel_error * 100:+9.2f}%   {pass_str}"
    )
if rec.convergence_rate is not None:
    print(f"\n  fitted rate: |err| ∝ n_dof^(-{rec.convergence_rate:.2f})")

# %%
# Render the deformed shape at the finest mesh
# --------------------------------------------

m = problem.build_model(nx=80, ny=3, nz=3)
res = m.solve_static()
u = np.asarray(res.displacement).reshape(-1, 3)

warped = m.grid.copy()
warped.points = m.grid.points + 1.0 * u
warped["|u|"] = np.linalg.norm(u, axis=1)

plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(m.grid, style="wireframe", color="#888888", opacity=0.3)
plotter.add_mesh(warped, scalars="|u|", cmap="viridis", show_edges=False)
plotter.view_isometric()
plotter.camera.zoom(1.1)
plotter.show()

# %%
# Acceptance
# ----------

finest = rec.results[-1]
assert finest.passed, (
    f"finest-mesh tip_deflection failed: {finest.rel_error:.2%} "
    f"above tolerance {finest.published.tolerance:.0%}"
)