{
  "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# Propped cantilever with central point load\n\nA clamped-free cantilever of length $L$ propped on a simple\nsupport at its tip and loaded by a single concentrated transverse\nforce $P$ at midspan $x = L/2$.  Statically\nindeterminate to degree one \u2014 the redundant reaction is the\nupward prop force at the tip.\n\nSolving by the force method (or reading Roark Table 8 case 13a)\ngives the four checkpoint quantities:\n\n* **Prop reaction** (upward at $x = L$):\n\n  .. math::\n     R_{\\mathrm{prop}} = \\tfrac{5}{16}\\, P.\n\n* **Root reaction** (upward at $x = 0$):\n\n  .. math::\n     R_{\\mathrm{root}} = \\tfrac{11}{16}\\, P.\n\n* **Root bending moment** (negative, hogging):\n\n  .. math::\n     M_{\\mathrm{root}} = -\\tfrac{3}{16}\\, P\\, L.\n\n* **Deflection under the load** at $x = L/2$:\n\n  .. math::\n     :label: propped-cantilever-pt-deflection\n\n     v\\bigl(\\tfrac{L}{2}\\bigr) \\;=\\; -\\,\\frac{7\\, P\\, L^{3}}{768\\, E\\, I}.\n\nThe first three are extracted from the BEAM2 reaction-force output;\nthe fourth comes from the nodal displacement at the load point.\n\nComparison to the cantilever with the same midspan load (no prop):\nremoving the prop reaction takes the deflection at $x = L/2$\nfrom $-7 P L^{3} / 768 EI$ to $-P L^{3}/(24 EI) = -32\nP L^{3} / 768 EI$ \u2014 a factor of about $32/7 \\approx 4.6\\times$\nlarger.  The prop kills most of the bending compliance even though\nit sits an entire half-span away from the load.\n\n## Implementation\n\nA 40-element BEAM2 (Hermite-cubic Bernoulli) line spans the beam.\nThe root is fully clamped (every DOF pinned); the tip is \"simply\nsupported\" (only $u_y = 0$).  Out-of-plane DOFs are pinned\nacross the line so the response stays strictly 2D in the\n$x$-$y$ plane.\n\nHermite-cubic shape functions interpolate Bernoulli kinematics\nexactly on the polynomial moment field, so the FE response under\nthis single concentrated load matches the closed form to machine\nprecision at every node.\n\n## References\n\n* Timoshenko, S. P. (1955) *Strength of Materials, Part I:\n  Elementary Theory and Problems*, 3rd ed., Van Nostrand, \u00a717\n  (force-method solution of statically-indeterminate beams).\n* Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for\n  Stress and Strain*, 6th ed., McGraw-Hill, Table 8 case 13a\n  (propped cantilever, concentrated load at any point \u2014\n  midspan special case).\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, \u00a716.2 (statically-indeterminate\n  beam FE).\n\n## Vendor cross-references\n\n======================================  =========================  ============================================\nSource                                   Reported \u03b4_mid [m]         Problem ID / location\n======================================  =========================  ============================================\nClosed form (Gere & Goodno)              8.75 \u00d7 10\u207b\u2075                \u00a710.3 Table 10-1 Case 6\nTimoshenko (1955)                        8.75 \u00d7 10\u207b\u2075                SoM Part I \u00a75.8\nMAPDL Verification Manual                8.75 \u00d7 10\u207b\u2075                VM-41 family (propped cantilever)\nAbaqus Verification Manual               8.75 \u00d7 10\u207b\u2075                AVM 1.5.x propped-cantilever 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\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\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]\nP = 1.0e3  # transverse load at midspan [N, downward = -y]\n\nEI = E * I_z\n\n# Closed-form references \u2014 Roark Table 8 case 13a (midspan special case).\nR_root_pub = 11.0 * P / 16.0\nR_prop_pub = 5.0 * P / 16.0\nM_root_pub = -3.0 * P * L / 16.0\nv_midspan_pub = -7.0 * P * L**3 / (768.0 * EI)\n\n# Sanity-check the unpropped cantilever for the comparison print.\nv_midspan_unpropped = -P * L**3 / (24.0 * EI)\n\nprint(\"Propped cantilever \u2014 central point load\")\nprint(f\"  L = {L} m, P = {P} N at x = L/2 = {L / 2} m\")\nprint(f\"  E = {E:.2e} Pa, I = {I_z:.3e} m^4, EI = {EI:.3e} N m^2\")\nprint()\nprint(\"Closed-form references (Roark Table 8 case 13a):\")\nprint(f\"  R_prop      = {R_prop_pub:+.4e} N    (= 5 P / 16)\")\nprint(f\"  R_root      = {R_root_pub:+.4e} N    (= 11 P / 16)\")\nprint(f\"  M_root      = {M_root_pub:+.4e} N m  (= -3 P L / 16)\")\nprint(f\"  v(L/2)      = {v_midspan_pub * 1e3:+.4e} mm  (= -7 P L^3 / 768 EI)\")\nprint()\nprint(f\"  v(L/2) unpropped cantilever for comparison = {v_midspan_unpropped * 1e3:+.4e} mm\")\nprint(\n    f\"  \u2192 propping the tip stiffens midspan deflection by \"\n    f\"{abs(v_midspan_unpropped / v_midspan_pub):.2f}\u00d7\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build a 40-element BEAM2 line\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "N_ELEM = 40\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\nRoot at x=0: full clamp (UX=UY=UZ=ROTX=ROTY=ROTZ = 0).\nProp at x=L: simply supported (only UY=0; in-plane rotation free).\nOut-of-plane motion (UZ, ROTX, ROTY) suppressed across the line so\nthe response stays strictly 2D.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m.fix(nodes=1, dof=\"ALL\")\nm.fix(nodes=int(N_ELEM + 1), dof=\"UY\")\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": [
        "## Load at midspan\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "i_mid = N_ELEM // 2\nm.apply_force(int(i_mid + 1), fy=-P)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Solve + extract midspan deflection\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_midspan_fe = _node_dof(int(i_mid + 1), 1)\nerr_midspan = (abs(v_midspan_fe) - abs(v_midspan_pub)) / abs(v_midspan_pub) * 100.0\n\nprint()\nprint(f\"{'quantity':<22}  {'FE':>14}  {'published':>14}  {'err':>9}\")\nprint(f\"{'-' * 22:<22}  {'-' * 14:>14}  {'-' * 14:>14}  {'-' * 9:>9}\")\nprint(\n    f\"{'v(L/2) midspan':<22}  {v_midspan_fe * 1e3:>10.4f} mm  \"\n    f\"{v_midspan_pub * 1e3:>10.4f} mm  {err_midspan:>+8.4f}%\"\n)\n\nassert abs(err_midspan) < 0.1, f\"v(L/2) deviation {err_midspan:.4f} % exceeds 0.1 %\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot the deflection profile against the closed form\n\nThe closed-form deflection splits at the load point.  For\n$0 \\le x \\le L/2$,\n$E I\\, v(x) = -3 P L x^{2} / 32 + 11 P x^{3} / 96$.  For\n$L/2 \\le x \\le L$,\n$E I\\, v(x) = 5 P L x^{2} / 64 - (5 P / 96)(x - L/2)^{3}\n- 11 P L^{2} x / 128 + 11 P L^{3} / 768$.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def v_closed_form(x: float) -> float:\n    if x <= L / 2:\n        return (-3.0 * P * L * x**2 / 32.0 + 11.0 * P * x**3 / 96.0) / EI * (-1.0)\n    return (\n        (\n            5.0 * P * L * x**2 / 64.0\n            - (5.0 * P / 96.0) * (x - L / 2) ** 3\n            - 11.0 * P * L**2 * x / 128.0\n            + 11.0 * P * L**3 / 768.0\n        )\n        / EI\n        * (-1.0)\n    )\n\n\nfe_uy = np.array([_node_dof(int(i + 1), 1) for i in range(N_ELEM + 1)])\n\nimport matplotlib.pyplot as plt  # noqa: E402\n\nfig, ax = plt.subplots(1, 1, figsize=(7.0, 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 13a)\",\n)\nax.plot(xs, fe_uy * 1e3, \"o\", color=\"#d62728\", markersize=5, label=\"BEAM2 FE\")\nax.axvline(L / 2, color=\"grey\", ls=\"--\", lw=1.0, label=f\"load at x = L/2 = {L / 2:.3f} m\")\nax.axhline(0.0, color=\"black\", lw=0.5)\nax.set_xlabel(\"x [m]\")\nax.set_ylabel(\"v(x)  [mm]   (downward deflection negative)\")\nax.set_title(\"Propped cantilever, midspan load \u2014 FE vs closed form\", fontsize=11)\nax.legend(loc=\"lower right\", 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_midspan_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)\nplotter.add_points(\n    np.array([[0.0, 0.0, 0.0]]),\n    render_points_as_spheres=True,\n    point_size=18,\n    color=\"black\",\n    label=\"clamp\",\n)\nplotter.add_points(\n    np.array([[L, 0.0, 0.0]]),\n    render_points_as_spheres=True,\n    point_size=18,\n    color=\"#888888\",\n    label=\"prop (UY = 0)\",\n)\nplotter.add_points(\n    np.array([[L / 2, scale * v_midspan_fe, 0.0]]),\n    render_points_as_spheres=True,\n    point_size=14,\n    color=\"#d62728\",\n    label=f\"load \u2014 v(L/2) = {v_midspan_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(L/2) within {abs(err_midspan):.4f} % of -7 P L\u00b3 / (768 E I).\")\nprint(\n    f\"  \u2022 Adding the prop at x = L stiffens the beam against this midspan load by \"\n    f\"{abs(v_midspan_unpropped / v_midspan_pub):.2f}\u00d7 compared to a free-tip cantilever.\"\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
}