.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/elements/solid186/example_hex20_shape_function_slices.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_elements_solid186_example_hex20_shape_function_slices.py: HEX20 shape-function slices on the bottom face ============================================== Renders the eight HEX20 serendipity shape functions that do *not* vanish on the :math:`\zeta = -1` face of the reference cube :math:`\hat\Omega = [-1, +1]^3`: four **corner** functions :math:`N^c_i` whose nodes lie on that face, plus four **mid-edge** functions :math:`N^m_j` whose mid-edge nodes lie on the same face. The other twelve basis functions (top-face corners, top-face mid-edges, and the four vertical-edge mid-edges with :math:`\zeta_j = 0`) all vanish at :math:`\zeta = -1` and are omitted from the figure. Serendipity basis ----------------- In reference coordinates :math:`(\xi, \eta, \zeta) \in [-1, +1]^3`, with corner-node signs :math:`(\xi_i, \eta_i, \zeta_i) \in \{-1, +1\}^3`, .. math:: N^c_i(\xi, \eta, \zeta) = \tfrac{1}{8} (1 + \xi_i\xi)(1 + \eta_i\eta)(1 + \zeta_i\zeta) (\xi_i\xi + \eta_i\eta + \zeta_i\zeta - 2), For mid-edges with :math:`\xi_j = 0`, .. math:: N^m_j(\xi, \eta, \zeta) = \tfrac{1}{4}(1 - \xi^{2}) (1 + \eta_j\eta)(1 + \zeta_j\zeta), with analogous forms for :math:`\eta_j = 0` and :math:`\zeta_j = 0` mid-edges (the missing-variable axis carries the :math:`(1 - \cdot^{2})` factor). Each corner function equals 1 at its own node and 0 at every other node — corner *and* mid-edge — and similarly for mid-edge functions. The figure visualises this Kronecker-delta property on the bottom-face slice. References ---------- * Ergatoudis, J. G., Irons, B. M. and Zienkiewicz, O. C. (1968) "Curved isoparametric quadrilateral elements for finite element analysis," *International Journal of Solids and Structures* 4 (1), 31–42. * Hughes, T. J. R. (2000) *The Finite Element Method — Linear Static and Dynamic Finite Element Analysis*, Dover, §3.7. * Zienkiewicz, O. C., Taylor, R. L. (2013) *The Finite Element Method: Its Basis and Fundamentals*, 7th ed., §8.7.3. Implementation: :class:`femorph_solver.elements.hex20.Hex20`. .. GENERATED FROM PYTHON SOURCE LINES 59-67 .. code-block:: Python from __future__ import annotations import itertools import matplotlib.pyplot as plt import numpy as np .. GENERATED FROM PYTHON SOURCE LINES 68-72 Corner-node signs and mid-edge node coordinates Corner ordering follows the VTK_QUADRATIC_HEXAHEDRON convention: bottom face (zeta=-1) ccw, then top face (zeta=+1) ccw. .. GENERATED FROM PYTHON SOURCE LINES 72-132 .. code-block:: Python corner_signs = np.array( [ [-1, -1, -1], [+1, -1, -1], [+1, +1, -1], [-1, +1, -1], [-1, -1, +1], [+1, -1, +1], [+1, +1, +1], [-1, +1, +1], ], dtype=float, ) # Mid-edge nodes: 12 entries, each with one zero coordinate (the # axis along which the edge runs). Order follows VTK_QUADRATIC_HEXAHEDRON: # four bottom-face mid-edges, four top-face mid-edges, four # vertical mid-edges. mid_edge_signs = np.array( [ [0, -1, -1], # 8 edge 0-1 (xi-axis, eta=-1, zeta=-1) [+1, 0, -1], # 9 edge 1-2 [0, +1, -1], # 10 edge 2-3 [-1, 0, -1], # 11 edge 3-0 [0, -1, +1], # 12 edge 4-5 (top face) [+1, 0, +1], # 13 [0, +1, +1], # 14 [-1, 0, +1], # 15 [-1, -1, 0], # 16 vertical edge 0-4 [+1, -1, 0], # 17 [+1, +1, 0], # 18 [-1, +1, 0], # 19 ], dtype=float, ) def n_corner(xi, eta, zeta, signs): """Serendipity corner shape function for a given (xi_i, eta_i, zeta_i).""" xi_i, eta_i, zeta_i = signs return ( 0.125 * (1 + xi_i * xi) * (1 + eta_i * eta) * (1 + zeta_i * zeta) * (xi_i * xi + eta_i * eta + zeta_i * zeta - 2) ) def n_mid_edge(xi, eta, zeta, signs): """Serendipity mid-edge shape function: the zero-axis carries (1 - .^2).""" xi_i, eta_i, zeta_i = signs if xi_i == 0: return 0.25 * (1 - xi**2) * (1 + eta_i * eta) * (1 + zeta_i * zeta) if eta_i == 0: return 0.25 * (1 - eta**2) * (1 + xi_i * xi) * (1 + zeta_i * zeta) return 0.25 * (1 - zeta**2) * (1 + xi_i * xi) * (1 + eta_i * eta) .. GENERATED FROM PYTHON SOURCE LINES 133-134 Sample the bottom face zeta = -1 on a (xi, eta) grid .. GENERATED FROM PYTHON SOURCE LINES 134-139 .. code-block:: Python n = 51 xi_grid, eta_grid = np.meshgrid(np.linspace(-1, 1, n), np.linspace(-1, 1, n), indexing="ij") zeta_slice = -1.0 .. GENERATED FROM PYTHON SOURCE LINES 140-147 Render — 4×2 grid of contour plots Top row: the four corner shape functions whose nodes lie on the bottom face (corners 0..3, all with zeta_i = -1). Bottom row: the four mid-edge functions whose nodes also lie on the bottom face (mid-edges 8..11). The other twelve basis functions vanish identically on this slice. .. GENERATED FROM PYTHON SOURCE LINES 147-181 .. code-block:: Python bottom_face_corners = [0, 1, 2, 3] bottom_face_mid_edges = [8, 9, 10, 11] fig, axes = plt.subplots(2, 4, figsize=(15, 7), constrained_layout=True) fig.suptitle( "HEX20 shape functions on the $\\zeta = -1$ face slice", fontsize=13, ) for col, ci in enumerate(bottom_face_corners): ax = axes[0, col] N = n_corner(xi_grid, eta_grid, zeta_slice, corner_signs[ci]) img = ax.contourf(xi_grid, eta_grid, N, levels=21, cmap="plasma") ax.set_aspect("equal") sx, sy, _ = corner_signs[ci].astype(int) ax.set_title(f"$N^c_{{{ci}}}$ $(\\xi_i,\\eta_i,\\zeta_i)=({sx:+d},{sy:+d},-1)$", fontsize=9) ax.set_xlabel(r"$\xi$") ax.set_ylabel(r"$\eta$") fig.colorbar(img, ax=ax, shrink=0.8) for col, mj in enumerate(bottom_face_mid_edges): ax = axes[1, col] N = n_mid_edge(xi_grid, eta_grid, zeta_slice, mid_edge_signs[mj - 8]) img = ax.contourf(xi_grid, eta_grid, N, levels=21, cmap="plasma") ax.set_aspect("equal") sx, sy, _ = mid_edge_signs[mj - 8].astype(int) ax.set_title(f"$N^m_{{{mj}}}$ $(\\xi_j,\\eta_j,\\zeta_j)=({sx:+d},{sy:+d},-1)$", fontsize=9) ax.set_xlabel(r"$\xi$") ax.set_ylabel(r"$\eta$") fig.colorbar(img, ax=ax, shrink=0.8) fig.show() .. image-sg:: /gallery/elements/solid186/images/sphx_glr_example_hex20_shape_function_slices_001.png :alt: HEX20 shape functions on the $\zeta = -1$ face slice, $N^c_{0}$ $(\xi_i,\eta_i,\zeta_i)=(-1,-1,-1)$, $N^c_{1}$ $(\xi_i,\eta_i,\zeta_i)=(+1,-1,-1)$, $N^c_{2}$ $(\xi_i,\eta_i,\zeta_i)=(+1,+1,-1)$, $N^c_{3}$ $(\xi_i,\eta_i,\zeta_i)=(-1,+1,-1)$, $N^m_{8}$ $(\xi_j,\eta_j,\zeta_j)=(+0,-1,-1)$, $N^m_{9}$ $(\xi_j,\eta_j,\zeta_j)=(+1,+0,-1)$, $N^m_{10}$ $(\xi_j,\eta_j,\zeta_j)=(+0,+1,-1)$, $N^m_{11}$ $(\xi_j,\eta_j,\zeta_j)=(-1,+0,-1)$ :srcset: /gallery/elements/solid186/images/sphx_glr_example_hex20_shape_function_slices_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 182-186 Verify partition of unity at every reference-element node Build the 20 nodal coordinates and check :math:`\sum_i N^c_i + \sum_j N^m_j \equiv 1` to machine precision. .. GENERATED FROM PYTHON SOURCE LINES 186-197 .. code-block:: Python all_nodes = np.vstack([corner_signs, mid_edge_signs]) # (20, 3) total = np.zeros(20) for k in range(20): xi_k, eta_k, zeta_k = all_nodes[k] total[k] = sum(n_corner(xi_k, eta_k, zeta_k, corner_signs[i]) for i in range(8)) + sum( n_mid_edge(xi_k, eta_k, zeta_k, mid_edge_signs[j]) for j in range(12) ) np.testing.assert_allclose(total, 1.0, atol=1e-14) print("OK — HEX20 partition-of-unity holds at every reference-element node.") .. rst-class:: sphx-glr-script-out .. code-block:: none OK — HEX20 partition-of-unity holds at every reference-element node. .. GENERATED FROM PYTHON SOURCE LINES 198-202 Verify Kronecker-delta interpolation The (20, 20) basis matrix evaluated at every reference-element node should be the identity. .. GENERATED FROM PYTHON SOURCE LINES 202-212 .. code-block:: Python basis = np.zeros((20, 20)) for k, (xi_k, eta_k, zeta_k) in enumerate(all_nodes): for i in range(8): basis[k, i] = n_corner(xi_k, eta_k, zeta_k, corner_signs[i]) for j in range(12): basis[k, 8 + j] = n_mid_edge(xi_k, eta_k, zeta_k, mid_edge_signs[j]) np.testing.assert_allclose(basis, np.eye(20), atol=1e-14) print("OK — HEX20 Kronecker-delta interpolation verified at every node.") .. rst-class:: sphx-glr-script-out .. code-block:: none OK — HEX20 Kronecker-delta interpolation verified at every node. .. GENERATED FROM PYTHON SOURCE LINES 213-219 Verify the bottom-face vanishing pattern On the zeta = -1 slice, all four top-face corners (4..7), all four top-face mid-edges (12..15), and all four vertical mid-edges (16..19) must be identically zero. The four bottom-face corners (0..3) and four bottom-face mid-edges (8..11) carry every value. .. GENERATED FROM PYTHON SOURCE LINES 219-228 .. code-block:: Python xi_test, eta_test = np.meshgrid(np.linspace(-1, 1, 5), np.linspace(-1, 1, 5), indexing="ij") for i in itertools.chain(range(4, 8)): Ni = n_corner(xi_test, eta_test, -1.0, corner_signs[i]) np.testing.assert_allclose(Ni, 0.0, atol=1e-14) for j in itertools.chain(range(4, 12)): Nj = n_mid_edge(xi_test, eta_test, -1.0, mid_edge_signs[j]) np.testing.assert_allclose(Nj, 0.0, atol=1e-14) print("OK — twelve off-face HEX20 basis functions vanish identically on zeta = -1.") .. rst-class:: sphx-glr-script-out .. code-block:: none OK — twelve off-face HEX20 basis functions vanish identically on zeta = -1. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.100 seconds) .. _sphx_glr_download_gallery_elements_solid186_example_hex20_shape_function_slices.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_hex20_shape_function_slices.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_hex20_shape_function_slices.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_hex20_shape_function_slices.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_