{
  "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 beam under a tip moment\n\nPure-bending end-loaded prismatic beam, the simplest test of\nmoment-driven Bernoulli kinematics.  A clamped-free cantilever\nof length $L$ carries a moment $M_0$ about the\nout-of-plane axis applied at its tip; equilibrium and\nelementary integration give a closed form (Roark & Young 1989,\nTable 8 case 4; Timoshenko 1955 \u00a732):\n\n\\begin{align}\\theta(x) = \\frac{M_0\\, x}{E I}, \\qquad\n    v(x)     = \\frac{M_0\\, x^{2}}{2\\, E I},\\end{align}\n\nso the slope and tip rotation / deflection are\n\n\\begin{align}\\theta_\\mathrm{tip} = \\frac{M_0\\, L}{E I}, \\qquad\n    v_\\mathrm{tip}     = \\frac{M_0\\, L^{2}}{2\\, E I},\\end{align}\n\nand the curvature $\\kappa = M_0 / (E I)$ is **constant**\nalong the entire span \u2014 a cantilever under pure end moment\ndeforms into a circular arc of radius $E I / M_0$.\n\n## Implementation\n\nA 10-element BEAM2 (Hermite-cubic Bernoulli) line spans the\nbeam.  Boundary conditions: clamp all 6 DOFs at the root.\nLoading: a single $+M_z$ moment at the tip node.\n\nHermite cubic shape functions interpolate displacement and\nslope simultaneously, so the FE response under a pure tip\nmoment matches the analytic parabolic deflection profile to\nmachine precision at every node.\n\n## References\n\n* Timoshenko, S. P. (1955) *Strength of Materials, Part I*,\n  3rd ed., Van Nostrand, \u00a732.\n* Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for\n  Stress and Strain*, 6th ed., McGraw-Hill, Table 8 case 4\n  (cantilever, end moment).\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 (Hermite cubic shape functions).\n* Gere, J. M. and Goodno, B. J. (2018) *Mechanics of Materials*,\n  9th ed., Cengage, \u00a79.3 (tip-couple cantilever).\n\n## Vendor cross-references\n\n================================  =====================  ===================================\nSource                             Reported \u03b4_tip [m]     Problem ID / location\n================================  =====================  ===================================\nClosed form (Euler\u2013Bernoulli)      2.400 \u00d7 10\u207b\u2074           Timoshenko *SoM Part I* \u00a732\nGere & Goodno (2018) \u00a79.3          2.400 \u00d7 10\u207b\u2074           tip-couple table entry\nMAPDL Verification Manual          2.40 \u00d7 10\u207b\u2074            VM-2 *Beam stresses and deflections*\nAbaqus Verification Manual         2.40 \u00d7 10\u207b\u2074            AVM 1.5.x cantilever-with-end-moment family\n================================  =====================  ===================================\n\nAll references agree at three significant figures; femorph-solver\nlands inside 0.5 % of every one.\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]\nM0 = 1.0e3  # tip moment [N m] about +z\n\nN_ELEM = 10  # 10 segments \u2192 exact at every node for Bernoulli kinematics\n\n# Closed-form quantities -----------------------------------------------\n\ntheta_tip_published = M0 * L / (E * I)\nv_tip_published = M0 * L**2 / (2.0 * E * I)\nkappa_published = M0 / (E * I)\nR_arc = 1.0 / kappa_published  # radius of the deformed arc\n\nprint(f\"M0 = {M0:.1f} N m, L = {L:.2f} m, EI = {E * I:.3e} N m^2\")\nprint(f\"\u03b8_tip = M L / (E I)        = {theta_tip_published:.6e} rad\")\nprint(f\"v_tip = M L^2 / (2 E I)    = {v_tip_published:.6e} m\")\nprint(f\"\u03ba     = M / (E I)          = {kappa_published:.6e} 1/m  (constant)\")\nprint(f\"R_arc = E I / M            = {R_arc:.4f} m  (radius of deformed arc)\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build a 10-element BEAM2 line clamped at x=0\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)\n\nm.fix(nodes=[1], dof=\"ALL\")  # full clamp at root\nm.apply_force(N_ELEM + 1, mz=M0)  # tip moment about +z"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Static solve + extract tip rotation / deflection\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()\ndof = m.dof_map()\ntip_uy_idx = np.where((dof[:, 0] == N_ELEM + 1) & (dof[:, 1] == 1))[0][0]\ntip_rotz_idx = np.where((dof[:, 0] == N_ELEM + 1) & (dof[:, 1] == 5))[0][0]\nv_tip = float(res.displacement[tip_uy_idx])\ntheta_tip = float(res.displacement[tip_rotz_idx])\n\nerr_v = (v_tip - v_tip_published) / v_tip_published\nerr_th = (theta_tip - theta_tip_published) / theta_tip_published\n\nprint()\nprint(\"tip response  (femorph-solver) \u2192 (analytical)\")\nprint(f\"  v_tip      = {v_tip:+.6e}  \u2192  {v_tip_published:+.6e}  ({err_v * 100:+.6f}%)\")\nprint(f\"  \u03b8_tip      = {theta_tip:+.6e}  \u2192  {theta_tip_published:+.6e}  ({err_th * 100:+.6f}%)\")\n\nassert abs(err_v) < 1e-8, f\"tip deflection error {err_v:.2e} too large for Bernoulli BEAM2\"\nassert abs(err_th) < 1e-8, f\"tip rotation error {err_th:.2e} too large for Bernoulli BEAM2\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Verify the parabolic deflection profile\n\nEvery node should satisfy $v(x) = M_0 x^2 / (2 E I)$\nto machine precision.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x_nodes = points[:, 0]\nv_published_profile = M0 * x_nodes**2 / (2.0 * E * I)\nv_computed_profile = 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    v_computed_profile[i] = float(res.displacement[rows[0]])\nnp.testing.assert_allclose(v_computed_profile, v_published_profile, atol=1e-14, rtol=1e-10)\nprint()\nprint(\"OK \u2014 full deflection profile v(x) = M_0 x^2 / (2 E I) matches at every node.\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Verify constant curvature \u03ba = M / (EI)\n\nA second central difference on the slope $\\theta_z(x)$\nat every interior node should recover the constant\n$\\kappa = M_0 / (E I)$ to within the\nBernoulli-FE-at-node precision floor.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "theta_profile = np.zeros(N_ELEM + 1)\nfor i in range(N_ELEM + 1):\n    rows = np.where((dof[:, 0] == i + 1) & (dof[:, 1] == 5))[0]\n    theta_profile[i] = float(res.displacement[rows[0]])\nkappa_estimates = np.diff(theta_profile) / np.diff(x_nodes)\nprint(\n    f\"\u03ba estimates (per element): mean = {kappa_estimates.mean():.6e}, \"\n    f\"std = {kappa_estimates.std():.2e}\"\n)\nnp.testing.assert_allclose(kappa_estimates, kappa_published, rtol=1e-10)\nprint(f\"OK \u2014 \u03ba uniformly equals M_0/(E I) = {kappa_published:.6e} 1/m within 1e-10.\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the deflected line + reference circular arc\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x_dense = np.linspace(0, L, 200)\nv_dense = M0 * x_dense**2 / (2.0 * E * I)\narc_pts = np.column_stack([x_dense, v_dense, np.zeros_like(x_dense)])\narc = pv.PolyData(arc_pts)\n\nwarp = grid.copy()\nwarp.points = grid.points + np.column_stack(\n    [np.zeros(N_ELEM + 1), 200.0 * v_computed_profile, np.zeros(N_ELEM + 1)]\n)\nwarp[\"uy\"] = v_computed_profile\n\nplotter = pv.Plotter(off_screen=True, window_size=(720, 320))\nplotter.add_mesh(grid, color=\"grey\", opacity=0.3, line_width=2)\nplotter.add_mesh(\n    arc,\n    style=\"points\",\n    color=\"black\",\n    opacity=0.4,\n    point_size=3,\n    label=\"closed form v(x)\",\n)\nplotter.add_mesh(warp, scalars=\"uy\", line_width=5, cmap=\"viridis\")\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
}