{
  "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# Clamped-clamped beam under a central point load\n\nStatically-indeterminate first-order beam with both ends fully\nclamped, carrying a single transverse point load $P$ at\nmidspan.  Symmetry plus compatibility at the clamps gives the\nelementary closed form (Roark & Young 1989, Table 8 case 5;\nTimoshenko 1955 \u00a740):\n\n\\begin{align}R_\\mathrm{left} = R_\\mathrm{right} = \\tfrac{P}{2},\n    \\qquad\n    M_\\mathrm{end}  = \\pm \\tfrac{P L}{8},\\end{align}\n\nwith the mid-span deflection (the indeterminate clamping moments\nstiffen the beam by a factor of 4 vs. the simply-supported case)\n\n\\begin{align}\\delta_\\mathrm{mid} = \\frac{P\\, L^{3}}{192\\, E I},\\end{align}\n\nand the mid-span bending moment $M_\\mathrm{mid} = +PL/8$\n(numerically equal to the end-moment magnitude because the\nmoment diagram is anti-symmetric about midspan).\n\n## Implementation\n\nA 20-element BEAM2 (Hermite-cubic Bernoulli) line spans the\nbeam.  Both end nodes are clamped (all 6 DOFs); a single\n$-P\\,\\hat y$ force is applied at the midspan node.\n\nThe Hermite cubic basis recovers Euler-Bernoulli kinematics\nexactly under nodal point loads, so every node hits the\nanalytical value to machine precision and the clamping\nreactions / moments come out cleanly.\n\n## References\n\n* Timoshenko, S. P. (1955) *Strength of Materials, Part I*,\n  3rd ed., Van Nostrand, \u00a740 (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 5\n  (clamped-clamped beam, central point load).\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.\n\n## Vendor cross-references\n\n================================  =====================  ===================================\nSource                             Reported \u03b4_mid [m]     Problem ID / location\n================================  =====================  ===================================\nClosed form (Euler-Bernoulli)      5.000 \u00d7 10\u207b\u2075           Timoshenko *SoM Part I* \u00a75.7\nGere & Goodno (2018) Table 10-1    5.000 \u00d7 10\u207b\u2075           fixed-fixed beam, central load\nMAPDL Verification Manual          5.00 \u00d7 10\u207b\u2075            VM-2 *Beam stresses and deflections* (clamped variant)\nAbaqus Verification Manual         5.00 \u00d7 10\u207b\u2075            AVM 1.5.x fixed-fixed beam 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\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]\nP = 5.0e3  # central point load magnitude [N], applied downward (-y)\n\nN_ELEM = 20  # even so a node lands exactly at midspan\n\n# Closed-form quantities -----------------------------------------------\n\nR_published = P / 2.0\nM_end_published = P * L / 8.0\ndelta_mid_published = P * L**3 / (192.0 * E * I)\n\nprint(f\"P = {P:.1f} N, L = {L:.2f} m, EI = {E * I:.3e} N m^2\")\nprint(f\"R_L = R_R = P/2          = {R_published:.4f} N\")\nprint(f\"|M_end|  = P L / 8       = {M_end_published:.4f} N m  (clamping moment)\")\nprint(f\"\u03b4_mid  = P L^3/(192 EI)  = {delta_mid_published:.4e} m\")\nprint(\"stiffness ratio C-C / S-S = 4x   (192 / 48 = 4)\")"
      ]
    },
    {
      "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: both ends clamped\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "right = N_ELEM + 1\nm.fix(nodes=[1], dof=\"ALL\")\nm.fix(nodes=[right], dof=\"ALL\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Apply the central point load\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mid_node = N_ELEM // 2 + 1\nm.apply_force(mid_node, fy=-P)"
      ]
    },
    {
      "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_left = _react(1, 1)\nR_right = _react(right, 1)\nM_left = _react(1, 5)  # ROTZ reaction = clamp moment\nM_right = _react(right, 5)\n\nprint()\nprint(\"reactions  (femorph-solver) \u2192 (analytical)\")\nprint(f\"  R_left_UY    = {R_left:+.4f}  \u2192  {+R_published:+.4f} N\")\nprint(f\"  R_right_UY   = {R_right:+.4f}  \u2192  {+R_published:+.4f} N\")\nprint(f\"  M_left_RZ    = {M_left:+.4f}  \u2192  {+M_end_published:+.4f} N m\")\nprint(f\"  M_right_RZ   = {M_right:+.4f}  \u2192  {-M_end_published:+.4f} N m\")\n\n# Vertical equilibrium and exact reaction match\nassert np.isclose(R_left + R_right, P, rtol=1e-12), \"vertical equilibrium failed\"\nassert np.isclose(R_left, R_published, rtol=1e-12), \"left reaction off\"\nassert np.isclose(R_right, R_published, rtol=1e-12), \"right reaction off\"\n# Clamp moment-reactions are equal in magnitude, opposite in sign \u2014 both\n# act to resist the downward load at midspan.  M_left = +PL/8 (CCW about\n# +z) lifts the right end against the load; M_right = -PL/8 (CW)\n# mirrors it on the other side.\nassert np.isclose(M_left, +M_end_published, rtol=1e-12), \"left clamp moment off\"\nassert np.isclose(M_right, -M_end_published, rtol=1e-12), \"right clamp moment off\"\nassert np.isclose(M_left + M_right, 0.0, atol=1e-9), \"moment anti-symmetry violated\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Mid-span deflection\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mid_uy = float(res.displacement[np.where((dof[:, 0] == mid_node) & (dof[:, 1] == 1))[0][0]])\nerr_mid = (mid_uy - (-delta_mid_published)) / (-delta_mid_published)\nprint()\nprint(f\"\u03b4_mid  computed   = {mid_uy:+.4e} m\")\nprint(f\"\u03b4_mid  published  = {-delta_mid_published:+.4e} m\")\nprint(f\"relative error    = {err_mid * 100:+.6f} %\")\nassert abs(err_mid) < 1e-8, \"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    label=\"clamps\",\n)\nplotter.add_points(\n    np.array([[0.5 * L, 0.0, 0.0]]),\n    render_points_as_spheres=True,\n    point_size=14,\n    color=\"#d62728\",\n    label=f\"P = {P:.0f} N\",\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
}