{
  "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# HEX8 reference geometry \u2014 nodes, shape functions, Gauss points\n\nRenders the 8-node trilinear hexahedron's **reference element**\n$\\hat\\Omega = [-1, 1]^3$ with three layers of annotation:\n\n1. **Corner nodes** with their natural-coordinate triplets\n   $(\\xi_i, \\eta_i, \\zeta_i) \\in \\{-1, +1\\}^3$.\n2. **Trilinear Lagrange shape functions**\n\n   .. math::\n\n      N_i(\\xi, \\eta, \\zeta) =\n        \\tfrac{1}{8}(1 + \\xi_i\\xi)(1 + \\eta_i\\eta)(1 + \\zeta_i\\zeta),\n        \\qquad i = 1, \\ldots, 8,\n\n   evaluated on a sampling lattice and visualised by colouring the\n   reference cube faces with $N_1(\\xi, \\eta, \\zeta)$.\n3. **2\u00d72\u00d72 Gauss-Legendre integration points** at\n   $(\\xi, \\eta, \\zeta) \\in \\{-1, +1\\}/\\sqrt{3}$, used by the\n   stiffness and consistent-mass integrals on this element.\n\n## References\n* Hughes, T. J. R. (2000) *The Finite Element Method \u2014 Linear\n  Static and Dynamic Finite Element Analysis*, Dover, \u00a73.1\n  (trilinear shape functions); \u00a73.10 (Gauss-Legendre quadrature).\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, \u00a76.2.\n* Zienkiewicz, O. C. and Taylor, R. L. (2013) *The Finite Element\n  Method: Its Basis and Fundamentals*, 7th ed., Butterworth-\n  Heinemann, \u00a76.\n\nImplementation: :class:`femorph_solver.elements.hex8.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"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Reference cube as a single-cell pyvista UnstructuredGrid\n\nNative HEX8 corner ordering ``(\u03be, \u03b7, \u03b6) \u2208 {-1, +1}^3`` follows the\nstandard isoparametric convention \u2014 VTK_HEXAHEDRON's node order\nmatches it, so we author the geometry once and visualise directly.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "corners = np.array(\n    [\n        [-1, -1, -1],\n        [+1, -1, -1],\n        [+1, +1, -1],\n        [-1, +1, -1],\n        [-1, -1, +1],\n        [+1, -1, +1],\n        [+1, +1, +1],\n        [-1, +1, +1],\n    ],\n    dtype=float,\n)\ncells = np.hstack([[8], np.arange(8, dtype=np.int64)])\ncell_types = np.array([pv.CellType.HEXAHEDRON], dtype=np.uint8)\nref_cube = pv.UnstructuredGrid(cells, cell_types, corners)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Shape function $N_1$ over a sampling lattice\n\nSample $N_1(\\xi, \\eta, \\zeta) = \\tfrac18 (1 - \\xi)(1 - \\eta)\n(1 - \\zeta)$ on a 21^3 lattice and store as a ``StructuredGrid``\npoint array.  Slicing through the volume gives 2D contour-friendly\nviews.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n = 21\nxi = eta = zeta = np.linspace(-1.0, 1.0, n)\nXI, ETA, ZETA = np.meshgrid(xi, eta, zeta, indexing=\"ij\")\nN1 = (1.0 - XI) * (1.0 - ETA) * (1.0 - ZETA) / 8.0\nsample = pv.StructuredGrid(XI, ETA, ZETA)\nsample[\"N_1\"] = N1.ravel(order=\"F\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 2\u00d72\u00d72 Gauss-Legendre integration points\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g = 1.0 / np.sqrt(3.0)\ngauss = np.array(\n    [(a, b, c) for a in (-g, +g) for b in (-g, +g) for c in (-g, +g)],\n    dtype=float,\n)\ngauss_pd = pv.PolyData(gauss)\ngauss_pd[\"weight\"] = np.full(len(gauss), 1.0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the figure\n\nThree glyphs share the same scene: the reference-cube wireframe\n(black), an $N_1$ slice through the centre (warm colormap),\nand the eight Gauss points as red spheres.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plotter = pv.Plotter(off_screen=True, window_size=(640, 480))\nplotter.add_mesh(\n    ref_cube,\n    style=\"wireframe\",\n    color=\"black\",\n    line_width=2,\n    label=\"reference element\",\n)\nslc = sample.slice_orthogonal(x=0.0, y=0.0, z=0.0)\nplotter.add_mesh(\n    slc,\n    scalars=\"N_1\",\n    cmap=\"plasma\",\n    show_edges=False,\n    opacity=0.85,\n    scalar_bar_args={\"title\": \"N_1(\u03be, \u03b7, \u03b6)\"},\n)\nplotter.add_points(\n    gauss,\n    render_points_as_spheres=True,\n    point_size=20,\n    color=\"#d62728\",\n    label=\"Gauss points (2\u00d72\u00d72)\",\n)\n# Corner-node spheres\nplotter.add_points(\n    corners,\n    render_points_as_spheres=True,\n    point_size=14,\n    color=\"black\",\n    label=\"corner nodes\",\n)\nplotter.add_axes(line_width=4, color=\"black\")\nplotter.view_isometric()\nplotter.camera.zoom(1.1)\nplotter.add_legend(face=None, size=(0.18, 0.10), bcolor=\"white\")\nplotter.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Verify the partition of unity\n\nAt any point inside $\\hat\\Omega$ the eight shape functions\nmust sum to one.  Sampling $N_i(\\xi_q, \\eta_q, \\zeta_q)$ at\nevery Gauss point and asserting $\\sum_i N_i = 1$ is a\nminimal sanity check \u2014 failure would point to a typo in the\nkernel's shape-function evaluation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "xi_i, eta_i, zeta_i = corners.T\nsums = np.array(\n    [\n        sum(\n            (1 + xi_i * a) * (1 + eta_i * b) * (1 + zeta_i * c) / 8.0\n            for a, b, c in [(g[0], g[1], g[2])]\n        ).sum()\n        for g in gauss\n    ]\n)\nprint(\"partition-of-unity \u03a3 N_i at each Gauss point:\")\nfor q, s in zip(gauss, sums, strict=True):\n    print(f\"  \u03be={q[0]:+.4f}  \u03b7={q[1]:+.4f}  \u03b6={q[2]:+.4f}  \u03a3={s:.15f}\")\nnp.testing.assert_allclose(sums, 1.0, atol=1.0e-14)\nprint(\"OK \u2014 all 8 shape functions sum to 1 at every Gauss point.\")"
      ]
    }
  ],
  "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
}