r"""
Cantilever natural frequency — Rao 2017 Table 8.1
==================================================

First transverse bending frequency of a clamped-free prismatic
beam:

.. math::

    f_{1} = \frac{(\beta_{1} L)^{2}}{2 \pi L^{2}}
           \sqrt{\frac{E I}{\rho A}}, \qquad \beta_{1} L = 1.8751.

Reference: Rao, S. S. *Mechanical Vibrations*, 6th ed.,
Pearson, 2017, §8.5, Table 8.1 first entry.

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

======================================  =====================  ============================================
Source                                   Reported f_1 [Hz]      Problem ID / location
======================================  =====================  ============================================
Closed form (Euler-Bernoulli)            208.6                  Rao §8.5 Table 8.1, Timoshenko VPE §5.3
NAFEMS Benchmark Tests Linear Elastic    208.6                  NAFEMS FV32 (cantilever modal)
MAPDL Verification Manual                208.6                  VM-55 (thin bar free vibration)
Abaqus Verification Manual               208.6                  AVM 1.4.2 (cantilever beam natural frequencies)
CalculiX ccx/test/beam1b.inp             208.0                  CalculiX 2.21 modal 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 CantileverNaturalFrequency

# %%
# Build + solve + convergence study
# ---------------------------------

problem = CantileverNaturalFrequency()
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()
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':<22s} {'n_dof':>8s} {'f_1 (Hz)':>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:<22s} {r.n_dof:>8d} {r.computed:+14.4f} "
        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 first bending mode shape
# -----------------------------------

m = problem.build_model(nx=40, ny=3, nz=3)
modal = m.solve_modal(n_modes=10)

shapes = np.asarray(modal.mode_shapes).reshape(-1, 3, modal.mode_shapes.shape[-1])
# Re-run the mode filter from the problem's extract to pick the
# actual first-bending mode (the square cross-section makes the
# UY/UZ bending pair degenerate).
selected = None
for k in range(shapes.shape[-1]):
    ux_rms = float(np.sqrt((shapes[:, 0, k] ** 2).mean()))
    uy_rms = float(np.sqrt((shapes[:, 1, k] ** 2).mean()))
    uz_rms = float(np.sqrt((shapes[:, 2, k] ** 2).mean()))
    if max(uy_rms, uz_rms) > 2.0 * ux_rms:
        selected = k
        break
if selected is None:
    selected = 0

phi = shapes[:, :, selected]
warped = m.grid.copy()
# Normalise the mode shape so the tip amplitude is ~10% of L,
# just for a visible figure.
scale = 0.1 * problem.L / max(np.linalg.norm(phi, axis=1).max(), 1.0e-30)
warped.points = m.grid.points + scale * phi
warped["|phi|"] = np.linalg.norm(phi, 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="|phi|", cmap="plasma", show_edges=False)
plotter.view_isometric()
plotter.camera.zoom(1.1)
plotter.show()

# %%
# Acceptance

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