{
  "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# TET10 shape-function slices on the unit tetrahedron\n\nRenders the ten quadratic-tetrahedron shape functions\n$N_i^c(L)$ (4 corners) and $N_{ij}^m(L)$ (6\nmid-edges) sliced through the $\\xi_3 = 0$ face of the\nunit tetrahedron.  In volume coordinates $L_i$ (where\n$L_0 = 1 - L_1 - L_2 - L_3$) the basis is\n\n\\begin{align}N_i^c = L_i\\, (2 L_i - 1), \\qquad i = 0, \\ldots, 3,\\end{align}\n\n\\begin{align}N_{ij}^m = 4\\, L_i\\, L_j, \\qquad\n    (i, j) \\in \\{(0,1), (1,2), (2,0), (0,3), (1,3), (2,3)\\}.\\end{align}\n\nEach corner shape function vanishes at every other corner and\nevery mid-edge node *not* on its incident edges; each\nmid-edge function equals 1 at its own mid-edge node and 0\neverywhere else.  The figure below visualises this property\non a 2D cross-section.\n\n## References\n\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J.\n  (2002) *Concepts and Applications of Finite Element\n  Analysis*, 4th ed., Wiley, Table 6.5-1.\n* Hughes, T. J. R. (2000) *The Finite Element Method \u2014\n  Linear Static and Dynamic Finite Element Analysis*, Dover,\n  \u00a73.8.\n* Zienkiewicz, O. C., Taylor, R. L. (2013) *The Finite\n  Element Method*, 7th ed., \u00a78.8.2.\n\nImplementation: :class:`femorph_solver.elements.tet10.Tet10`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from __future__ import annotations\n\nimport matplotlib.pyplot as plt\nimport numpy as np"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Unit-tet corner volume coordinates.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "corner_coords = np.eye(4)  # L_0 = (1,0,0,0); L_1 = (0,1,0,0); ...\nedge_pairs = [(0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)]\n\n\ndef n_corner(L, i):\n    \"\"\"Quadratic corner shape function L_i (2 L_i - 1).\"\"\"\n    return L[..., i] * (2 * L[..., i] - 1)\n\n\ndef n_mid_edge(L, i, j):\n    \"\"\"Quadratic mid-edge shape function 4 L_i L_j.\"\"\"\n    return 4 * L[..., i] * L[..., j]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Sample a triangular slice through the xi_3 = 0 face\n\nThat face is the triangle with corners (L_0, L_1, L_2, L_3) =\n(1,0,0,0), (0,1,0,0), (0,0,1,0).  Sample (L_1, L_2) on a\ntriangular grid; L_0 = 1 - L_1 - L_2 fills in.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n = 31\nL1, L2 = np.meshgrid(np.linspace(0, 1, n), np.linspace(0, 1, n), indexing=\"ij\")\nmask = L1 + L2 <= 1\nL0 = 1.0 - L1 - L2\nL3 = np.zeros_like(L0)\n# Stack into shape (n, n, 4); set values outside the triangle to NaN\n# so they don't render.\nL = np.stack([L0, L1, L2, L3], axis=-1)\nL[~mask] = np.nan"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Render \u2014 4\u00d73 grid: rows are corner shape functions (3) +\nthe three mid-edges that lie on the xi_3=0 face (3 more);\nthe six mid-edges off the face (those that include node 3)\nvanish on this slice and aren't worth plotting.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "face_face_pairs = [(0, 1), (1, 2), (2, 0)]\nn_corners_on_face = 3\nfig, axes = plt.subplots(2, 3, figsize=(12, 7), constrained_layout=True)\nfig.suptitle(\n    \"TET10 shape functions on the $\\\\xi_3 = 0$ face slice\",\n    fontsize=13,\n)\n\nfor col in range(3):\n    ax = axes[0, col]\n    N = n_corner(L, col)\n    img = ax.contourf(L1, L2, N, levels=21, cmap=\"plasma\")\n    ax.set_aspect(\"equal\")\n    ax.set_title(f\"$N_{{{col}}}^c = L_{{{col}}} (2 L_{{{col}}} - 1)$\", fontsize=10)\n    ax.set_xlabel(r\"$L_1$\")\n    ax.set_ylabel(r\"$L_2$\")\n    fig.colorbar(img, ax=ax, shrink=0.85)\n\nfor col, (i, j) in enumerate(face_face_pairs):\n    ax = axes[1, col]\n    N = n_mid_edge(L, i, j)\n    img = ax.contourf(L1, L2, N, levels=21, cmap=\"plasma\")\n    ax.set_aspect(\"equal\")\n    ax.set_title(f\"$N_{{{i}{j}}}^m = 4 L_{{{i}}} L_{{{j}}}$\", fontsize=10)\n    ax.set_xlabel(r\"$L_1$\")\n    ax.set_ylabel(r\"$L_2$\")\n    fig.colorbar(img, ax=ax, shrink=0.85)\n\nfig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Verify partition of unity at every node of the tetrahedron\n\nAll 10 nodes \u2014 4 corners + 6 mid-edges \u2014 collected into a\nsingle (10, 4) volume-coordinate array.  At every one,\nsum_i N_i^c + sum_{ij} N_{ij}^m must equal 1.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mid_coords = np.array([0.5 * (corner_coords[i] + corner_coords[j]) for i, j in edge_pairs])\nall_node_coords = np.vstack([corner_coords, mid_coords])  # (10, 4)\ntotal = sum(n_corner(all_node_coords, i) for i in range(4)) + sum(\n    n_mid_edge(all_node_coords, i, j) for i, j in edge_pairs\n)\nnp.testing.assert_allclose(total, 1.0, atol=1e-14)\nprint(\"OK \u2014 TET10 partition-of-unity holds at every reference-element node.\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Verify Kronecker-delta interpolation at every node\n\n- corner i: N_i^c(corner_j) = delta_{ij}; all six mid-edge\n  functions vanish at every corner.\n- mid-edge (i, j): N_{ij}^m(mid_kl) = delta_{(i,j),(k,l)};\n  all four corner functions vanish at every mid-edge.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Stack everything for one big check: the (10, 10) basis matrix\n# evaluated at every reference-element node should be identity.\nn_corner_at_all = np.stack([n_corner(all_node_coords, i) for i in range(4)], axis=-1)\nn_mid_at_all = np.stack([n_mid_edge(all_node_coords, i, j) for i, j in edge_pairs], axis=-1)\nbasis = np.hstack([n_corner_at_all, n_mid_at_all])  # (10, 10)\nnp.testing.assert_allclose(basis, np.eye(10), atol=1e-14)\nprint(\"OK \u2014 TET10 Kronecker-delta interpolation verified at every node.\")"
      ]
    }
  ],
  "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
}