HEX20 shape-function slices on the bottom face#

Renders the eight HEX20 serendipity shape functions that do not vanish on the $\zeta = -1$ face of the reference cube $\hat\Omega = [-1, +1]^3$: four corner functions $N^c_i$ whose nodes lie on that face, plus four mid-edge functions $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 $\zeta_j = 0$) all vanish at $\zeta = -1$ and are omitted from the figure.

Serendipity basis#

In reference coordinates $(\xi, \eta, \zeta) \in [-1, +1]^3$, with corner-node signs $(\xi_i, \eta_i, \zeta_i) \in \{-1, +1\}^3$,

\[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 $\xi_j = 0$,

\[N^m_j(\xi, \eta, \zeta) = \tfrac{1}{4}(1 - \xi^{2}) (1 + \eta_j\eta)(1 + \zeta_j\zeta),\]

with analogous forms for $\eta_j = 0$ and $\zeta_j = 0$ mid-edges (the missing-variable axis carries the $(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: femorph_solver.elements.hex20.Hex20.

from __future__ import annotations

import itertools

import matplotlib.pyplot as plt
import numpy as np

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.

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)

Sample the bottom face zeta = -1 on a (xi, eta) grid

n = 51
xi_grid, eta_grid = np.meshgrid(np.linspace(-1, 1, n), np.linspace(-1, 1, n), indexing="ij")
zeta_slice = -1.0

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.

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()

$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)$$

Verify partition of unity at every reference-element node

Build the 20 nodal coordinates and check $\sum_i N^c_i + \sum_j N^m_j \equiv 1$ to machine precision.

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.")

OK — HEX20 partition-of-unity holds at every reference-element node.

Verify Kronecker-delta interpolation

The (20, 20) basis matrix evaluated at every reference-element node should be the identity.

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.")

OK — HEX20 Kronecker-delta interpolation verified at every node.

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.

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.")

OK — twelve off-face HEX20 basis functions vanish identically on zeta = -1.

Total running time of the script: (0 minutes 1.100 seconds)

Gallery generated by Sphinx-Gallery