{ "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 reference geometry \u2014 corner + mid-edge nodes, Gauss points\n\nRenders the 20-node serendipity hexahedron's reference element\n$\\hat\\Omega = [-1, 1]^3$:\n\n* 8 **corner nodes** at $(\\xi, \\eta, \\zeta) \\in \\{-1, +1\\}^3$.\n* 12 **mid-edge nodes** at the centre of every cube edge.\n* 14-point **Irons (1971)** quadrature points used for the\n consistent-mass integral on this kernel.\n* 2\u00d72\u00d72 reduced-integration Gauss points\n $(\\xi, \\eta, \\zeta) \\in \\{\\pm 1/\\sqrt{3}\\}^3$ used for\n the stiffness integral by default in femorph-solver\n (the femorph-solver default for HEX20).\n\n## Shape functions\n\nThe serendipity HEX20 basis (Hughes 2000 \u00a73.7; Zienkiewicz &\nTaylor 2013 \u00a76.6) splits the 20 nodal functions into corner\nand mid-edge groups:\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\n\\begin{align}N^m_j(\\xi, \\eta, \\zeta) = \\tfrac{1}{4}\n (1 - \\xi^{2})(1 + \\eta_j\\eta)(1 + \\zeta_j\\zeta)\n \\quad \\text{(for the 4 mid-edges with } \\xi_j = 0\\text{)},\\end{align}\n\nwith analogous expressions for the four mid-edges along each\nof $\\eta_j = 0$ and $\\zeta_j = 0$.\n\n## References\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* 4\n (1), 31\u201342 \u2014 original serendipity basis.\n* Hughes, T. J. R. (2000) *The Finite Element Method*, Dover, \u00a73.7.\n* Irons, B. M. (1971) \"Quadrature rules for brick based finite\n elements,\" *International Journal for Numerical Methods in\n Engineering* 3 (2), 293\u2013294 \u2014 14-point degree-5 rule.\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 numpy as np\nimport pyvista as pv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Node layout\n\n8 corners + 12 mid-edges = 20 nodes. We label corner indices\n0..7 and mid-edge indices 8..19 to match the standard\nisoparametric ordering used by VTK_QUADRATIC_HEXAHEDRON.\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)\nmid_edges = np.array(\n [\n # bottom (z = -1) ring: 8, 9, 10, 11\n [0, -1, -1],\n [+1, 0, -1],\n [0, +1, -1],\n [-1, 0, -1],\n # top (z = +1) ring: 12, 13, 14, 15\n [0, -1, +1],\n [+1, 0, +1],\n [0, +1, +1],\n [-1, 0, +1],\n # vertical mid-edges (z = 0 ring): 16, 17, 18, 19\n [-1, -1, 0],\n [+1, -1, 0],\n [+1, +1, 0],\n [-1, +1, 0],\n ],\n dtype=float,\n)\nall_nodes = np.vstack([corners, mid_edges])\n\ncells = np.hstack([[20], np.arange(20, dtype=np.int64)])\ncell_types = np.array([pv.CellType.QUADRATIC_HEXAHEDRON], dtype=np.uint8)\nref_cube = pv.UnstructuredGrid(cells, cell_types, all_nodes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Quadrature points \u2014 the two rules HEX20 uses\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# 2x2x2 Gauss-Legendre \u2014 stiffness (reduced).\ng = 1.0 / np.sqrt(3.0)\ngauss_8 = 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)\n\n# Irons 14-point \u2014 consistent mass. Six face-centre points at\n# (\u00b11, 0, 0), (0, \u00b11, 0), (0, 0, \u00b11) and 8 corner-octant points at\n# \u00b1sqrt(19/30)\u00b7(1, 1, 1) with sign permutations.\nb_ = float(np.sqrt(19.0 / 30.0))\nfaces = np.array(\n [[+1, 0, 0], [-1, 0, 0], [0, +1, 0], [0, -1, 0], [0, 0, +1], [0, 0, -1]],\n dtype=float,\n)\noctants = np.array(\n [[a * b_, b * b_, c * b_] for a in (-1, +1) for b in (-1, +1) for c in (-1, +1)],\n dtype=float,\n)\nirons_14 = np.vstack([faces, octants])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render\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 HEX20\",\n)\nplotter.add_points(\n corners,\n render_points_as_spheres=True,\n point_size=18,\n color=\"black\",\n label=\"corner nodes (8)\",\n)\nplotter.add_points(\n mid_edges,\n render_points_as_spheres=True,\n point_size=12,\n color=\"#1f77b4\",\n label=\"mid-edge nodes (12)\",\n)\nplotter.add_points(\n gauss_8,\n render_points_as_spheres=True,\n point_size=14,\n color=\"#d62728\",\n label=\"2\u00d72\u00d72 Gauss (stiffness)\",\n)\nplotter.add_points(\n irons_14,\n render_points_as_spheres=True,\n point_size=10,\n color=\"#2ca02c\",\n label=\"Irons 14-pt (mass)\",\n)\nplotter.add_axes(line_width=4, color=\"black\")\nplotter.view_isometric()\nplotter.camera.zoom(1.05)\nplotter.add_legend(face=None, size=(0.24, 0.16), bcolor=\"white\")\nplotter.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sanity: every Irons-14 point is inside $\\hat\\Omega$\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "assert (np.abs(irons_14) <= 1.0 + 1e-12).all(), \"Irons-14 point left the cube\"\nprint(f\"Irons 14-pt rule: {len(irons_14)} points, all inside [-1, 1]^3.\")\nprint(\" 6 face-centre points at \u00b11 along one axis\")\nprint(\" 8 corner-octant points at \u00b1sqrt(19/30)\u00b7(1,1,1) sign-perm\")" ] } ], "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 }