{
  "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 under a central point load\n\nClassical statically-determinate first-order beam: simple\nsupports (pin + roller) at $x = 0$ and $x = L$,\nwith a single transverse point load $P$ at midspan.\nSymmetry plus equilibrium gives the elementary closed form\n\n\\begin{align}R_\\mathrm{left} = R_\\mathrm{right} = \\tfrac{P}{2},\n    \\qquad\n    M_\\mathrm{max} = \\tfrac{P L}{4}\\;\\text{at}\\;x = L/2,\\end{align}\n\nwith mid-span and quarter-span deflections (Roark & Young\n1989, Table 8 case 1; Timoshenko 1955 \u00a732)\n\n\\begin{align}\\delta(L/2) = \\frac{P L^{3}}{48\\, E I},\n    \\qquad\n    \\delta(L/4) = \\frac{11\\, P L^{3}}{768\\, E I}.\\end{align}\n\n## Implementation\n\nA 20-element BEAM2 (Hermite-cubic Bernoulli) line spans the\nbeam.  Boundary conditions:\n\n* node 1 (pin): UY = UZ = UX = 0; ROTX, ROTY suppressed; ROTZ\n  free (slope rotates).\n* node $N_e + 1$ (roller): UY = UZ = 0; ROTX, ROTY\n  suppressed; UX and ROTZ free.\n* tip load $-P\\,\\hat y$ at the central node $N_e/2 + 1$.\n\nThe Hermite-cubic basis recovers Euler-Bernoulli kinematics\nexactly under nodal point loads, so the FE solution at every\nnode equals the analytical value to machine precision.\n\n## References\n\n* Timoshenko, S. P. (1955) *Strength of Materials, Part I*,\n  3rd ed., Van Nostrand, \u00a732 (concentrated transverse loads).\n* Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for\n  Stress and Strain*, 6th ed., McGraw-Hill, Table 8 case 1\n  (simply-supported beam, central point load).\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002)\n  *Concepts and Applications of Finite Element Analysis*, 4th\n  ed., Wiley, \u00a72.5.\n* Gere, J. M. and Goodno, B. J. (2018) *Mechanics of Materials*,\n  9th ed., Cengage, \u00a79.3 Table 9-2 case 5.\n\n## Vendor cross-references\n\n======================================  =====================  ============================================\nSource                                   Reported \u03b4_mid [m]     Problem ID / location\n======================================  =====================  ============================================\nClosed form (Euler-Bernoulli)            2.000 \u00d7 10\u207b\u2074           Timoshenko SoM Part I \u00a75.6\nGere & Goodno (2018) Table 9-2           2.000 \u00d7 10\u207b\u2074           SS beam with concentrated mid-load\nMAPDL Verification Manual                2.00 \u00d7 10\u207b\u2074            VM-2 Beam stresses and deflections (SS)\nAbaqus Verification Manual               2.00 \u00d7 10\u207b\u2074            AVM 1.5.x SS beam family\nNAFEMS Background to Benchmarks          2.00 \u00d7 10\u207b\u2074            \u00a73.2 SS beam with central load\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\nfrom vtkmodules.util.vtkConstants import VTK_LINE\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Problem data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "E = 2.0e11  # Pa (steel)\nNU = 0.30\nRHO = 7850.0  # kg/m^3\nb = 0.050  # square cross-section side [m]\nA = b * b\nI = b**4 / 12.0  # noqa: E741\nJ = 2.0 * I\nL = 1.0  # span [m]\nP = 5.0e3  # tip load magnitude [N], applied downward (-y)\n\nN_ELEM = 20  # must be even so a node lands at midspan\n\n# Closed-form quantities -----------------------------------------------\n\nR_left_published = P / 2.0\nR_right_published = P / 2.0\ndelta_mid_published = P * L**3 / (48.0 * E * I)\ndelta_quarter_published = 11.0 * P * L**3 / (768.0 * E * I)\nM_max_published = P * L / 4.0\n\nprint(f\"P = {P:.1f} N, L = {L:.2f} m, EI = {E * I:.3e} N m^2\")\nprint(f\"R_L = R_R = P/2          = {R_left_published:.4f} N\")\nprint(f\"M_max  = P L / 4         = {M_max_published:.4f} N m at midspan\")\nprint(f\"\u03b4(L/2) = P L^3/(48 EI)   = {delta_mid_published:.4e} m\")\nprint(f\"\u03b4(L/4) = 11 P L^3/(768EI) = {delta_quarter_published:.4e} m\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build a 20-element BEAM2 line\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "points = np.array(\n    [[i * L / N_ELEM, 0.0, 0.0] for i in range(N_ELEM + 1)],\n    dtype=np.float64,\n)\ncells_list: list[int] = []\nfor i in range(N_ELEM):\n    cells_list.extend([2, i, i + 1])\ncells = np.asarray(cells_list, dtype=np.int64)\ncell_types = np.full(N_ELEM, VTK_LINE, dtype=np.uint8)\ngrid = pv.UnstructuredGrid(cells, cell_types, points)\n\nm = femorph_solver.Model.from_grid(grid)\nm.assign(\n    ELEMENTS.BEAM2,\n    material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n    real=(A, I, I, J),\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Boundary conditions\n\nStandard simply-supported beam:\n\n* left pin  \u2014 UY, UZ, UX = 0; ROTX, ROTY suppressed; ROTZ free.\n* right roller \u2014 UY, UZ = 0; ROTX, ROTY suppressed; UX & ROTZ free.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Left pin \u2014 pin axial too so the system is statically determinate.\nm.fix(nodes=[1], dof=\"UX\")\nm.fix(nodes=[1], dof=\"UY\")\nm.fix(nodes=[1], dof=\"UZ\")\nm.fix(nodes=[1], dof=\"ROTX\")\nm.fix(nodes=[1], dof=\"ROTY\")\n\n# Right roller \u2014 only UY pinned (vertical reaction); UX free.\nright = N_ELEM + 1\nm.fix(nodes=[right], dof=\"UY\")\nm.fix(nodes=[right], dof=\"UZ\")\nm.fix(nodes=[right], dof=\"ROTX\")\nm.fix(nodes=[right], dof=\"ROTY\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Apply the central point load\n\nA single $-P\\hat y$ force at the midspan node.  No\ndistributed-load consistent vector needed \u2014 the Hermite cubic\nbasis captures the response of a point load on a Bernoulli\nbeam exactly at every node.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mid_node = N_ELEM // 2 + 1\nm.apply_force(mid_node, fy=-P)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Static solve + reaction extraction\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()\ndof = m.dof_map()\n\n\ndef _react(node: int, dof_idx: int) -> float:\n    rows = np.where((dof[:, 0] == node) & (dof[:, 1] == dof_idx))[0]\n    return float(res.reaction[rows[0]]) if len(rows) else 0.0\n\n\nR_left = _react(1, 1)\nR_right = _react(right, 1)\n\nprint()\nprint(\"reactions  (femorph-solver) \u2192 (analytical)\")\nprint(f\"  R_left_UY   = {R_left:+.4f}  \u2192  {+R_left_published:+.4f} N\")\nprint(f\"  R_right_UY  = {R_right:+.4f}  \u2192  {+R_right_published:+.4f} N\")\n\n# Vertical equilibrium and exact agreement.\nassert np.isclose(R_left + R_right, P, rtol=1e-12), \"vertical equilibrium failed\"\nassert np.isclose(R_left, R_left_published, rtol=1e-12), \"left reaction off\"\nassert np.isclose(R_right, R_right_published, rtol=1e-12), \"right reaction off\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Mid-span and quarter-span deflection\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mid_uy = float(res.displacement[np.where((dof[:, 0] == mid_node) & (dof[:, 1] == 1))[0][0]])\nerr_mid = (mid_uy - (-delta_mid_published)) / (-delta_mid_published)\nprint()\nprint(f\"\u03b4(L/2) computed   = {mid_uy:+.4e} m\")\nprint(f\"\u03b4(L/2) published  = {-delta_mid_published:+.4e} m\")\nprint(f\"relative error    = {err_mid * 100:+.6f} %\")\nassert abs(err_mid) < 1e-8, \"mid-span deflection error too large for Bernoulli BEAM2\"\n\nquarter_node = N_ELEM // 4 + 1\nquarter_uy = float(res.displacement[np.where((dof[:, 0] == quarter_node) & (dof[:, 1] == 1))[0][0]])\nerr_quarter = (quarter_uy - (-delta_quarter_published)) / (-delta_quarter_published)\nprint()\nprint(f\"\u03b4(L/4) computed   = {quarter_uy:+.4e} m\")\nprint(f\"\u03b4(L/4) published  = {-delta_quarter_published:+.4e} m\")\nprint(f\"relative error    = {err_quarter * 100:+.6f} %\")\nassert abs(err_quarter) < 1e-8, \"quarter-span deflection error too large for Bernoulli BEAM2\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the deflected line\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pts = np.asarray(grid.points)\ndisp_y = np.zeros(N_ELEM + 1)\nfor i in range(N_ELEM + 1):\n    rows = np.where((dof[:, 0] == i + 1) & (dof[:, 1] == 1))[0]\n    if len(rows):\n        disp_y[i] = float(res.displacement[rows[0]])\nwarped = grid.copy()\nwarped.points = pts + np.column_stack([np.zeros(N_ELEM + 1), 200.0 * disp_y, np.zeros(N_ELEM + 1)])\nwarped[\"uy\"] = disp_y\n\nplotter = pv.Plotter(off_screen=True, window_size=(720, 280))\nplotter.add_mesh(grid, color=\"grey\", line_width=2, opacity=0.5)\nplotter.add_mesh(warped, scalars=\"uy\", line_width=5, cmap=\"viridis\")\nplotter.add_points(\n    np.array([[0.0, 0.0, 0.0], [L, 0.0, 0.0]]),\n    render_points_as_spheres=True,\n    point_size=18,\n    color=\"black\",\n)\nplotter.add_points(\n    np.array([[0.5 * L, 0.0, 0.0]]),\n    render_points_as_spheres=True,\n    point_size=14,\n    color=\"#d62728\",\n    label=f\"P = {P:.0f} N\",\n)\nplotter.view_xy()\nplotter.camera.zoom(1.05)\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
}