Source code for femorph_solver.validation._problems.cantilever_nf

"""Clamped cantilever — first transverse bending natural frequency.

The fundamental bending frequency of a clamped-free prismatic
beam (Euler-Bernoulli kinematics) is

.. math::

    f_1 = \\frac{(\\beta_1 L)^2}{2 \\pi L^2}
          \\sqrt{\\frac{E I}{\\rho A}}, \\qquad \\beta_1 L = 1.8751,

where :math:`\\beta_1 L` is the first root of
:math:`1 + \\cos(\\beta L) \\cosh(\\beta L) = 0`.

References
----------
* Rao, S. S.  *Mechanical Vibrations*, 6th ed., Pearson, 2017,
  §8.5, Table 8.1 first entry.
* Timoshenko, S. P.  *Vibration Problems in Engineering*, 4th
  ed., Wiley, 1974, §5.3.
"""

from __future__ import annotations

import math
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,
)

#: First root of :math:`1 + \\cos(\\beta L)\\cosh(\\beta L) = 0`.
BETA_1_L = 1.8751040687119611




@dataclass
class CantileverNaturalFrequency(BenchmarkProblem):
    """Cantilever first bending natural frequency (Rao 2017 Table 8.1)."""

    name: str = "cantilever_nf"
    description: str = (
        "First transverse bending frequency of a clamped-free prismatic "
        "beam — Rao 2017 §8.5 Table 8.1."
    )
    analysis: str = "modal"
    #: Request 10 modes so the extract can filter for the first
    #: *bending* mode even if axial / torsional modes appear below
    #: it.
    default_n_modes: int = 10

    L: float = 1.0
    width: float = 0.05
    height: float = 0.05
    E: float = 2.0e11
    nu: float = 0.3
    rho: float = 7850.0

    @property
    def published_values(self) -> tuple[PublishedValue, ...]:
        A = self.width * self.height
        I = self.width * self.height**3 / 12.0  # noqa: E741
        omega1 = BETA_1_L**2 * math.sqrt(self.E * I / (self.rho * A * self.L**4))
        f1 = omega1 / (2.0 * math.pi)
        return (
            PublishedValue(
                name="f1_hz",
                value=f1,
                unit="Hz",
                source="Rao 2017 §8.5 Table 8.1",
                formula="f1 = (beta1 L)^2 / (2 pi L^2) sqrt(E I / (rho A))",
                # HEX8 enhanced-strain hits f1 to within ~0.3 % on the
                # default 40 × 3 × 3 mesh.
                tolerance=0.01,
            ),
        )



    def build_model(self, **mesh_params: Any) -> femorph_solver.Model:
        nx = int(mesh_params.get("nx", 40))
        ny = int(mesh_params.get("ny", 3))
        nz = int(mesh_params.get("nz", 3))

        xs = np.linspace(0.0, self.L, nx + 1)
        ys = np.linspace(0.0, self.width, ny + 1)
        zs = np.linspace(0.0, self.height, nz + 1)
        grid = pv.StructuredGrid(
            *np.meshgrid(xs, ys, zs, indexing="ij")
        ).cast_to_unstructured_grid()

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

        pts = np.asarray(m.grid.points)
        clamped = np.where(pts[:, 0] < 1e-9)[0]
        m.fix(nodes=(clamped + 1).tolist(), dof="ALL")

        # Stash for extract().
        m._bench_axis_x = pts[:, 0]  # type: ignore[attr-defined]
        return m




    def extract(self, model: femorph_solver.Model, result: Any, name: str) -> float:
        if name != "f1_hz":
            raise KeyError(f"unknown quantity {name!r}")
        freqs = np.asarray(result.frequency, dtype=np.float64)
        shapes = np.asarray(result.mode_shapes)  # (n_dof, n_modes)
        shapes = shapes.reshape(-1, 3, shapes.shape[-1])  # (n_nodes, 3, n_modes)

        # Pick the first mode dominated by **transverse bending**
        # (non-axial, non-torsional).  The square cross-section
        # (width == height) means the UY and UZ bending modes are
        # degenerate and the eigensolver may return either
        # orientation at the same frequency — accept either.
        n_modes = shapes.shape[-1]
        for k in range(n_modes):
            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()))
            transverse_rms = max(uy_rms, uz_rms)
            if transverse_rms > 2.0 * ux_rms:
                return float(freqs[k])
        # Fallback — first mode (caller gets a loose match and the
        # test will flag it as non-bending).
        return float(freqs[0])