r"""
Free-free axial rod — natural frequencies
=========================================

Companion to
:ref:`sphx_glr_gallery_verification_example_verify_axial_rod_nf.py`
(fixed-free).  A uniform rod with both ends free vibrates
longitudinally with cosine mode shapes
:math:`u_n(x) = \cos(n \pi x / L)` and the canonical
free-free natural frequencies

.. math::

    f_n = \frac{n}{2 L} \sqrt{\frac{E}{\rho}},
    \qquad n = 1, 2, 3, \ldots

The mode at :math:`n = 0` is the rigid-body axial translation
— a zero-eigenvalue mode produced by the unconstrained boundary
condition.  Every elastic frequency :math:`f_n, n \ge 1` follows
the formula above.

Implementation
--------------

Drives the existing
:class:`~femorph_solver.validation.problems.FreeFreeRodModes`
problem class on a 40-element TRUSS2 mesh — the cleanest 1D
discretisation of the axial-wave operator
(:doc:`/reference/elements/truss2`).  Transverse DOFs at every
node are pinned to UY = UZ = 0 so the element's zero transverse
stiffness doesn't admit arbitrary side-sway modes; the axial DOF
is unconstrained at both ends, producing the rigid-body /
elastic spectrum the closed form predicts.

The extraction filter walks the spectrum starting at 1 Hz to
skip the (analytically) zero-frequency rigid-body mode that the
solver returns at machine-precision noise level.

References
----------

* Rao, S. S. (2017) *Mechanical Vibrations*, 6th ed., Pearson,
  §8.2 Table 8.1 (axial wave equation eigenvalues).
* Meirovitch, L. (2010) *Fundamentals of Vibrations*, Long Grove,
  §6.6 (free-free rod).
* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002)
  *Concepts and Applications of Finite Element Analysis*, 4th
  ed., Wiley, §11 (modal analysis with rigid-body modes).

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

======================================  =====================  ============================================
Source                                   Reported f_1 [Hz]      Problem ID / location
======================================  =====================  ============================================
Closed form (wave equation)              2 523.77               Rao 2017 §8.2 Table 8.1
Meirovitch (2010) §6.6                   2 523.77               free-free rod
MAPDL Verification Manual                ≈ 2 524                VM-47 (free-free torsion analogue)
Abaqus Verification Manual               ≈ 2 524                AVM 1.6.x free-free rod NF family
======================================  =====================  ============================================
"""

from __future__ import annotations

import numpy as np
import pyvista as pv

from femorph_solver.validation.problems import FreeFreeRodModes

# %%
# Build the model from the validation problem class
# --------------------------------------------------

problem = FreeFreeRodModes()
m = problem.build_model()
print(
    f"free-free rod mesh: {m.grid.n_points} nodes, {m.grid.n_cells} TRUSS2 cells; "
    f"L = {problem.L} m, E = {problem.E / 1e9:.0f} GPa, ρ = {problem.rho:.1f} kg/m³"
)

c = np.sqrt(problem.E / problem.rho)
f1_published = c / (2.0 * problem.L)
f2_published = 2.0 * f1_published
print()
print(f"Wave speed c = √(E/ρ)        = {c:.3f} m/s")
print(f"f_1 = c / (2 L) (published)  = {f1_published:.3f} Hz")
print(f"f_2 = 2 f_1     (published)  = {f2_published:.3f} Hz")

# %%
# Modal solve + frequency extraction
# ----------------------------------
#
# A free-free rod exhibits one rigid-body axial-translation mode
# at :math:`f = 0` (machine-precision noise after numeric
# discretisation).  ``problem.extract`` skips it via a 1 Hz floor;
# the next two finite frequencies are the analytical :math:`f_1,
# f_2`.

res = m.solve_modal(n_modes=problem.default_n_modes if hasattr(problem, "default_n_modes") else 6)
freqs = np.asarray(res.frequency, dtype=np.float64)

print()
print("  All recovered frequencies (lowest first):")
for k, f in enumerate(freqs):
    label = "rigid-body" if f < 1.0 else "elastic"
    print(f"    mode {k + 1}: f = {f:>10.4f} Hz   ({label})")

f1_computed = problem.extract(m, res, "f1_axial")
f2_computed = problem.extract(m, res, "f2_axial")
err1 = (f1_computed - f1_published) / f1_published * 100.0
err2 = (f2_computed - f2_published) / f2_published * 100.0
print()
print(
    f"  f_1 computed = {f1_computed:.3f} Hz   published = {f1_published:.3f} Hz   err = {err1:+.4f} %"
)
print(
    f"  f_2 computed = {f2_computed:.3f} Hz   published = {f2_published:.3f} Hz   err = {err2:+.4f} %"
)

assert abs(err1) < 2.0, f"f_1 deviation {err1:.3f}% exceeds 2 %"
assert abs(err2) < 5.0, f"f_2 deviation {err2:.3f}% exceeds 5 %"

# %%
# Render the first elastic axial mode
# -----------------------------------
#
# The first elastic mode has :math:`u(x) = \cos(\pi x / L)` —
# antisymmetric about the rod's midspan, with the two ends moving
# in opposite directions and a node (zero-displacement) at
# :math:`x = L/2`.

shapes = np.asarray(res.mode_shapes).reshape(-1, 3, len(freqs))
elastic_idx = next((k for k, f in enumerate(freqs) if f > 1.0), 0)
phi = shapes[:, :, elastic_idx]

scale = 0.05 * problem.L / max(np.abs(phi).max(), 1e-30)
warped = m.grid.copy()
warped.points = m.grid.points + scale * phi
warped["|phi|"] = np.linalg.norm(phi, axis=1)
warped["ux"] = phi[:, 0]

plotter = pv.Plotter(off_screen=True, window_size=(720, 280))
plotter.add_mesh(m.grid, style="wireframe", color="grey", opacity=0.3)
plotter.add_mesh(warped, scalars="ux", cmap="coolwarm", line_width=4, show_edges=False)
plotter.add_text(
    f"first elastic mode  —  f = {freqs[elastic_idx]:.2f} Hz",
    position="upper_left",
    font_size=11,
)
plotter.view_xy()
plotter.camera.zoom(1.05)
plotter.show()