{ "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 shape-function surfaces on the reference cube\n\nRenders all eight trilinear Lagrange shape functions\n$N_i(\\xi, \\eta, \\zeta)$ for the 8-node hexahedron on\nthe reference cube $\\hat\\Omega = [-1, +1]^3$. The\nbasis functions are\n\n\\begin{align}N_i(\\xi, \\eta, \\zeta) = \\tfrac{1}{8}\n (1 + \\xi_i \\xi)(1 + \\eta_i \\eta)(1 + \\zeta_i \\zeta),\n \\qquad i = 1, \\ldots, 8,\\end{align}\n\nwith $(\\xi_i, \\eta_i, \\zeta_i) \\in \\{-1, +1\\}^3$ the\nnodal coordinates. Each $N_i$ is a tensor product of\n1-D linear \"hat\" functions, equal to 1 at corner $i$\nand 0 at every other corner \u2014 the Kronecker-delta\ninterpolation property.\n\nThe plot below shows $N_i$ evaluated on the\n$\\zeta = -1$ cube face ($i = 1, \\ldots, 4$) and\nthe $\\zeta = +1$ face ($i = 5, \\ldots, 8$).\n\n## References\n\n* Hughes, T. J. R. (2000) *The Finite Element Method \u2014\n Linear Static and Dynamic Finite Element Analysis*, Dover,\n \u00a73.1 \u2014 trilinear shape functions.\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, \u00a76.2.\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 matplotlib.pyplot as plt\nimport numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Corner natural-coordinate map (matches VTK_HEXAHEDRON ordering)\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)\n\n\ndef n_i(xi, eta, zeta, i):\n \"\"\"Trilinear Lagrange shape function for corner i.\"\"\"\n xi_i, eta_i, zeta_i = corners[i]\n return 0.125 * (1 + xi_i * xi) * (1 + eta_i * eta) * (1 + zeta_i * zeta)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sample on a 41\u00d741 grid on each \u03b6 = \u00b11 face\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n = 41\nxi = eta = np.linspace(-1.0, 1.0, n)\nXI, ETA = np.meshgrid(xi, eta, indexing=\"ij\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "4\u00d72 surface plot \u2014 one column per \u03b6-face, 4 rows for the 4 corners\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axes = plt.subplots(\n 4,\n 2,\n figsize=(10.0, 14.0),\n subplot_kw={\"projection\": \"3d\"},\n constrained_layout=True,\n)\nfig.suptitle(\"HEX8 shape functions $N_i(\\\\xi, \\\\eta, \\\\zeta)$ on cube faces\", fontsize=13)\n\nface_zeta = (-1.0, +1.0)\nnode_indices_per_face = ((0, 1, 2, 3), (4, 5, 6, 7))\nfor col, (zeta, idxs) in enumerate(zip(face_zeta, node_indices_per_face, strict=True)):\n for row, idx in enumerate(idxs):\n ax = axes[row, col]\n N = n_i(XI, ETA, zeta, idx)\n ax.plot_surface(XI, ETA, N, cmap=\"plasma\", linewidth=0, antialiased=True)\n x_i, y_i, _ = corners[idx]\n ax.set_title(\n f\"$N_{{{idx + 1}}}$ at corner $({int(x_i):+d}, {int(y_i):+d}, {int(zeta):+d})$\",\n fontsize=10,\n )\n ax.set_xlabel(r\"$\\xi$\")\n ax.set_ylabel(r\"$\\eta$\")\n ax.set_zlabel(r\"$N$\")\n ax.set_xticks([-1, 0, 1])\n ax.set_yticks([-1, 0, 1])\n ax.set_zticks([0, 0.5, 1])\n ax.view_init(elev=22, azim=-55)\n\nfig.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Verify partition of unity at every interior sample point\n\n$\\sum_{i=1}^8 N_i(\\xi, \\eta, \\zeta) = 1$ for any\n$(\\xi, \\eta, \\zeta) \\in \\hat\\Omega$. Sample on a\n21\u00d721\u00d721 lattice and assert.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n3 = 21\nxi3 = eta3 = zeta3 = np.linspace(-1.0, 1.0, n3)\nXI3, ETA3, ZETA3 = np.meshgrid(xi3, eta3, zeta3, indexing=\"ij\")\ntotal = sum(n_i(XI3, ETA3, ZETA3, idx) for idx in range(8))\nmax_dev = float(np.max(np.abs(total - 1.0)))\nassert max_dev < 1.0e-14, f\"partition-of-unity violated: max |\u03a3 N_i - 1| = {max_dev:.3e}\"\nprint(f\"max |\u03a3 N_i - 1| over a 21^3 grid = {max_dev:.3e}\")\nprint(\"OK \u2014 eight trilinear HEX8 shape functions form a partition of unity.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Verify Kronecker-delta interpolation at every corner\n\n$N_i(\\boldsymbol{\\xi}_j) = \\delta_{ij}$: each shape\nfunction equals 1 at its own corner and 0 at all other\ncorners \u2014 required so the nodal-displacement DOFs are\nliterally the displacement values.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print()\nfor i in range(8):\n row = [n_i(*corners[j], i) for j in range(8)]\n expected = [1.0 if j == i else 0.0 for j in range(8)]\n np.testing.assert_allclose(row, expected, atol=1.0e-15)\nprint(\"OK \u2014 Kronecker-delta interpolation verified at every corner.\")" ] } ], "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 }