{
  "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# HEX20 \u2014 uniaxial tension on a 20-node hex\n\nSingle HEX20 brick (a unit cube, 20 nodes including the mid-edge\nnodes of the serendipity family) loaded in uniaxial tension.  The\nconsistent-load vector for a serendipity 8-node face has corner weight\n$-1/12$ and mid-edge weight $+1/3$; applying these produces\na uniform \u03c3xx field, so ``\u03b5xx = \u03c3/E`` and ``\u03b5yy = \u03b5zz = \u2212\u03bd\u00b7\u03b5xx``\nexactly at every node.\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_QUADRATIC_HEXAHEDRON\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Reference 20-node unit cube\nCorners 1-8 in VTK_QUADRATIC_HEXAHEDRON order, then mid-edge nodes 9-20 on the\nbottom, top, and vertical edges.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "E = 2.1e11  # Pa\nNU = 0.30\nF_TOTAL = 4.2e4  # N\n\ncoords = np.array(\n    [\n        [0.0, 0.0, 0.0],\n        [1.0, 0.0, 0.0],\n        [1.0, 1.0, 0.0],\n        [0.0, 1.0, 0.0],\n        [0.0, 0.0, 1.0],\n        [1.0, 0.0, 1.0],\n        [1.0, 1.0, 1.0],\n        [0.0, 1.0, 1.0],\n        [0.5, 0.0, 0.0],\n        [1.0, 0.5, 0.0],\n        [0.5, 1.0, 0.0],\n        [0.0, 0.5, 0.0],\n        [0.5, 0.0, 1.0],\n        [1.0, 0.5, 1.0],\n        [0.5, 1.0, 1.0],\n        [0.0, 0.5, 1.0],\n        [0.0, 0.0, 0.5],\n        [1.0, 0.0, 0.5],\n        [1.0, 1.0, 0.5],\n        [0.0, 1.0, 0.5],\n    ],\n    dtype=np.float64,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build the model\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cells = np.concatenate([[20], np.arange(20, dtype=np.int64)])\ncell_types = np.array([VTK_QUADRATIC_HEXAHEDRON], dtype=np.uint8)\ngrid = pv.UnstructuredGrid(cells, cell_types, coords)\n\nm = femorph_solver.Model.from_grid(grid)\nm.assign(ELEMENTS.HEX20, material={\"EX\": E, \"PRXY\": NU})\n\n# Symmetry BC: x=0 face UX, y=0 face UY, z=0 face UZ.\nx0 = [1, 4, 5, 8, 12, 16, 17, 20]\ny0 = [1, 2, 5, 6, 9, 13, 17, 18]\nz0 = [1, 2, 3, 4, 9, 10, 11, 12]\nfor nn in x0:\n    m.fix(nodes=[nn], dof=\"UX\", value=0.0)\nfor nn in y0:\n    m.fix(nodes=[nn], dof=\"UY\", value=0.0)\nfor nn in z0:\n    m.fix(nodes=[nn], dof=\"UZ\", value=0.0)\n\n# Consistent-load on x=1 face: 4 corners \u00d7 \u2212F/12 + 4 mid-edges \u00d7 +F/3\n# (integrates to F_TOTAL exactly for a uniform traction).\nfor nn in (2, 3, 6, 7):\n    m.apply_force(nn, fx=-F_TOTAL / 12.0)\nfor nn in (10, 14, 18, 19):\n    m.apply_force(nn, fx=F_TOTAL / 3.0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Static solve and strain recovery\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()\neps = res.elastic_strain(model=m)\n\neps_xx_expected = F_TOTAL / E\neps_yy_expected = -NU * eps_xx_expected\nprint(f\"\u03b5xx expected = {eps_xx_expected:.3e}  |  recovered (mean) = {eps[:, 0].mean():.3e}\")\nprint(f\"\u03b5yy expected = {eps_yy_expected:.3e}  |  recovered (mean) = {eps[:, 1].mean():.3e}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot the deformed cube, coloured by \u03b5xx\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "grid = femorph_solver.io.static_result_to_grid(m, res)\ngrid.point_data[\"eps_xx\"] = eps[:, 0]\n\nplotter = pv.Plotter(off_screen=True)\nplotter.add_mesh(grid, style=\"wireframe\", color=\"gray\", line_width=1, opacity=0.4)\nplotter.add_mesh(\n    grid.warp_by_vector(\"displacement\", factor=2.0e5),\n    scalars=\"eps_xx\",\n    show_edges=True,\n    cmap=\"viridis\",\n    scalar_bar_args={\"title\": \"\u03b5xx\"},\n)\nplotter.add_axes()\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
}