Source code for femorph_solver.validation._problems.free_free_rod_nf

"""Free-free axial-rod natural frequencies.

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

.. math::

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

The :math:`n = 0` mode is the rigid-body axial translation at
zero frequency.  The first elastic mode (:math:`n = 1`) is at
:math:`f_1 = \\sqrt{E/\\rho} / (2 L)` — twice the fixed-free
fundamental.

Useful as a sanity check for the eigen solver's handling of
near-zero / rigid-body modes: the solver should return one
near-zero rigid-body mode followed by the elastic family.

References
----------
* Rao, S. S.  *Mechanical Vibrations*, 6th ed., Pearson, 2017,
  §8.2 Table 8.1 — free-free rod axial modes.
* Meirovitch, L.  *Fundamentals of Vibrations*, Waveland, 2010,
  §6.6 — axial rod with mixed boundary conditions.

Cross-references (public verification manuals — fair-use
problem-ID citations only; no vendor content vendored):

* **Abaqus AVM 1.6.x** free-free-rod NF family.
"""

from __future__ import annotations

from dataclasses import dataclass
from typing import Any

import numpy as np
import pyvista as pv

import femorph_solver
from femorph_solver import ELEMENTS
from femorph_solver.validation._benchmark import (
    BenchmarkProblem,
    PublishedValue,
)




@dataclass
class FreeFreeRodModes(BenchmarkProblem):
    """Free-free axial rod first elastic + second elastic frequencies."""

    name: str = "free_free_rod_modes"
    description: str = (
        "Free-free uniform rod axial modes — f_n = n / (2L) sqrt(E/ρ) (Rao 2017 §8.2 Table 8.1)."
    )
    analysis: str = "modal"
    default_n_modes: int = 4

    L: float = 1.0
    area: float = 1.0e-4
    E: float = 2.0e11
    nu: float = 0.3
    rho: float = 7850.0

    @property
    def published_values(self) -> tuple[PublishedValue, ...]:
        c = np.sqrt(self.E / self.rho)
        f1 = c / (2.0 * self.L)
        f2 = 2.0 * f1
        return (
            PublishedValue(
                name="f1_axial",
                value=f1,
                unit="Hz",
                source="Rao 2017 §8.2 Table 8.1",
                formula="f1 = sqrt(E/ρ) / (2 L)",
                tolerance=0.02,
            ),
            PublishedValue(
                name="f2_axial",
                value=f2,
                unit="Hz",
                source="Rao 2017 §8.2 Table 8.1",
                formula="f2 = 2 f1 = sqrt(E/ρ) / L",
                tolerance=0.05,
            ),
        )



    def build_model(self, **mesh_params: Any) -> femorph_solver.Model:
        n_elem = int(mesh_params.get("n_elem", 40))
        xs = np.linspace(0.0, self.L, n_elem + 1)
        pts = np.column_stack((xs, np.zeros_like(xs), np.zeros_like(xs)))
        cells = np.empty((n_elem, 3), dtype=np.int64)
        cells[:, 0] = 2
        cells[:, 1] = np.arange(n_elem)
        cells[:, 2] = np.arange(1, n_elem + 1)
        grid = pv.UnstructuredGrid(
            cells.ravel(),
            np.full(n_elem, 3, dtype=np.int64),
            pts,
        )

        m = femorph_solver.Model.from_grid(grid)
        m.assign(
            ELEMENTS.TRUSS2,
            material={"EX": self.E, "PRXY": self.nu, "DENS": self.rho},
            real=(self.area,),
        )

        # Free-free axially: no axial constraint anywhere.  Pin
        # transverse DOFs (UY, UZ) on every node so the truss
        # element's zero transverse stiffness doesn't admit
        # arbitrary side-sway modes.
        n_nodes = n_elem + 1
        for i in range(1, n_nodes + 1):
            m.fix(nodes=i, dof="UY", value=0.0)
            m.fix(nodes=i, dof="UZ", value=0.0)
        return m




    def extract(self, model: femorph_solver.Model, result: Any, name: str) -> float:
        freqs = np.asarray(result.frequency, dtype=np.float64)
        if name == "f1_axial":
            # First elastic mode = first finite frequency past the
            # rigid-body translation (which sits at f = 0).  Skip
            # any mode below 1 Hz to filter the rigid-body family.
            for f in freqs:
                if f > 1.0:
                    return float(f)
            return float(freqs[0])
        if name == "f2_axial":
            # Second elastic mode.
            elastic = [f for f in freqs if f > 1.0]
            if len(elastic) >= 2:
                return float(elastic[1])
            return float(freqs[-1])
        raise KeyError(f"unknown quantity {name!r}")