{
  "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# Reaction forces \u2014 global equilibrium and root moment from a tip load\n\nEvery linear static solve produces two displacement fields: the\n``(N,)`` solved-for ``displacement`` at the free DOFs, and the\n``(N,)`` ``reaction`` field nonzero only at the constrained DOFs \u2014\nthe force the supports must carry to hold the prescribed boundary\ndisplacements.  Sum the reactions and you reconstruct the support\nloads engineers actually need: bolt loads, weld loads, foundation\nreactions.\n\nTwo textbook checks every static solve must pass (Hibbeler 2017 \u00a79.5,\nCook et al. 2002 \u00a711.6):\n\n\\begin{align}\\mathbf{F}_\\text{ext, free} + \\mathbf{F}_\\text{reaction} = \\mathbf{0},\n    \\qquad\n    \\sum (\\mathbf{r} \\times \\mathbf{F})_\\text{free}\n      + \\sum (\\mathbf{r} \\times \\mathbf{F})_\\text{reaction}\n      = \\mathbf{0}.\\end{align}\n\nA tip-loaded cantilever is the cleanest demo because every term has a\nsimple closed form:\n\n* Total clamped-face reaction in $y$ equals $-P$\n  (shear balance).\n* Net $x$-reaction is zero, but the *individual* node-by-node\n  $R_x$ values are nonzero \u2014 they form an internal couple whose\n  arm \u00d7 force product gives the root bending moment $M_z = -PL$.\n\nThe second point is what surprises new users: the clamp face is in\n*compression on the bottom, tension on the top*, and that couple is\nexactly what carries the bending moment back into the support.\n\n## Implementation\n\nA slender HEX8-EAS cantilever (L = 2 m, 0.05 \u00d7 0.05 m, 40 \u00d7 3 \u00d7 3\nmesh) clamped at $x = 0$ and loaded with $P = -1000$ N\nin $y$ distributed uniformly across the tip face nodes.\n:attr:`StaticResult.reaction` is the DOF-indexed reaction-force\nvector \u2014 index it via the model's DOF map to read the components per\nnode.\n\n## References\n\n* Hibbeler, R. C. (2017) *Mechanics of Materials*, 10th ed., Pearson,\n  \u00a79.5 \u2014 internal stress resultants from cross-section equilibrium.\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002)\n  *Concepts and Applications of Finite Element Analysis*, 4th ed.,\n  Wiley, \u00a711.6 \u2014 reaction-force recovery in displacement-based FEM.\n* Bathe, K.-J. (2014) *Finite Element Procedures*, 2nd ed.,\n  Prentice Hall, \u00a74.2.2 \u2014 partitioned static system and the\n  recovered reactions.\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": [
        "## Build the cantilever\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "E = 2.0e11\nNU = 0.30\nRHO = 7850.0\nL = 2.0\nWIDTH = 0.05\nHEIGHT = 0.05\nP_TOTAL = -1000.0  # N \u2014 total tip load in -y\n\nNX, NY, NZ = 40, 3, 3\nxs = np.linspace(0.0, L, NX + 1)\nys = np.linspace(0.0, WIDTH, NY + 1)\nzs = np.linspace(0.0, HEIGHT, NZ + 1)\ngrid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing=\"ij\")).cast_to_unstructured_grid()\n\nm = femorph_solver.Model.from_grid(grid)\nm.assign(\n    ELEMENTS.HEX8(integration=\"enhanced_strain\"),\n    material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n)\n\npts = np.asarray(m.grid.points)\nnode_nums = np.asarray(m.grid.point_data[\"ansys_node_num\"], dtype=np.int64)\nclamped_pt_idx = np.where(pts[:, 0] < 1e-9)[0]\ntip_pt_idx = np.where(pts[:, 0] > L - 1e-9)[0]\n\nm.fix(nodes=node_nums[clamped_pt_idx].tolist(), dof=\"ALL\")\n\n# Distribute the total tip load uniformly across the tip-face nodes.\nload_per_node = P_TOTAL / float(tip_pt_idx.size)\nfor k in tip_pt_idx:\n    m.apply_force(int(node_nums[k]), fy=load_per_node)\n\nprint(\"Tip-loaded cantilever \u2014 reaction-force recovery\")\nprint(f\"  L = {L} m, cross = {WIDTH} \u00d7 {HEIGHT} m, {NX}\u00d7{NY}\u00d7{NZ} HEX8-EAS\")\nprint(f\"  P_tip = {P_TOTAL:.1f} N in -y, distributed across {tip_pt_idx.size} face nodes\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Static solve\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()\nprint(f\"\\n  {res!r}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Sum reactions at the clamped face\n\n``StaticResult.reaction`` is ``(n_dof,)`` and nonzero only at the\nconstrained DOFs.  ``Model.dof_map()`` gives ``(node_number, dof_idx)``\nper row so we can fold the flat vector into per-node force vectors.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dof_map = m.dof_map()\nreaction = np.asarray(res.reaction, dtype=np.float64)\nclamped_node_set = set(int(n) for n in node_nums[clamped_pt_idx])\n\nn_points = m.grid.n_points\nnode_to_pt_idx = {int(node_nums[i]): i for i in range(n_points)}\n\n# Gather (R_x, R_y, R_z) per clamped node.\nclamp_nodes = sorted(clamped_node_set)\nclamp_reactions = np.zeros((len(clamp_nodes), 3), dtype=np.float64)\nfor row, (node, dof_idx) in enumerate(dof_map.tolist()):\n    if int(node) in clamped_node_set and dof_idx < 3:\n        i = clamp_nodes.index(int(node))\n        clamp_reactions[i, int(dof_idx)] += reaction[row]\n\nr_total = clamp_reactions.sum(axis=0)\nprint()\nprint(f\"  \u03a3R_x  = {r_total[0]:+.3f} N   (expected  0)\")\nprint(f\"  \u03a3R_y  = {r_total[1]:+.3f} N   (expected {-P_TOTAL:+.3f})\")\nprint(f\"  \u03a3R_z  = {r_total[2]:+.3f} N   (expected  0)\")\n\nresidual = np.abs(r_total + np.array([0.0, P_TOTAL, 0.0]))\nprint(f\"\\n  global force residual  \u2016\u03a3R + \u03a3F_applied\u2016 = {np.linalg.norm(residual):.3e} N\")\nassert np.linalg.norm(residual) < 1e-6 * abs(P_TOTAL), \"global force balance violated\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Bending-moment recovery from the root couple\n\nEven though \u03a3R_x = 0, the per-node R_x values are *not* zero \u2014 they\nform a couple about the neutral axis:\n\n    M_z(root) = \u03a3 over clamp nodes of  y \u00b7 R_x  +  (something)\n\nFor a transverse load in y on a tip-loaded cantilever the dominant\ncouple is the ``y * R_x`` integral.  Equivalently, the root moment\nis $-P \\cdot L$ (textbook), so we expect\n$M_z \\approx 2000\\,\\text{N\u00b7m}$ for our configuration.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "clamp_pts = pts[[node_to_pt_idx[n] for n in clamp_nodes]]\ny = clamp_pts[:, 1] - 0.5 * WIDTH  # offset to neutral axis\n\n# \u03a3 (r \u00d7 R)_z = \u03a3 (x \u00b7 R_y \u2212 y \u00b7 R_x).  At the clamp, x = 0, so\n# M_z(at clamp) = \u2212 \u03a3 y \u00b7 R_x.  Add the contribution from R_y across\n# x = 0 (which is zero by construction).\nm_z_clamp = -float(np.sum(y * clamp_reactions[:, 0]))\nm_z_pub = -P_TOTAL * L\nprint()\nprint(f\"  M_z(at clamp) recovered from \u03a3yR_x  = {m_z_clamp:+.3f} N\u00b7m\")\nprint(f\"  M_z(at clamp) closed form  \u2212P\u00b7L     = {m_z_pub:+.3f} N\u00b7m\")\nprint(\n    f\"  relative error                       = {abs(m_z_clamp - m_z_pub) / abs(m_z_pub) * 100:.3f}%\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Per-node table for the four extreme corners\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print()\nprint(f\"  {'node':>6}  {'y':>8}  {'z':>8}  {'R_x [N]':>12}  {'R_y [N]':>12}  {'R_z [N]':>12}\")\nprint(\"  \" + \"-\" * 64)\n# Sort by descending |R_x| to highlight the bending couple.\norder = np.argsort(-np.abs(clamp_reactions[:, 0]))\nfor idx in order[:8]:\n    node = clamp_nodes[idx]\n    p = clamp_pts[idx]\n    r = clamp_reactions[idx]\n    print(\n        f\"  {node:>6}  {p[1]:>8.4f}  {p[2]:>8.4f}  {r[0]:>+12.3f}  {r[1]:>+12.3f}  {r[2]:>+12.3f}\"\n    )"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the reaction-force vectors at the clamped face\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reaction_grid = pv.PolyData(clamp_pts)\nreaction_grid[\"reaction\"] = clamp_reactions\nreaction_grid[\"|R|\"] = np.linalg.norm(clamp_reactions, axis=1)\n\n# Scale glyphs so the longest reaction arrow reaches 25 % of beam length.\nmax_reaction = float(reaction_grid[\"|R|\"].max())\narrow_scale = 0.25 * L / max_reaction if max_reaction > 0 else 1.0\n\nreaction_grid.set_active_vectors(\"reaction\")\narrows = reaction_grid.glyph(orient=\"reaction\", scale=\"|R|\", factor=arrow_scale)\n\nplotter = pv.Plotter(off_screen=True, window_size=(900, 540))\nplotter.add_mesh(m.grid, style=\"wireframe\", color=\"grey\", opacity=0.4, line_width=1)\nplotter.add_mesh(\n    arrows,\n    scalars=\"|R|\",\n    cmap=\"plasma\",\n    show_scalar_bar=True,\n    scalar_bar_args={\"title\": \"|R|  [N]\"},\n)\nplotter.add_text(\n    f\"clamp-face reactions  \u2014  \u03a3R_y = {r_total[1]:.1f} N\", position=\"upper_left\", font_size=11\n)\nplotter.view_xy()\nplotter.camera.zoom(1.1)\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(\n    \"  \u2022 StaticResult.reaction is a DOF-indexed (N,) array \u2014 nonzero only at \"\n    \"constrained DOFs.  Index it via Model.dof_map() to get per-node force \"\n    \"vectors at the supports.\"\n)\nprint(\n    \"  \u2022 Global force balance: \u03a3F_applied + \u03a3F_reaction = 0 to machine precision.  \"\n    \"Use this as a sanity check on every solve before trusting downstream \"\n    \"stress / strain recovery.\"\n)\nprint(\n    \"  \u2022 Internal moments are carried by NODE-TO-NODE couples at the \"\n    \"support: \u03a3R_x = 0 net, but R_x varies linearly with the lever-arm \"\n    \"coordinate.  \u03a3 y\u00b7R_x reconstructs the root bending moment to within \"\n    \"FEM discretisation error.\"\n)\nprint(\n    \"  \u2022 Reaction-force visualisation at the clamp is the engineer's first \"\n    \"check that loads are flowing through the structure as intended \u2014 \"\n    \"asymmetries here usually point at a misapplied BC or a missing fix.\"\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
}