{
  "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 with off-center point load\n\nA prismatic beam of length $L$ is simply supported at both\nends and loaded by a single concentrated transverse force\n$P$ at $x = a$ (with $b = L - a$).  The\nEuler-Bernoulli closed form (Timoshenko 1955 \u00a728; Roark Table 8\ncase 2) gives reactions\n\n\\begin{align}R_{\\mathrm{left}} = \\frac{P\\, b}{L},\n    \\qquad\n    R_{\\mathrm{right}} = \\frac{P\\, a}{L},\\end{align}\n\na peak bending moment under the load\n\n\\begin{align}M_{\\mathrm{load}} = \\frac{P\\, a\\, b}{L},\\end{align}\n\nand a deflection profile that splits at $x = a$:\n\n\\begin{align}:label: ss-beam-off-center-deflection\n\n    v(x) = \\begin{cases}\n      \\dfrac{P\\, b\\, x\\,(L^{2} - b^{2} - x^{2})}{6\\, L\\, E\\, I},\n        & 0 \\le x \\le a \\\\[8pt]\n      \\dfrac{P\\, a\\,(L - x)\\,(2 L\\, x - a^{2} - x^{2})}{6\\, L\\, E\\, I},\n        & a \\le x \\le L.\n    \\end{cases}\\end{align}\n\nSo the deflection at the load and the absolute maximum deflection\nfollow as\n\n\\begin{align}v(a)\n    \\;=\\;\n    \\frac{P\\, a^{2}\\, b^{2}}{3\\, E\\, I\\, L},\n    \\qquad\n    v_{\\max}\n    \\;=\\;\n    \\frac{P\\, b\\, (L^{2} - b^{2})^{3/2}}{9\\sqrt{3}\\, E\\, I\\, L},\\end{align}\n\nwith the maximum occurring in the **longer** span at\n$x_{\\max} = \\sqrt{(L^{2} - b^{2}) / 3}$ from the closer\nsupport (Roark Table 8 case 2; both formulas valid for\n$a \\le L/2$).\n\nFor the centred case $a = L/2$ the formula collapses to\n$v(L/2) = P L^{3} / (48\\, E\\, I)$, recovering the\n`sphx_glr_gallery_verification_example_verify_ss_beam_central_load.py`\nresult.  Off-centring the load reduces both $v_{\\max}$ and\n$M_{\\max}$ because the moment-arm distribution changes \u2014 at\n$a = L / 3$, the peak moment drops to\n$2 P L / 9$ (vs $P L / 4$ for the centred case) and\nthe peak deflection falls to roughly 88 % of the centred value.\n\n## Implementation\n\nA 60-element BEAM2 (Hermite-cubic Bernoulli) line spans the beam.\nThe two ends are simply supported (only $u_y = 0$); the load\nis applied at the closest mesh node to $x = a$.  Hermite-cubic\nshape functions interpolate Bernoulli kinematics exactly on the\npiecewise-linear moment field this point load creates, so the FE\ndeflection at the loaded node and at the analytical peak position\nmatch the closed form to machine precision.\n\n## References\n\n* Timoshenko, S. P. (1955) *Strength of Materials, Part I:\n  Elementary Theory and Problems*, 3rd ed., Van Nostrand, \u00a728\n  (simply-supported beam, concentrated load at any point).\n* Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for\n  Stress and Strain*, 6th ed., McGraw-Hill, Table 8 case 2.\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J.\n  (2002) *Concepts and Applications of Finite Element\n  Analysis*, 4th ed., Wiley, \u00a72.5 (Hermite-cubic shape\n  functions).\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\nimport femorph_solver\nfrom femorph_solver import ELEMENTS"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Problem data \u2014 point load at L / 3\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "E = 2.0e11\nNU = 0.30\nRHO = 7850.0\nWIDTH = 0.05\nHEIGHT = 0.05\nA_section = WIDTH * HEIGHT\nI_z = WIDTH * HEIGHT**3 / 12.0\nI_y = HEIGHT * WIDTH**3 / 12.0\nJ = (1.0 / 3.0) * min(WIDTH, HEIGHT) ** 3 * max(WIDTH, HEIGHT)\n\nL = 1.0  # span [m]\na_load = L / 3.0  # load position from the left support\nb_load = L - a_load\nP = 1.0e3  # transverse load [N, downward = -y]\n\nEI = E * I_z\n\n# Closed-form references (Roark Table 8 case 2).\nR_left_pub = P * b_load / L\nR_right_pub = P * a_load / L\nM_load_pub = P * a_load * b_load / L\nv_load_pub = -P * a_load**2 * b_load**2 / (3.0 * EI * L)\n\n# Peak deflection: in the longer span (from a_load to L since a < L/2).\n# x_peak measured FROM the right support; convert to absolute x.\nx_peak_from_right = np.sqrt((L**2 - b_load**2) / 3.0)\nx_peak_abs = L - x_peak_from_right\nv_peak_pub = -P * b_load * (L**2 - b_load**2) ** 1.5 / (9.0 * np.sqrt(3.0) * EI * L)\n\nprint(\"Simply-supported beam, off-center point load\")\nprint(f\"  L = {L} m, a = L/3 = {a_load:.4f} m, b = 2L/3 = {b_load:.4f} m\")\nprint(f\"  P = {P} N (downward), E I = {EI:.3e} N m\u00b2\")\nprint()\nprint(\"Closed-form references (Roark Table 8 case 2):\")\nprint(f\"  R_left  = {R_left_pub:+.4e} N      (= P b / L = 2P/3)\")\nprint(f\"  R_right = {R_right_pub:+.4e} N      (= P a / L = P/3)\")\nprint(f\"  M_load  = {M_load_pub:+.4e} N m    (= P a b / L = 2 P L / 9)\")\nprint(f\"  v(a)    = {v_load_pub * 1e3:+.4e} mm   (= P a\u00b2 b\u00b2 / (3 E I L))\")\nprint(f\"  v_max   = {v_peak_pub * 1e3:+.4e} mm   at x = {x_peak_abs:.4f} m\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build a 60-element BEAM2 line\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "N_ELEM = 60\nxs = np.linspace(0.0, L, N_ELEM + 1)\npts = np.column_stack((xs, np.zeros_like(xs), np.zeros_like(xs)))\n\ncells = np.empty((N_ELEM, 3), dtype=np.int64)\ncells[:, 0] = 2\ncells[:, 1] = np.arange(N_ELEM)\ncells[:, 2] = np.arange(1, N_ELEM + 1)\ngrid = pv.UnstructuredGrid(\n    cells.ravel(),\n    np.full(N_ELEM, 3, dtype=np.int64),\n    pts,\n)\n\nm = femorph_solver.Model.from_grid(grid)\nm.assign(\n    ELEMENTS.BEAM2,\n    material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n    real=(A_section, I_z, I_y, J),\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Boundary conditions\n\nSimply supported at both ends \u2014 pin UY only.  Suppress out-of-\nplane motion + axial rigid body across the line; pin UX at the\nleft support to remove the axial rigid-body mode.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m.fix(nodes=1, dof=\"UY\")\nm.fix(nodes=int(N_ELEM + 1), dof=\"UY\")\nm.fix(nodes=1, dof=\"UX\")\n\nfor i in range(N_ELEM + 1):\n    m.fix(nodes=int(i + 1), dof=\"UZ\")\n    m.fix(nodes=int(i + 1), dof=\"ROTX\")\n    m.fix(nodes=int(i + 1), dof=\"ROTY\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Apply the off-center point load\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "i_load = int(np.argmin(np.abs(xs - a_load)))\nm.apply_force(int(i_load + 1), fy=-P)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Solve + extract\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()\ndof_map = m.dof_map()\ndisp = np.asarray(res.displacement, dtype=np.float64)\n\n\ndef _node_dof(node_id: int, dof_idx: int) -> float:\n    \"\"\"0=UX, 1=UY, 2=UZ, 3=ROTX, 4=ROTY, 5=ROTZ.\"\"\"\n    rows = np.where(dof_map[:, 0] == node_id)[0]\n    for r in rows:\n        if int(dof_map[r, 1]) == dof_idx:\n            return float(disp[r])\n    return 0.0\n\n\nv_load_fe = _node_dof(int(i_load + 1), 1)\n\n\ndef v_closed_form(x: float) -> float:\n    \"\"\"Piecewise-cubic deflection profile (Roark Table 8 case 2).\"\"\"\n    if x <= a_load:\n        return -P * b_load * x * (L**2 - b_load**2 - x**2) / (6.0 * L * EI)\n    return -P * a_load * (L - x) * (2.0 * L * x - a_load**2 - x**2) / (6.0 * L * EI)\n\n\n# The Bernoulli-Hermite-cubic FE response is exact at every node\n# for this problem.  The analytical absolute peak at\n# ``x_peak_abs`` falls between mesh nodes \u2014 comparing FE vs that\n# would conflate FE accuracy with mesh-sampling resolution.\n# Cleaner check: at the FE peak node, the FE value should match\n# the closed form *evaluated at the same x* to machine precision.\n\nfe_uy_all = np.array([_node_dof(int(i + 1), 1) for i in range(N_ELEM + 1)])\ni_peak = int(np.argmax(np.abs(fe_uy_all)))\nv_peak_fe = float(fe_uy_all[i_peak])\nv_at_peak_node_pub = v_closed_form(float(xs[i_peak]))\n\nerr_load = (abs(v_load_fe) - abs(v_load_pub)) / abs(v_load_pub) * 100.0\nerr_peak_node = (abs(v_peak_fe) - abs(v_at_peak_node_pub)) / abs(v_at_peak_node_pub) * 100.0\n\nprint()\nprint(f\"{'quantity':<28}  {'FE':>14}  {'published':>14}  {'err':>9}\")\nprint(f\"{'-' * 28:<28}  {'-' * 14:>14}  {'-' * 14:>14}  {'-' * 9:>9}\")\nprint(\n    f\"{'v(a) at load':<28}  {v_load_fe * 1e3:>10.4f} mm  \"\n    f\"{v_load_pub * 1e3:>10.4f} mm  {err_load:>+8.4f}%\"\n)\nprint(\n    f\"{'v at FE peak node':<28}  {v_peak_fe * 1e3:>10.4f} mm  \"\n    f\"{v_at_peak_node_pub * 1e3:>10.4f} mm  {err_peak_node:>+8.4f}%\"\n)\nprint(f\"  FE peak-node x      = {xs[i_peak]:.4f} m  vs  analytical peak x = {x_peak_abs:.4f} m\")\nprint(\n    f\"  analytical |v_max|  = {abs(v_peak_pub) * 1e3:.4f} mm  \"\n    f\"(off-grid; FE-peak-node value above is on the same cubic curve)\"\n)\n\nassert abs(err_load) < 0.1, f\"v(a) deviation {err_load:.4f} % exceeds 0.1 %\"\nassert abs(err_peak_node) < 0.1, f\"v at FE peak node deviation {err_peak_node:.4f} % exceeds 0.1 %\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot the deflection profile against the closed form\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fe_uy = fe_uy_all\n\nimport matplotlib.pyplot as plt  # noqa: E402\n\nfig, ax = plt.subplots(1, 1, figsize=(7.5, 4.0))\nxs_dense = np.linspace(0.0, L, 401)\nax.plot(\n    xs_dense,\n    [v_closed_form(x) * 1e3 for x in xs_dense],\n    \"-\",\n    color=\"#1f77b4\",\n    lw=2.0,\n    label=\"closed form (Roark Table 8 case 2)\",\n)\nax.plot(xs, fe_uy * 1e3, \"o\", color=\"#d62728\", markersize=4, label=\"BEAM2 FE\")\nax.axvline(a_load, color=\"grey\", ls=\"--\", lw=1.0, label=f\"load at x = a = {a_load:.4f} m\")\nax.axvline(x_peak_abs, color=\"black\", ls=\":\", lw=1.0, label=f\"v_max at x = {x_peak_abs:.4f} m\")\nax.axhline(0.0, color=\"black\", lw=0.5)\nax.set_xlabel(\"x [m]\")\nax.set_ylabel(\"v(x)  [mm]   (downward negative)\")\nax.set_title(\"SS beam, off-center point load \u2014 FE vs closed form\", fontsize=11)\nax.legend(loc=\"lower left\", fontsize=9)\nax.grid(True, ls=\":\", alpha=0.5)\nfig.tight_layout()\nfig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the deformed beam\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g = m.grid.copy()\ndisp_xyz = np.zeros((g.n_points, 3))\nfor ni in range(g.n_points):\n    disp_xyz[ni, 1] = _node_dof(int(ni + 1), 1)\ng.point_data[\"displacement\"] = disp_xyz\ng.point_data[\"UY\"] = disp_xyz[:, 1]\n\nscale = 1.0 / max(abs(v_peak_fe), 1e-12) * 0.05\nwarped = g.copy()\nwarped.points = np.asarray(g.points) + scale * disp_xyz\n\nplotter = pv.Plotter(off_screen=True, window_size=(900, 320))\nplotter.add_mesh(g, color=\"grey\", opacity=0.4, line_width=2, label=\"undeformed\")\nplotter.add_mesh(warped, scalars=\"UY\", cmap=\"coolwarm\", line_width=4)\nsupport_pts = np.array([[0.0, 0.0, 0.0], [L, 0.0, 0.0]])\nplotter.add_points(\n    support_pts,\n    render_points_as_spheres=True,\n    point_size=18,\n    color=\"black\",\n    label=\"simple supports\",\n)\nplotter.add_points(\n    np.array([[a_load, scale * v_load_fe, 0.0]]),\n    render_points_as_spheres=True,\n    point_size=14,\n    color=\"#d62728\",\n    label=f\"load (P = {P:.0f} N at x = {a_load:.3f} m)\",\n)\nplotter.add_points(\n    np.array([[xs[i_peak], scale * v_peak_fe, 0.0]]),\n    render_points_as_spheres=True,\n    point_size=12,\n    color=\"#9467bd\",\n    label=f\"v_max FE = {v_peak_fe * 1e3:.4f} mm\",\n)\nplotter.view_xy()\nplotter.add_legend()\nplotter.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Take-aways\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print()\nprint(\"Take-aways:\")\nprint(f\"  \u2022 v(a) at the load within {abs(err_load):.4f} % of P a\u00b2 b\u00b2 / (3 E I L).\")\nprint(\n    f\"  \u2022 v at the FE peak node within {abs(err_peak_node):.4f} % of the closed-form \"\n    \"deflection at the same x \u2014 Hermite-cubic FE is exact at every node for this load.\"\n)\nprint(\n    \"  \u2022 Off-centering the load (a = L/3) reduces both peak moment \"\n    \"(2PL/9 vs PL/4 centred) and peak deflection (~88 % of centred) \u2014 \"\n    \"the moment-arm asymmetry dominates.\"\n)"
      ]
    }
  ],
  "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
}