{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import pyvista\npyvista.OFF_SCREEN = True\npyvista.set_jupyter_backend('static')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Simply-supported beam \u2014 first three bending natural frequencies\n\nFree-vibration verification for a slender simply-supported beam.\nThe transverse-bending eigenvalue problem on the Euler-Bernoulli\nslender-beam operator has a closed-form spectrum\n\n\\begin{align}f_n = \\frac{n^{2}\\, \\pi}{2 L^{2}}\n          \\sqrt{\\frac{E\\, I}{\\rho\\, A}},\n    \\qquad n = 1, 2, 3, \\ldots\\end{align}\n\nwith mode shapes $\\phi_n(x) = \\sin(n \\pi x / L)$.  The\nn-th mode has $n$ antinodes along the span and\n$n - 1$ interior nodes (zero-crossings).\n\nCompanion to\n`sphx_glr_gallery_verification_example_verify_cantilever_higher_modes.py`\n(clamped-free) and\n`sphx_glr_gallery_verification_example_verify_cc_beam_modes.py`\n(clamped-clamped) \u2014 same physics, different end conditions, and\neach end-condition family produces a distinct eigenvalue\ncharacteristic equation.\n\n## Implementation\n\nDrives the existing\n:class:`~femorph_solver.validation.problems.SimplySupportedBeamModes`\nproblem class on its default 80 \u00d7 3 \u00d7 3 HEX8-EAS mesh at\nslenderness $h / L = 1/80$.  At that aspect ratio the\nTimoshenko shear correction the 3D solid truthfully captures\nstays under 0.5 % at mode 3, well below the 0.5 % tolerance the\nbenchmark asserts.\n\nThe mode-selector logic in\n:meth:`~femorph_solver.validation.problems.SimplySupportedBeamModes.extract`\nwalks the spectrum picking successive transverse-dominant (UZ)\nmodes whose antinode counts match $n = 1, 2, 3$ \u2014 a robust\nfilter against bi-degenerate transverse pairs and stray axial /\ntorsional modes.\n\n## References\n\n* Rao, S. S. (2017) *Mechanical Vibrations*, 6th ed., Pearson,\n  \u00a78.5 (transverse-bending eigenvalue problem).\n* Timoshenko, S. P. (1974) *Vibration Problems in Engineering*,\n  4th ed., Wiley, \u00a75.3 (simply-supported beam).\n* Meirovitch, L. (2010) *Fundamentals of Vibrations*, Long Grove,\n  \u00a77.4.\n\n## Vendor cross-references\n\n======================================  =====================  ============================================\nSource                                   Reported f_1 [Hz]      Problem ID / location\n======================================  =====================  ============================================\nClosed form (Euler-Bernoulli)            114.44                 Rao \u00a78.5, Timoshenko VPE \u00a75.3\nMeirovitch (2010) \u00a77.4                   114.44                 SS beam transverse vibration\nMAPDL Verification Manual                \u2248 114.4                VM-89 Natural frequencies of a SS beam\nAbaqus Verification Manual               \u2248 114.4                AVM 1.6.x SS beam natural-frequency family\n======================================  =====================  ============================================\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from __future__ import annotations\n\nimport numpy as np\nimport pyvista as pv\n\nfrom femorph_solver.validation.problems import SimplySupportedBeamModes"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build the model from the validation problem class\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "problem = SimplySupportedBeamModes()\nm = problem.build_model()\nprint(\n    f\"SS beam modes mesh: {m.grid.n_points} nodes, {m.grid.n_cells} HEX8 cells; \"\n    f\"L = {problem.L} m, cross = {problem.width} \u00d7 {problem.height} m, \"\n    f\"h/L = {problem.height / problem.L:.4f}\"\n)\n\nI = problem.width * problem.height**3 / 12.0  # noqa: E741\nA = problem.width * problem.height\nbase = np.sqrt(problem.E * I / (problem.rho * A))\nprint()\nprint(\"First three transverse-bending natural frequencies (closed form):\")\nfor n in (1, 2, 3):\n    fn = (n * n * np.pi) / (2.0 * problem.L**2) * base\n    print(f\"  f_{n} = {fn:.4f} Hz   (n\u00b2\u03c0 / (2 L\u00b2)\u00b7\u221a(EI/\u03c1A))\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Modal solve + per-mode extraction\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_modal(n_modes=problem.default_n_modes)\nprint()\nprint(f\"  {'mode':>4}  {'f computed [Hz]':>16}  {'f published [Hz]':>17}  {'rel err [%]':>12}\")\nprint(\"  \" + \"-\" * 56)\nfor n in (1, 2, 3):\n    name = f\"f{n}_bending\"\n    f_computed = problem.extract(m, res, name)\n    f_published = (n * n * np.pi) / (2.0 * problem.L**2) * base\n    err = (f_computed - f_published) / f_published * 100.0\n    print(f\"  {n:>4}  {f_computed:>14.4f}    {f_published:>15.4f}    {err:>+10.4f}\")\n    # 1 % tolerance: at h/L = 1/80 the Timoshenko shear correction\n    # the 3D solid truthfully captures rises smoothly with mode\n    # number, sitting at ~0.5 % on mode 3 \u2014 comfortably below the\n    # 1 % gate but above the 0.5 % the asymptotic Bernoulli closed\n    # form would suggest.  The problem class itself uses the same\n    # 1 % tolerance for the same reason.\n    assert abs(err) < 1.0, f\"mode {n} deviation {err:.3f}% exceeds 1 %\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the first three transverse modes\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "shapes = np.asarray(res.mode_shapes).reshape(-1, 3, len(res.frequency))\nfreqs = np.asarray(res.frequency, dtype=np.float64)\n\n# Same antinode-count selector ``problem.extract`` uses, replicated\n# inline for the visualisation so the rendered modes are exactly\n# the ones the assertions just passed.\npts = np.asarray(m.grid.points)\nline_mask = (pts[:, 1] < 1e-9) & (pts[:, 2] > pts[:, 2].max() - 1e-9)\n\n\ndef antinode_count(uz_along_x: np.ndarray) -> int:\n    sign_changes = np.sum(np.diff(np.signbit(uz_along_x)).astype(int))\n    return int(sign_changes + 1)\n\n\nx_line = pts[line_mask, 0]\norder = np.argsort(x_line)\n\npicked = {}\nfor k in range(len(freqs)):\n    u = shapes[:, :, k]\n    uz_frac = (u[:, 2] ** 2).sum() / ((u**2).sum() + 1e-30)\n    if uz_frac < 0.5:\n        continue\n    uz_top = u[line_mask, 2][order]\n    n = antinode_count(uz_top)\n    if n in (1, 2, 3) and n not in picked:\n        picked[n] = (k, uz_top)\n    if len(picked) == 3:\n        break\n\nplotter = pv.Plotter(shape=(1, 3), off_screen=True, window_size=(1300, 280), border=False)\nfor col, n in enumerate((1, 2, 3)):\n    plotter.subplot(0, col)\n    k, _ = picked[n]\n    grid_warped = m.grid.copy()\n    disp = shapes[:, :, k]\n    peak = float(np.linalg.norm(disp, axis=1).max())\n    scale = (0.06 * problem.L) / peak if peak > 0 else 1.0\n    grid_warped.points = np.asarray(m.grid.points) + scale * disp\n    grid_warped[\"uz\"] = disp[:, 2]\n    plotter.add_mesh(m.grid, style=\"wireframe\", color=\"grey\", opacity=0.3)\n    plotter.add_mesh(\n        grid_warped, scalars=\"uz\", cmap=\"coolwarm\", show_edges=False, show_scalar_bar=False\n    )\n    plotter.add_text(f\"mode {n}  \u2014  f = {freqs[k]:.2f} Hz\", position=\"upper_left\", font_size=10)\n    plotter.view_xz()\n    plotter.camera.zoom(1.1)\nplotter.link_views()\nplotter.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}