{
  "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# Cantilever with off-tip point load \u2014 load-position shape function\n\nA clamped-free prismatic beam of length $L$ with a single\ntransverse point load $P$ applied at an interior position\n$x = a \\le L$.  The Euler-Bernoulli closed form (Timoshenko\n1955 \u00a732; Roark Table 8 case 1) splits at the load:\n\n\\begin{align}v(x) = \\begin{cases}\n        \\dfrac{P\\, x^{2}\\,(3 a - x)}{6\\, E\\, I}, & 0 \\le x \\le a \\\\[8pt]\n        \\dfrac{P\\, a^{2}\\,(3 x - a)}{6\\, E\\, I}, & a \\le x \\le L\n    \\end{cases}\\end{align}\n\nso the deflection grows as $x^{2}\\,(3 a - x)$ while the load\nis \"above\" $x$ and as $a^{2}\\,(3 x - a)$ past the load\n(a straight ramp of slope $P\\,a^{2} / (2\\, E\\, I)$ away from\nthe load).  Three closed-form check points fall out:\n\n* **Deflection at the load**\n\n  .. math::\n     v(a) = \\frac{P\\, a^{3}}{3\\, E\\, I}.\n\n* **Tip deflection** $v(L)$\n\n  .. math::\n     v(L) = \\frac{P\\, a^{2}\\,(3 L - a)}{6\\, E\\, I}.\n\n* **Tip slope** $\\theta(L)$ (= the slope past the load)\n\n  .. math::\n     \\theta(L) = \\frac{P\\, a^{2}}{2\\, E\\, I}.\n\nFor the tip-load limit $a = L$ the formula collapses to\n$v(L) = P L^{3} / (3 E I)$, recovering the\n`sphx_glr_gallery_verification_example_verify_cantilever_eb.py`\nresult.  For $a = L/2$ (midspan loading) the tip deflection\nis $5 P L^{3} / (48 E I) = 5/16$ of the tip-loaded value \u2014\nload position matters far more than load magnitude on a slender\ncantilever.\n\n## Implementation\n\nA 40-element BEAM2 (Hermite-cubic Bernoulli) line spans the beam.\nThe load is applied at the closest mesh node to $x = a$;\nthe example reports the FE deflection / slope at three check\npoints along with the closed-form expectation for each.\n\nHermite-cubic shape functions interpolate Bernoulli kinematics\nexactly on the polynomial moment field, so the FE response under\nthis point load matches the closed form to machine precision at\nthe loaded node and the tip \u2014 well within a 0.1 % tolerance.\n\n## References\n\n* Timoshenko, S. P. (1955) *Strength of Materials, Part I:\n  Elementary Theory and Problems*, 3rd ed., Van Nostrand, \u00a732\n  (cantilever with intermediate load).\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  (cantilever, concentrated load, intermediate position).\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 slender beam, midspan load\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 = L / 2.0  # load position from the clamp [m]\nP = 1.0e3  # transverse load [N, downward = -y]\n\nEI = E * I_z\n\n\ndef v_closed_form(x: float) -> float:\n    \"\"\"Closed-form deflection v(x) for the off-tip cantilever load.\"\"\"\n    if x <= a:\n        return P * x**2 * (3.0 * a - x) / (6.0 * EI)\n    return P * a**2 * (3.0 * x - a) / (6.0 * EI)\n\n\ndef theta_closed_form(x: float) -> float:\n    \"\"\"Closed-form slope dv/dx(x) for the off-tip cantilever load.\"\"\"\n    if x <= a:\n        return P * x * (2.0 * a - x) / (2.0 * EI)\n    return P * a**2 / (2.0 * EI)\n\n\nv_at_load_pub = v_closed_form(a)\nv_tip_pub = v_closed_form(L)\ntheta_tip_pub = theta_closed_form(L)\n\nprint(\"Cantilever with off-tip point load\")\nprint(f\"  L = {L} m, load at x = a = {a} m  (a / L = {a / L:.3f})\")\nprint(f\"  E = {E:.2e} Pa, I = {I_z:.3e} m^4, EI = {EI:.3e} N m^2, P = {P} N\")\nprint()\nprint(\"Closed-form references (Roark Table 8 case 1):\")\nprint(f\"  v(a)     = P a^3 / (3 E I)             = {v_at_load_pub * 1e3:+.4e} mm\")\nprint(f\"  v(L)     = P a^2 (3 L - a) / (6 E I)   = {v_tip_pub * 1e3:+.4e} mm\")\nprint(f\"  theta(L) = P a^2 / (2 E I)             = {theta_tip_pub:+.4e} rad\")"
      ]
    },
    {
      "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: full clamp at x = 0\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m.fix(nodes=1, dof=\"ALL\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Apply the load at the closest node to x = a\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "i_load = int(np.argmin(np.abs(xs - a)))\nm.apply_force(int(i_load + 1), fy=-P)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Suppress out-of-plane motion at every node so the strictly-2D\nresponse matches the closed-form Bernoulli kinematics.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for 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": [
        "## 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_at_load_fe = _node_dof(i_load + 1, 1)\nv_tip_fe = _node_dof(N_ELEM + 1, 1)\ntheta_tip_fe = _node_dof(N_ELEM + 1, 5)  # ROTZ \u2014 in-plane rotation\n\n# The applied load points in -y (gravity-like), so FE deflections\n# come out negative.  The closed-form formulas above return the\n# magnitude.  Compare on absolute value across the board so the\n# sign convention drops out.\nerr_load = (abs(v_at_load_fe) - v_at_load_pub) / v_at_load_pub * 100.0\nerr_tip = (abs(v_tip_fe) - v_tip_pub) / v_tip_pub * 100.0\nerr_theta = (abs(theta_tip_fe) - theta_tip_pub) / theta_tip_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(a)':<22}  {v_at_load_fe * 1e3:>10.4f} mm  \"\n    f\"{v_at_load_pub * 1e3:>10.4f} mm  {err_load:>+8.4f}%\"\n)\nprint(f\"{'v(L)':<22}  {v_tip_fe * 1e3:>10.4f} mm  {v_tip_pub * 1e3:>10.4f} mm  {err_tip:>+8.4f}%\")\nprint(\n    f\"{'|theta(L)|':<22}  {abs(theta_tip_fe):>11.6e}  {theta_tip_pub:>11.6e}  {err_theta:>+8.4f}%\"\n)\n\nassert abs(err_load) < 0.1, f\"v(a) deviation {err_load:.4f} % exceeds 0.1 %\"\nassert abs(err_tip) < 0.1, f\"v(L) deviation {err_tip:.4f} % exceeds 0.1 %\"\nassert abs(err_theta) < 0.1, f\"theta(L) deviation {err_theta:.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 = np.array([_node_dof(int(i + 1), 1) for i in range(N_ELEM + 1)])\nanalytical_uy = np.array([v_closed_form(x) for x in xs])\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 1)\",\n)\nax.plot(xs, fe_uy * 1e3, \"o\", color=\"#d62728\", markersize=5, label=\"BEAM2 FE\")\nax.axvline(a, color=\"grey\", ls=\"--\", lw=1.0, label=f\"load at x = a = {a:.3f} m\")\nax.set_xlabel(\"x [m]\")\nax.set_ylabel(\"v(x)  [mm]   (downward deflection negative)\")\nax.set_title(\"Cantilever with off-tip 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 shape\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_tip_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([[a, scale * v_at_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:.2f} m)\",\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) within {abs(err_load):.4f} % of P a\u00b3 / (3 E I).\")\nprint(f\"  \u2022 v(L) within {abs(err_tip):.4f} % of P a\u00b2 (3 L \u2212 a) / (6 E I).\")\nprint(f\"  \u2022 |theta(L)| within {abs(err_theta):.4f} % of P a\u00b2 / (2 E I).\")\nratio = v_tip_pub / (P * L**3 / (3.0 * EI))\nprint(\n    f\"  \u2022 Compared to a tip-loaded cantilever of the same beam, the midspan-loaded \"\n    f\"tip deflects only {ratio:.4f} of the tip-loaded value (= 5/16 for a = L/2).\"\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
}