{
  "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# HEX20 shape-function slices on the bottom face\n\nRenders the eight HEX20 serendipity shape functions that do\n*not* vanish on the $\\zeta = -1$ face of the reference\ncube $\\hat\\Omega = [-1, +1]^3$: four **corner** functions\n$N^c_i$ whose nodes lie on that face, plus four **mid-edge**\nfunctions $N^m_j$ whose mid-edge nodes lie on the same\nface.  The other twelve basis functions (top-face corners,\ntop-face mid-edges, and the four vertical-edge mid-edges with\n$\\zeta_j = 0$) all vanish at $\\zeta = -1$ and are\nomitted from the figure.\n\n## Serendipity basis\n\nIn reference coordinates $(\\xi, \\eta, \\zeta) \\in [-1, +1]^3$,\nwith corner-node signs\n$(\\xi_i, \\eta_i, \\zeta_i) \\in \\{-1, +1\\}^3$,\n\n\\begin{align}N^c_i(\\xi, \\eta, \\zeta) = \\tfrac{1}{8}\n      (1 + \\xi_i\\xi)(1 + \\eta_i\\eta)(1 + \\zeta_i\\zeta)\n      (\\xi_i\\xi + \\eta_i\\eta + \\zeta_i\\zeta - 2),\\end{align}\n\nFor mid-edges with $\\xi_j = 0$,\n\n\\begin{align}N^m_j(\\xi, \\eta, \\zeta) = \\tfrac{1}{4}(1 - \\xi^{2})\n      (1 + \\eta_j\\eta)(1 + \\zeta_j\\zeta),\\end{align}\n\nwith analogous forms for $\\eta_j = 0$ and $\\zeta_j = 0$\nmid-edges (the missing-variable axis carries the\n$(1 - \\cdot^{2})$ factor).\n\nEach corner function equals 1 at its own node and 0 at every\nother node \u2014 corner *and* mid-edge \u2014 and similarly for mid-edge\nfunctions.  The figure visualises this Kronecker-delta property\non the bottom-face slice.\n\n## References\n\n* Ergatoudis, J. G., Irons, B. M. and Zienkiewicz, O. C. (1968)\n  \"Curved isoparametric quadrilateral elements for finite element\n  analysis,\" *International Journal of Solids and Structures*\n  4 (1), 31\u201342.\n* Hughes, T. J. R. (2000) *The Finite Element Method \u2014\n  Linear Static and Dynamic Finite Element Analysis*, Dover,\n  \u00a73.7.\n* Zienkiewicz, O. C., Taylor, R. L. (2013) *The Finite Element\n  Method: Its Basis and Fundamentals*, 7th ed., \u00a78.7.3.\n\nImplementation: :class:`femorph_solver.elements.hex20.Hex20`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from __future__ import annotations\n\nimport itertools\n\nimport matplotlib.pyplot as plt\nimport numpy as np"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Corner-node signs and mid-edge node coordinates\n\nCorner ordering follows the VTK_QUADRATIC_HEXAHEDRON convention: bottom face\n(zeta=-1) ccw, then top face (zeta=+1) ccw.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "corner_signs = 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)\n\n# Mid-edge nodes: 12 entries, each with one zero coordinate (the\n# axis along which the edge runs).  Order follows VTK_QUADRATIC_HEXAHEDRON:\n# four bottom-face mid-edges, four top-face mid-edges, four\n# vertical mid-edges.\nmid_edge_signs = np.array(\n    [\n        [0, -1, -1],  # 8  edge 0-1 (xi-axis, eta=-1, zeta=-1)\n        [+1, 0, -1],  # 9  edge 1-2\n        [0, +1, -1],  # 10 edge 2-3\n        [-1, 0, -1],  # 11 edge 3-0\n        [0, -1, +1],  # 12 edge 4-5 (top face)\n        [+1, 0, +1],  # 13\n        [0, +1, +1],  # 14\n        [-1, 0, +1],  # 15\n        [-1, -1, 0],  # 16 vertical edge 0-4\n        [+1, -1, 0],  # 17\n        [+1, +1, 0],  # 18\n        [-1, +1, 0],  # 19\n    ],\n    dtype=float,\n)\n\n\ndef n_corner(xi, eta, zeta, signs):\n    \"\"\"Serendipity corner shape function for a given (xi_i, eta_i, zeta_i).\"\"\"\n    xi_i, eta_i, zeta_i = signs\n    return (\n        0.125\n        * (1 + xi_i * xi)\n        * (1 + eta_i * eta)\n        * (1 + zeta_i * zeta)\n        * (xi_i * xi + eta_i * eta + zeta_i * zeta - 2)\n    )\n\n\ndef n_mid_edge(xi, eta, zeta, signs):\n    \"\"\"Serendipity mid-edge shape function: the zero-axis carries (1 - .^2).\"\"\"\n    xi_i, eta_i, zeta_i = signs\n    if xi_i == 0:\n        return 0.25 * (1 - xi**2) * (1 + eta_i * eta) * (1 + zeta_i * zeta)\n    if eta_i == 0:\n        return 0.25 * (1 - eta**2) * (1 + xi_i * xi) * (1 + zeta_i * zeta)\n    return 0.25 * (1 - zeta**2) * (1 + xi_i * xi) * (1 + eta_i * eta)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Sample the bottom face zeta = -1 on a (xi, eta) grid\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n = 51\nxi_grid, eta_grid = np.meshgrid(np.linspace(-1, 1, n), np.linspace(-1, 1, n), indexing=\"ij\")\nzeta_slice = -1.0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Render \u2014 4\u00d72 grid of contour plots\n\nTop row: the four corner shape functions whose nodes lie on\nthe bottom face (corners 0..3, all with zeta_i = -1).\nBottom row: the four mid-edge functions whose nodes also lie\non the bottom face (mid-edges 8..11).  The other twelve basis\nfunctions vanish identically on this slice.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "bottom_face_corners = [0, 1, 2, 3]\nbottom_face_mid_edges = [8, 9, 10, 11]\n\nfig, axes = plt.subplots(2, 4, figsize=(15, 7), constrained_layout=True)\nfig.suptitle(\n    \"HEX20 shape functions on the $\\\\zeta = -1$ face slice\",\n    fontsize=13,\n)\n\nfor col, ci in enumerate(bottom_face_corners):\n    ax = axes[0, col]\n    N = n_corner(xi_grid, eta_grid, zeta_slice, corner_signs[ci])\n    img = ax.contourf(xi_grid, eta_grid, N, levels=21, cmap=\"plasma\")\n    ax.set_aspect(\"equal\")\n    sx, sy, _ = corner_signs[ci].astype(int)\n    ax.set_title(f\"$N^c_{{{ci}}}$  $(\\\\xi_i,\\\\eta_i,\\\\zeta_i)=({sx:+d},{sy:+d},-1)$\", fontsize=9)\n    ax.set_xlabel(r\"$\\xi$\")\n    ax.set_ylabel(r\"$\\eta$\")\n    fig.colorbar(img, ax=ax, shrink=0.8)\n\nfor col, mj in enumerate(bottom_face_mid_edges):\n    ax = axes[1, col]\n    N = n_mid_edge(xi_grid, eta_grid, zeta_slice, mid_edge_signs[mj - 8])\n    img = ax.contourf(xi_grid, eta_grid, N, levels=21, cmap=\"plasma\")\n    ax.set_aspect(\"equal\")\n    sx, sy, _ = mid_edge_signs[mj - 8].astype(int)\n    ax.set_title(f\"$N^m_{{{mj}}}$  $(\\\\xi_j,\\\\eta_j,\\\\zeta_j)=({sx:+d},{sy:+d},-1)$\", fontsize=9)\n    ax.set_xlabel(r\"$\\xi$\")\n    ax.set_ylabel(r\"$\\eta$\")\n    fig.colorbar(img, ax=ax, shrink=0.8)\n\nfig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Verify partition of unity at every reference-element node\n\nBuild the 20 nodal coordinates and check\n$\\sum_i N^c_i + \\sum_j N^m_j \\equiv 1$ to machine precision.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "all_nodes = np.vstack([corner_signs, mid_edge_signs])  # (20, 3)\ntotal = np.zeros(20)\nfor k in range(20):\n    xi_k, eta_k, zeta_k = all_nodes[k]\n    total[k] = sum(n_corner(xi_k, eta_k, zeta_k, corner_signs[i]) for i in range(8)) + sum(\n        n_mid_edge(xi_k, eta_k, zeta_k, mid_edge_signs[j]) for j in range(12)\n    )\nnp.testing.assert_allclose(total, 1.0, atol=1e-14)\nprint(\"OK \u2014 HEX20 partition-of-unity holds at every reference-element node.\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Verify Kronecker-delta interpolation\n\nThe (20, 20) basis matrix evaluated at every reference-element\nnode should be the identity.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "basis = np.zeros((20, 20))\nfor k, (xi_k, eta_k, zeta_k) in enumerate(all_nodes):\n    for i in range(8):\n        basis[k, i] = n_corner(xi_k, eta_k, zeta_k, corner_signs[i])\n    for j in range(12):\n        basis[k, 8 + j] = n_mid_edge(xi_k, eta_k, zeta_k, mid_edge_signs[j])\nnp.testing.assert_allclose(basis, np.eye(20), atol=1e-14)\nprint(\"OK \u2014 HEX20 Kronecker-delta interpolation verified at every node.\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Verify the bottom-face vanishing pattern\n\nOn the zeta = -1 slice, all four top-face corners (4..7), all\nfour top-face mid-edges (12..15), and all four vertical mid-edges\n(16..19) must be identically zero.  The four bottom-face corners\n(0..3) and four bottom-face mid-edges (8..11) carry every value.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "xi_test, eta_test = np.meshgrid(np.linspace(-1, 1, 5), np.linspace(-1, 1, 5), indexing=\"ij\")\nfor i in itertools.chain(range(4, 8)):\n    Ni = n_corner(xi_test, eta_test, -1.0, corner_signs[i])\n    np.testing.assert_allclose(Ni, 0.0, atol=1e-14)\nfor j in itertools.chain(range(4, 12)):\n    Nj = n_mid_edge(xi_test, eta_test, -1.0, mid_edge_signs[j])\n    np.testing.assert_allclose(Nj, 0.0, atol=1e-14)\nprint(\"OK \u2014 twelve off-face HEX20 basis functions vanish identically on zeta = -1.\")"
      ]
    }
  ],
  "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
}