{
  "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 under uniformly-distributed load\n\nStatically-indeterminate first-order beam: clamped at\n$x = 0$, simply supported at $x = L$, with a\nuniform downward distributed load $q$ over the full\nspan.  The closed-form Euler-Bernoulli solution (Roark &\nYoung 1989, Table 8 case 7; Timoshenko 1955 \u00a740) gives\n\n\\begin{align}R_\\mathrm{prop} = \\tfrac{3}{8} q L,\n    \\qquad\n    R_\\mathrm{clamp} = \\tfrac{5}{8} q L,\n    \\qquad\n    M_\\mathrm{clamp} = \\tfrac{1}{8} q L^{2},\\end{align}\n\nand the maximum (downward) deflection occurs at\n$x^{*} \\approx 0.5785\\, L$ with magnitude\n\n\\begin{align}\\delta_\\mathrm{max}\n        \\;=\\; \\frac{q\\, L^{4}}{185\\, E I},\n    \\qquad\n    \\delta(L/2) = \\frac{q L^{4}}{192\\, E I}.\\end{align}\n\n## Implementation\n\nStandard 3D BEAM2 (Hermite-cubic Bernoulli kinematics, Cook\nTable 16.3-1) line of $N_e = 20$ elements along the\nspan.  Both end-section restraints use the native :meth:`Model.fix`\nAPI:\n\n* node 1 (clamp): all 6 DOFs pinned.\n* node $N_e + 1$ (prop): UY = 0; out-of-plane DOFs\n  (UZ, ROTX, ROTY, ROTZ) suppressed to keep the response purely\n  in-plane.\n\nThe UDL is applied as the consistent nodal-force vector for a\nuniform load on a Hermite beam (Cook \u00a716.5,\n$\\{f_e\\} = q L_e \\{1/2,\\, L_e/12,\\, 1/2,\\, -L_e/12\\}^{T}$),\nwhich is what an APDL ``F,FY,\u2026`` over the line equivalents to.\n\n## References\n\n* Timoshenko, S. P. (1955) *Strength of Materials, Part I*,\n  3rd ed., Van Nostrand, \u00a740.\n* Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for\n  Stress and Strain*, 6th ed., McGraw-Hill, Table 8 case 7.\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, \u00a716.5 (consistent loads on beam elements).\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  # second moment of area  # noqa: E741\nJ = 2.0 * I  # thin-square Saint-Venant J (no torsion in this BC set)\nL = 1.0  # span [m]\nq = 1.0e3  # uniform line load [N/m] (downward \u2192 -y)\n\nN_ELEM = 20\n\n# Closed-form quantities -----------------------------------------------\n\nR_clamp_published = 5.0 / 8.0 * q * L\nR_prop_published = 3.0 / 8.0 * q * L\nM_clamp_published = q * L**2 / 8.0\ndelta_mid_published = q * L**4 / (192.0 * E * I)\ndelta_max_published = q * L**4 / (185.0 * E * I)\nx_star = 0.5785  # location of max deflection / L\n\nprint(f\"q = {q:.1f} N/m, L = {L:.2f} m, EI = {E * I:.3e} N m^2\")\nprint(f\"R_clamp  = 5/8 q L          = {R_clamp_published:.4f} N\")\nprint(f\"R_prop   = 3/8 q L          = {R_prop_published:.4f} N\")\nprint(f\"M_clamp  = 1/8 q L^2        = {M_clamp_published:.4f} N m\")\nprint(f\"\u03b4(L/2)   = q L^4 / (192 EI) = {delta_mid_published:.4e} m\")\nprint(f\"\u03b4_max    = q L^4 / (185 EI) = {delta_max_published:.4e} m  (at x \u2248 {x_star} L)\")"
      ]
    },
    {
      "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\nNode 1 \u2014 clamp (all 6 DOFs pinned).\nNode N_e + 1 \u2014 pin only UY (the prop) and suppress\nout-of-plane DOFs so the response is in the x-y plane.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m.fix(nodes=[1], dof=\"ALL\")\nprop = N_ELEM + 1\nm.fix(nodes=[prop], dof=\"UY\")\n# Suppress out-of-plane / torsional DOFs so the response stays\n# in the x-y plane.  ROTZ is **free** at the prop \u2014 pinning it\n# would turn the simple support into a clamp and the problem\n# would be a clamped-clamped beam (50/50 reaction split) instead\n# of a propped cantilever (5/8 - 3/8 split).\nm.fix(nodes=[prop], dof=\"UZ\")\nm.fix(nodes=[prop], dof=\"ROTX\")\nm.fix(nodes=[prop], dof=\"ROTY\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Consistent UDL on a Hermite beam\n\nFor a single Bernoulli element of length $L_e$ under\nuniform line load $q$, the consistent nodal forces /\nmoments at its two endpoints are\n\n\\begin{align}\\{f_e\\} = q\\,L_e\n    \\begin{Bmatrix} 1/2 \\\\ L_e/12 \\\\ 1/2 \\\\ -L_e/12 \\end{Bmatrix},\\end{align}\n\ni.e. half the segment load goes to each endpoint as a\ntransverse force, and a fixing moment of magnitude\n$q L_e^{2}/12$ accumulates at each end with opposite\nsigns (Cook \u00a716.5).  We assemble these onto the model nodes.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "L_e = L / N_ELEM\nfy_per_node = np.zeros(N_ELEM + 1)\nmz_per_node = np.zeros(N_ELEM + 1)\nfor e in range(N_ELEM):\n    f_half = -q * L_e / 2.0  # half of -q L_e per endpoint (downward)\n    m_left = -q * L_e**2 / 12.0  # = -q L_e^2 / 12 from {-qL/2,-qL^2/12,-qL/2,+qL^2/12}\n    m_right = +q * L_e**2 / 12.0\n    fy_per_node[e] += f_half\n    fy_per_node[e + 1] += f_half\n    mz_per_node[e] += m_left\n    mz_per_node[e + 1] += m_right\n\nfor i in range(N_ELEM + 1):\n    if abs(fy_per_node[i]) > 0:\n        m.apply_force(i + 1, fy=float(fy_per_node[i]))\n    if abs(mz_per_node[i]) > 0:\n        m.apply_force(i + 1, mz=float(mz_per_node[i]))"
      ]
    },
    {
      "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_clamp_y = _react(1, 1)  # UY reaction at clamp\nM_clamp_z = _react(1, 5)  # ROTZ reaction at clamp\nR_prop_y = _react(prop, 1)  # UY reaction at prop\n\nprint()\nprint(\"reactions  (femorph-solver) \u2192 (analytical)\")\nprint(f\"  R_clamp_UY  = {R_clamp_y:+.4f}  \u2192  {+R_clamp_published:+.4f} N\")\nprint(f\"  R_prop_UY   = {R_prop_y:+.4f}  \u2192  {+R_prop_published:+.4f} N\")\nprint(f\"  M_clamp_RZ  = {M_clamp_z:+.4f}  \u2192  {+M_clamp_published:+.4f} N m\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Equilibrium check + tolerance assertions\n\nSum of vertical reactions must balance the applied UDL\n(5/8 + 3/8 = 1) and the moment about the clamp from the\ndownward load + prop reaction must be zero.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "assert np.isclose(R_clamp_y + R_prop_y, q * L, rtol=1e-6), \"vertical equilibrium failed\"\nassert np.isclose(R_clamp_y, R_clamp_published, rtol=1e-6), \"clamp UY reaction off\"\nassert np.isclose(R_prop_y, R_prop_published, rtol=1e-6), \"prop UY reaction off\"\nassert np.isclose(M_clamp_z, M_clamp_published, rtol=1e-6), \"clamp moment off\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Mid-span deflection\n\nHermite cubic bases recover Euler-Bernoulli kinematics\nexactly for prismatic beams under nodal loads, so the\nmid-span UY equals the analytical $\\delta(L/2) = q L^4\n/ (192 E I)$ to machine precision on a uniform mesh.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mid_node = N_ELEM // 2 + 1\nmid_uy_val = float(res.displacement[np.where((dof[:, 0] == mid_node) & (dof[:, 1] == 1))[0][0]])\nerr_mid = (mid_uy_val - (-delta_mid_published)) / (-delta_mid_published)\nprint()\nprint(f\"\u03b4(L/2) computed   = {mid_uy_val:+.4e} m\")\nprint(f\"\u03b4(L/2) published  = {-delta_mid_published:+.4e} m\")\nprint(f\"relative error    = {err_mid * 100:+.4f} %\")\nassert abs(err_mid) < 1.0e-4, \"mid-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.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
}