{
  "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# Tutorial 1 \u2014 Cantilever beam under combined loading\n\nYou're tasked with sizing a steel cantilever bracket for an\noverhead-mounted equipment console.  The bracket carries the\nconsole's dead weight (a uniform load along its span), an\nattached cable that pulls down at the tip, and a small\napplied moment from an off-axis fastener at the tip.  A\npreliminary cross-section is on the table; you need to confirm\nthat\n\n1. the maximum bending stress stays below the material's design\n   allowable, and\n2. the answer doesn't drift as the mesh is refined \u2014 i.e., the\n   solution has converged.\n\nThis is a textbook design-by-analysis exercise.  This tutorial\nwalks through it end-to-end with femorph-solver:\n\n* **Step 1** \u2014 define the geometry and material.\n* **Step 2** \u2014 apply the boundary conditions.\n* **Step 3** \u2014 apply the three-component combined load.\n* **Step 4** \u2014 run a static solve.\n* **Step 5** \u2014 recover and inspect the stress field.\n* **Step 6** \u2014 check the result against an engineering allowable.\n* **Step 7** \u2014 refine the mesh and confirm convergence.\n\nBy the end you should be comfortable with the entire static-\nanalysis loop on a 3D solid model.\n\n## Theoretical reference\n\nThe Euler-Bernoulli closed form for a slender clamped cantilever\nunder a combined tip force $P$, tip moment $M_0$,\nand uniformly distributed transverse load $q$ is the\nlinear superposition of three textbook cases (Timoshenko 1955\n\u00a75.4, Cook et al. 2002 \u00a72.4):\n\n\\begin{align}:label: combined-tip-deflection\n\n    \\delta_\\text{tip}\n    \\;=\\;\n    \\underbrace{\\frac{P L^{3}}{3 E I}}_\\text{tip-shear}\n    \\;+\\;\n    \\underbrace{\\frac{M_0 L^{2}}{2 E I}}_\\text{tip-moment}\n    \\;+\\;\n    \\underbrace{\\frac{q L^{4}}{8 E I}}_\\text{UDL}.\\end{align}\n\nThe corresponding **maximum bending moment at the root** is the\nsum of contributions:\n\n\\begin{align}M_\\text{root} = P\\, L \\;+\\; M_0 \\;+\\; \\tfrac{1}{2} q\\, L^{2},\\end{align}\n\nand the extreme-fibre bending stress $\\sigma_\\text{max}\n= M_\\text{root}\\, c / I$ (with $c = h/2$) gives the\nquantity we'll compare against the design allowable.\n\nA 3-D HEX8 EAS solid mesh recovers these closed-form values\nexactly to within mesh-discretisation error \u2014 and shows you the\n*spatial distribution* of stress that a 1-D Bernoulli line\nelement cannot.\n\nCompanion theory in :doc:`/reference/theory/variational_form`\nand :doc:`/reference/elements/hex8`.\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\nfrom femorph_solver.recover import compute_nodal_stress, stress_invariants"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 1 \u2014 geometry and material\n\nWe model a slender prismatic steel cantilever, $L = 1\\,\\text{m}$\nlong with a square $50 \\times 50\\,\\text{mm}$ cross-section.\nThe slenderness ratio $L/h = 20$ is well into the Euler-\nBernoulli regime, so the closed forms above are an excellent\naccuracy reference.\n\nWe pick HEX8 with **enhanced-assumed-strain** (Simo & Rifai 1990)\nbecause plain bilinear hexes shear-lock badly on slender bending.\nThe EAS variant adds nine static-condensed internal modes that\nfix the locking \u2014 the natural choice for any 3-D solid bending\nproblem.  See :doc:`/reference/elements/hex8` for the full\ntreatment.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "E = 2.0e11  # Young's modulus, Pa (steel)\nNU = 0.30  # Poisson's ratio\nRHO = 7850.0  # density, kg/m\u00b3 (steel)\nL = 1.0  # length, m\nWIDTH = 0.05  # cross-section width, m\nHEIGHT = 0.05  # cross-section height, m\n\nI_section = WIDTH * HEIGHT**3 / 12.0  # second moment of area\nc_extreme = HEIGHT / 2.0  # extreme-fibre distance\nsigma_allowable = 165.0e6  # design allowable, 165 MPa (steel ~SF 1.5 vs \u03c3_y 250 MPa)\n\nprint(\"Cantilever bracket \u2014 combined-loading design check\")\nprint(f\"  L = {L} m, cross = {WIDTH} \u00d7 {HEIGHT} m  (L/h = {L / HEIGHT:.0f})\")\nprint(f\"  E = {E / 1e9:.0f} GPa, \u03bd = {NU}, \u03c1 = {RHO} kg/m\u00b3\")\nprint(f\"  I = {I_section:.3e} m\u2074\")\nprint(f\"  \u03c3_allowable = {sigma_allowable / 1e6:.0f} MPa  (design target)\")\n\n\ndef build_cantilever(nx: int, ny: int, nz: int) -> femorph_solver.Model:\n    \"\"\"Build a HEX8-EAS cantilever mesh of size ``nx \u00d7 ny \u00d7 nz``.\"\"\"\n    xs = np.linspace(0.0, L, nx + 1)\n    ys = np.linspace(0.0, WIDTH, ny + 1)\n    zs = np.linspace(0.0, HEIGHT, nz + 1)\n    grid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing=\"ij\")).cast_to_unstructured_grid()\n\n    m = femorph_solver.Model.from_grid(grid)\n    m.assign(\n        ELEMENTS.HEX8(integration=\"enhanced_strain\"),\n        material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n    )\n    return m"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 2 \u2014 boundary conditions\n\nA cantilever is *fully clamped* at the root.  In a 3-D solid\nmodel that means every node on the $x = 0$ face has\nall three translations pinned.  This is what\n:meth:`Model.fix(nodes=..., dof=\"ALL\") <femorph_solver.Model.fix>`\nexpresses \u2014 the ``\"ALL\"`` keyword is the natural shorthand for\n\"lock every translational DOF on these nodes\".\n\nIdentifying the right node set is a recurring task in\npre-processing.  We read coordinates off ``Model.grid.points``\nand use NumPy-style boolean masks to pick the face we want.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m = build_cantilever(nx=40, ny=3, nz=3)\npts = np.asarray(m.grid.points)\nnode_nums = np.asarray(m.grid.point_data[\"ansys_node_num\"], dtype=np.int64)\n\nclamped = np.where(pts[:, 0] < 1e-9)[0]\nm.fix(nodes=node_nums[clamped].tolist(), dof=\"ALL\")\nprint(f\"\\n  Clamped {len(clamped)} root-face nodes (full-fix).\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 3 \u2014 combined loading\n\nThree load components, applied to the tip face at $x = L$:\n\n* a downward dead-weight equivalent shear $P = -300\\,\\text{N}$\n  (the cable + console reaction force at the bracket tip),\n* a moment $M_0 = -10\\,\\text{N\u00b7m}$ about the $z$ axis\n  from the off-axis fastener,\n* a uniformly distributed line load $q = -200\\,\\text{N/m}$\n  along the span \u2014 the bracket's self-weight per unit length plus\n  the console's distributed reaction.\n\nDistributing each load across the appropriate node set is the\n3-D-solid analogue of \"applying a moment at a point\" in a 1-D\nbeam model.  For the tip force we lump $P$ evenly across\nthe tip face; for the tip moment we apply equal-and-opposite\nx-forces on the top and bottom edges of the tip face (the\n**mechanical couple** form of an applied moment); for the UDL\nwe distribute $q$ along the span using the trapezoid-\nrule convention from\n`sphx_glr_gallery_verification_example_verify_cantilever_udl.py`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "P_tip = -300.0  # N, downward in -z (tip shear, bending about y)\nM_y_tip = -10.0  # N\u00b7m, about y axis (bends in xz plane, same as P)\nq_udl = -200.0  # N/m, downward distributed (UDL, also -z, also bends about y)\n\n# All three loads now bend the cantilever in the **xz plane**, so the\n# Euler-Bernoulli superposition formula applies directly.\n\n# --- tip face: shear + moment -----------------------------------------\ntip_face = np.where(pts[:, 0] > L - 1e-9)[0]\n# For a moment about y, the couple acts on the +z and -z extreme\n# fibres \u2014 tension on top (+z = HEIGHT), compression on bottom (z = 0).\ntip_top_edge = tip_face[pts[tip_face, 2] > HEIGHT - 1e-9]\ntip_bot_edge = tip_face[pts[tip_face, 2] < 1e-9]\n\n# Shear: distribute P_tip uniformly over the tip-face nodes (in -z).\nfz_per_node = P_tip / float(tip_face.size)\nfor n in tip_face:\n    m.apply_force(int(node_nums[n]), fz=fz_per_node, accumulate=True)\n\n# Moment about y applied as an x-direction tension/compression couple\n# at the extreme fibres: total resultant force = M_y / HEIGHT at each\n# edge, signs opposite.  Net z-force from this couple is zero so it\n# doesn't contaminate the shear balance check.\ncouple_force = M_y_tip / HEIGHT\nfor n in tip_top_edge:\n    m.apply_force(int(node_nums[n]), fx=+couple_force / float(tip_top_edge.size))\nfor n in tip_bot_edge:\n    m.apply_force(int(node_nums[n]), fx=-couple_force / float(tip_bot_edge.size))\n\n# --- span: distributed UDL on the top face ----------------------------\ntop_face = np.where(pts[:, 2] > HEIGHT - 1e-9)[0]\nnx_steps = 40\ndx = L / nx_steps\nny_steps = 3\nfor n in top_face:\n    x = pts[n, 0]\n    col_w = 0.5 if (x < 1e-9 or x > L - 1e-9) else 1.0\n    y = pts[n, 1]\n    y_w = 0.5 if (y < 1e-9 or y > WIDTH - 1e-9) else 1.0\n    fz = (q_udl * dx * col_w * y_w) / ny_steps\n    m.apply_force(int(node_nums[n]), fz=fz, accumulate=True)\n\n# Closed-form predictions (Euler-Bernoulli):\ndelta_pred = (\n    abs(P_tip) * L**3 / (3.0 * E * I_section)\n    + abs(M_y_tip) * L**2 / (2.0 * E * I_section)\n    + abs(q_udl) * L**4 / (8.0 * E * I_section)\n)\nM_root_pred = abs(P_tip) * L + abs(M_y_tip) + 0.5 * abs(q_udl) * L**2\nsigma_max_pred = M_root_pred * c_extreme / I_section\n\nprint()\nprint(\"  Closed-form Euler-Bernoulli predictions:\")\nprint(f\"    \u03b4_tip    (super-position)        = {delta_pred * 1e3:.4f} mm\")\nprint(f\"    M_root   (P L + M\u2080 + q L\u00b2/2)     = {M_root_pred:.2f} N\u00b7m\")\nprint(f\"    \u03c3_max    (M c / I)               = {sigma_max_pred / 1e6:.2f} MPa\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 4 \u2014 run the static solve\n\nFor a linear-elastic model with the BCs and loads now fully\nspecified, :meth:`Model.solve <femorph_solver.Model.solve>`\nassembles the global stiffness and force, eliminates the\nconstrained DOFs (:doc:`/reference/theory/bc_elimination`), and\ndispatches to the registered linear backend (Pardiso / CHOLMOD\n/ MUMPS / SuperLU \u2014 see\n:doc:`/reference/solvers/linear_backends`).  The default chain\npicks the fastest backend installed on your system.\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": [
        "## Step 5 \u2014 recover and inspect the stress field\n\nDisplacement-method FEM produces nodal displacements as the\nprimary unknown.  Stress is a *derived* quantity recovered by\ndifferentiating the displacement field via the element's\nstrain-displacement matrix $\\mathbf{B}$, then applying\nthe elasticity matrix $\\mathbf{C}$ (Cook \u00a76.14).\n\n:func:`femorph_solver.recover.compute_nodal_stress` does this\nat every grid node, averaging across every element that\ntouches the node.\n\n:func:`femorph_solver.recover.stress_invariants` then derives\nthe standard scalar fields (von Mises, principal stresses,\nhydrostatic, max shear) from the per-node Voigt array.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "u = np.asarray(res.displacement, dtype=np.float64).ravel()\nsigma = compute_nodal_stress(m, u)\ninv = stress_invariants(sigma)\n\n# Mid-span, top-fibre is where the bending stress peaks for this\n# loading.  Read the value off the recovered field.\n# Bending in the xz-plane peaks \u03c3_xx at the +z (top) extreme fibre\n# on the root face \u2014 that's where the largest bending tension lives.\nroot_top_idx = np.where((pts[:, 0] < 1e-9) & (pts[:, 2] > HEIGHT - 1e-9))[0]\nsigma_xx_at_root = float(sigma[root_top_idx, 0].mean())  # \u03c3_xx is bending\nsigma_vm_at_root = float(inv[\"von_mises\"][root_top_idx].mean())\n\n# Tip deflection sampled at a tip-face corner; read the z-component\n# (the bending direction in the xz plane).\ntip_corner = tip_face[(pts[tip_face, 1] < 1e-9) & (pts[tip_face, 2] < 1e-9)][0]\ndelta_tip_uz = float(np.asarray(res.displacement).reshape(-1, 3)[tip_corner, 2])\n\nprint()\nprint(\n    f\"  \u03b4_tip computed  = {delta_tip_uz * 1e3:+.4f} mm  (closed form: {-delta_pred * 1e3:+.4f} mm)\"\n)\nerr_delta = (abs(delta_tip_uz) - delta_pred) / delta_pred * 100.0\nprint(f\"  \u03b4_tip rel error vs 1D Euler-Bernoulli = {err_delta:+.3f} %\")\n# A few % drift below the 1-D Bernoulli closed form is expected on a\n# 3-D HEX8-EAS mesh: the solid model captures Poisson cross-section\n# coupling that the plane-sections-remain-plane assumption ignores.\n# Convergence on the mesh-refinement ladder below tells us whether\n# the FE answer has stabilised.\n\nprint()\nprint(\n    f\"  \u03c3_xx at root, top fibre  = {sigma_xx_at_root / 1e6:+.2f} MPa  \"\n    f\"(closed form \u2248 {sigma_max_pred / 1e6:+.2f} MPa)\"\n)\nprint(f\"  \u03c3_VM at root, top fibre  = {sigma_vm_at_root / 1e6:+.2f} MPa\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 6 \u2014 engineering check against the allowable\n\nWith $\\sigma_\\text{max}$ recovered, the design check is a\nsimple ratio.  We use $\\sigma_\\text{VM}$ (von Mises) as\nthe consolidated failure metric \u2014 it accounts for biaxial\nstress contributions automatically, even though they're small\nin pure bending.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "safety_factor = sigma_allowable / sigma_vm_at_root\nprint()\nprint(\n    f\"  Safety factor: \u03c3_allowable / \u03c3_VM_max = \"\n    f\"{sigma_allowable / 1e6:.0f} MPa / {sigma_vm_at_root / 1e6:.2f} MPa = {safety_factor:.2f}\"\n)\n\nif safety_factor >= 1.5:\n    print(f\"  PASS: SF = {safety_factor:.2f} \u2265 1.5 design target.\")\nelse:\n    print(f\"  FAIL: SF = {safety_factor:.2f} < 1.5 design target \u2014 increase section.\")\nassert safety_factor >= 1.5, \"design fails the 1.5 SF check\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 7 \u2014 mesh-refinement convergence study\n\nA single-mesh result is never enough.  Even on a problem with\na clean closed form, the recovered stress depends on the mesh\n\u2014 the 3-D solid captures Poisson contributions that the 1-D\nclosed form ignores, and the stress recovery itself smooths\nacross element boundaries.  A converged answer is one whose\nvalue stops changing as we refine.\n\nWe sweep three meshes (coarse / medium / fine) and tabulate\nthe recovered $\\sigma_\\text{xx}$ at the root-top fibre.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print()\nprint(f\"  {'mesh':>16}  {'\u03c3_xx at root [MPa]':>20}  {'\u03b4 tip [mm]':>11}  {'safety factor':>14}\")\nprint(\"  \" + \"-\" * 70)\n\n\ndef recover_at_root(model: femorph_solver.Model, result):\n    pts_local = np.asarray(model.grid.points)\n    u_local = np.asarray(result.displacement, dtype=np.float64).ravel()\n    sigma_local = compute_nodal_stress(model, u_local)\n    inv_local = stress_invariants(sigma_local)\n    root_top_local = np.where((pts_local[:, 0] < 1e-9) & (pts_local[:, 2] > HEIGHT - 1e-9))[0]\n    sxx = float(sigma_local[root_top_local, 0].mean())\n    vm = float(inv_local[\"von_mises\"][root_top_local].mean())\n    tip_face_local = np.where(pts_local[:, 0] > L - 1e-9)[0]\n    tip_corner_local = tip_face_local[\n        (pts_local[tip_face_local, 1] < 1e-9) & (pts_local[tip_face_local, 2] < 1e-9)\n    ][0]\n    delta_local = float(np.asarray(result.displacement).reshape(-1, 3)[tip_corner_local, 2])\n    return sxx, vm, delta_local\n\n\nfor tag, (nx_r, ny_r, nz_r) in (\n    (\"coarse  20\u00d73\u00d73\", (20, 3, 3)),\n    (\"medium  40\u00d73\u00d73\", (40, 3, 3)),\n    (\"fine    80\u00d73\u00d73\", (80, 3, 3)),\n):\n    m_r = build_cantilever(nx=nx_r, ny=ny_r, nz=nz_r)\n    pts_r = np.asarray(m_r.grid.points)\n    nn_r = np.asarray(m_r.grid.point_data[\"ansys_node_num\"], dtype=np.int64)\n    clamped_r = np.where(pts_r[:, 0] < 1e-9)[0]\n    m_r.fix(nodes=nn_r[clamped_r].tolist(), dof=\"ALL\")\n\n    tip_face_r = np.where(pts_r[:, 0] > L - 1e-9)[0]\n    tip_top_r = tip_face_r[pts_r[tip_face_r, 2] > HEIGHT - 1e-9]\n    tip_bot_r = tip_face_r[pts_r[tip_face_r, 2] < 1e-9]\n    fz_r = P_tip / float(tip_face_r.size)\n    for n in tip_face_r:\n        m_r.apply_force(int(nn_r[n]), fz=fz_r, accumulate=True)\n    couple_r = M_y_tip / HEIGHT\n    for n in tip_top_r:\n        m_r.apply_force(int(nn_r[n]), fx=+couple_r / float(tip_top_r.size))\n    for n in tip_bot_r:\n        m_r.apply_force(int(nn_r[n]), fx=-couple_r / float(tip_bot_r.size))\n    top_r = np.where(pts_r[:, 2] > HEIGHT - 1e-9)[0]\n    dx_r = L / nx_r\n    for n in top_r:\n        x = pts_r[n, 0]\n        col_w = 0.5 if (x < 1e-9 or x > L - 1e-9) else 1.0\n        y = pts_r[n, 1]\n        y_w = 0.5 if (y < 1e-9 or y > WIDTH - 1e-9) else 1.0\n        fz = (q_udl * dx_r * col_w * y_w) / ny_r\n        m_r.apply_force(int(nn_r[n]), fz=fz, accumulate=True)\n\n    res_r = m_r.solve_static()\n    sxx_r, vm_r, delta_r = recover_at_root(m_r, res_r)\n    sf_r = sigma_allowable / vm_r\n    print(f\"  {tag:>16}  {sxx_r / 1e6:>+18.3f}  {delta_r * 1e3:>+9.3f}  {sf_r:>13.2f}\")\n\nprint()\nprint(\n    f\"  Closed form (EB):   \u03c3_xx \u2248 {-sigma_max_pred / 1e6:+.2f} MPa, \"\n    f\"\u03b4 \u2248 {-delta_pred * 1e3:+.3f} mm\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Take-aways\n\nYou just walked the full design-by-analysis loop on a 3-D solid\ncantilever:\n\n- **Geometry and material** are 1-line operations on a\n  ``Model`` built from a pyvista grid.\n- **Boundary conditions** use coordinate-based masks plus\n  ``Model.fix(nodes=..., dof=\"ALL\")`` for clamps.\n- **Combined loading** decomposes the design load case into\n  nodal-force components \u2014 moments via the mechanical couple\n  form, distributed loads via the trapezoid-rule convention.\n- **The static solve** is a single ``Model.solve_static()`` call;\n  the linear-backend chain picks the fastest direct factor\n  installed on the system.\n- **Stress recovery** runs through the public\n  ``compute_nodal_stress`` + ``stress_invariants`` pair \u2014\n  covered in detail by the post-processing recipes\n  `sphx_glr_gallery_post-processing_example_nodal_stress_recovery.py`\n  and\n  `sphx_glr_gallery_post-processing_example_principal_stress.py`.\n- **The design check** reduces to a safety-factor ratio against\n  the material allowable.\n- **Mesh-refinement convergence** is non-negotiable on any 3-D\n  solid model \u2014 single-mesh stress is suspect until you've\n  shown the value stops changing under refinement.\n\nNext steps:\n\n- Tutorial 2 (planned) \u2014 modal survey of a bracket and\n  response-spectrum analysis.\n- Tutorial 3 (planned) \u2014 pressure vessel design-by-analysis\n  (Lam\u00e9 / thick cylinder).\n\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
}