{
  "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\n# HEX8 \u2014 cantilever plate static analysis\n\nStatic analysis of a cantilever plate under a distributed tip load.\nThis example walks through the full\n:meth:`femorph_solver.Model.solve` \u2192 :class:`StaticResult` \u2192\n:func:`femorph_solver.io.static_result_to_grid` \u2192 pyvista rendering loop and\nchecks static equilibrium via :attr:`StaticResult.reaction`.\n\nEuler\u2013Bernoulli beam theory is included as a back-of-envelope\nreference; HEX8 is a first-order hex with full 2 \u00d7 2 \u00d7 2 Gauss\nintegration, which exhibits well-known shear locking in thin-bending\nproblems unless many elements are used through the thickness.  The\ndifference between the two is a feature of the element, not a bug \u2014\nswap in ``HEX20`` (quadratic) to recover EB to a few percent with\nthe same mesh.\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_HEXAHEDRON\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Geometry + material\nSteel cantilever, 1 m \u00d7 0.1 m \u00d7 0.05 m, meshed 40 \u00d7 4 \u00d7 4 hex\n(640 HEX8 elements, 1 025 nodes).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "E = 2.0e11  # Pa\nNU = 0.30\nRHO = 7850.0\nLX, LY, LZ = 1.0, 0.1, 0.05\nNX, NY, NZ = 40, 4, 4\nF_TIP = -5.0e3  # N (downward)\n\nxs = np.linspace(0.0, LX, NX + 1)\nys = np.linspace(0.0, LY, NY + 1)\nzs = np.linspace(0.0, LZ, NZ + 1)\nxx, yy, zz = np.meshgrid(xs, ys, zs, indexing=\"ij\")\npoints = np.stack([xx.ravel(), yy.ravel(), zz.ravel()], axis=1)\n\n\ndef _node_idx(i: int, j: int, k: int) -> int:\n    \"\"\"0-based VTK point index for the structured mesh.\"\"\"\n    return (i * (NY + 1) + j) * (NZ + 1) + k\n\n\ncells_flat: list[int] = []\nfor i in range(NX):\n    for j in range(NY):\n        for k in range(NZ):\n            cells_flat.extend(\n                [\n                    8,\n                    _node_idx(i, j, k),\n                    _node_idx(i + 1, j, k),\n                    _node_idx(i + 1, j + 1, k),\n                    _node_idx(i, j + 1, k),\n                    _node_idx(i, j, k + 1),\n                    _node_idx(i + 1, j, k + 1),\n                    _node_idx(i + 1, j + 1, k + 1),\n                    _node_idx(i, j + 1, k + 1),\n                ]\n            )\n\nn_cells = NX * NY * NZ\ncell_types = np.full(n_cells, VTK_HEXAHEDRON, dtype=np.uint8)\ngrid = pv.UnstructuredGrid(np.asarray(cells_flat, dtype=np.int64), cell_types, points)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build the model\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m = femorph_solver.Model.from_grid(grid)\nm.assign(ELEMENTS.HEX8, material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO})\n\nnode_nums = np.asarray(m.grid.point_data[\"ansys_node_num\"])\npts = np.asarray(m.grid.points)\n\n# Clamp the ``x = 0`` face in all 3 DOFs.\nx0_mask = pts[:, 0] < 1e-9\nx0_nodes = node_nums[x0_mask].tolist()\nm.fix(nodes=x0_nodes, dof=\"UX\")\nm.fix(nodes=x0_nodes, dof=\"UY\")\nm.fix(nodes=x0_nodes, dof=\"UZ\")\n\n# Distributed downward tip load.\ntip_mask = pts[:, 0] > LX - 1e-9\ntip_nodes = node_nums[tip_mask].tolist()\nfz_each = F_TIP / len(tip_nodes)\nfor nn in tip_nodes:\n    m.apply_force(int(nn), fz=fz_each)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Solve + reaction check\n:meth:`Model.solve` returns a :class:`StaticResult` with\n``displacement``, ``reaction``, and ``free_mask``.  Reactions are\nnonzero only at constrained DOFs; summing ``FZ`` at the clamp must\nequal ``-F_TIP`` to machine precision for a well-posed static solve.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()\n\ndof = m.dof_map()\nfz_clamp = 0.0\nfor nn in x0_nodes:\n    rows = np.where((dof[:, 0] == nn) & (dof[:, 1] == 2))[0]\n    for r in rows:\n        fz_clamp += float(res.reaction[r])\nprint(f\"\u03a3 FZ reaction at clamp = {fz_clamp:.4e} N   (expected {-F_TIP:.4e})\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Tip deflection vs Euler\u2013Bernoulli\n$\\\\delta_\\\\mathrm{EB} = F L^3 / (3 E I)$ with\n$I = b h^3 / 12$ is the slender-beam estimate.  With 4 elements\nthrough the thickness, HEX8's shear locking gives a few percent\nerror \u2014 swap in ``HEX20`` (see `ref_solid186_example`) to\nremove it entirely.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "I_y = LY * LZ**3 / 12.0\ndelta_eb = F_TIP * LX**3 / (3.0 * E * I_y)\n\ngrid = femorph_solver.io.static_result_to_grid(m, res)\ntip_mask = grid.points[:, 0] > LX - 1e-9\nw_tip_femorph_solver = grid.point_data[\"displacement\"][tip_mask, 2].min()\nprint(f\"Euler-Bernoulli tip deflection = {delta_eb:.3e} m\")\nprint(f\"femorph-solver tip deflection (min UZ)   = {w_tip_femorph_solver:.3e} m\")\nprint(\n    f\"relative error                  = {abs(w_tip_femorph_solver - delta_eb) / abs(delta_eb):.2%}\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the deformed plate, coloured by displacement magnitude\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "warped = grid.warp_by_vector(\"displacement\", factor=20.0)\n\nplotter = pv.Plotter(off_screen=True)\nplotter.add_mesh(\n    m.grid,\n    style=\"wireframe\",\n    color=\"gray\",\n    opacity=0.35,\n    label=\"undeformed\",\n)\nplotter.add_mesh(\n    warped,\n    scalars=\"displacement_magnitude\",\n    show_edges=True,\n    cmap=\"viridis\",\n    scalar_bar_args={\"title\": \"|u| [m]\"},\n    label=\"deformed \u00d720\",\n)\nplotter.add_legend()\nplotter.add_axes()\nplotter.camera_position = [(2.4, -1.6, 1.0), (0.5, 0.05, 0.0), (0.0, 0.0, 1.0)]\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
}