HEX8 shape-function surfaces on the reference cube

Renders all eight trilinear Lagrange shape functions \(N_i(\xi, \eta, \zeta)\) for the 8-node hexahedron on the reference cube \(\hat\Omega = [-1, +1]^3\). The basis functions are

\[N_i(\xi, \eta, \zeta) = \tfrac{1}{8} (1 + \xi_i \xi)(1 + \eta_i \eta)(1 + \zeta_i \zeta), \qquad i = 1, \ldots, 8,\]

with \((\xi_i, \eta_i, \zeta_i) \in \{-1, +1\}^3\) the nodal coordinates. Each \(N_i\) is a tensor product of 1-D linear “hat” functions, equal to 1 at corner \(i\) and 0 at every other corner — the Kronecker-delta interpolation property.

The plot below shows \(N_i\) evaluated on the \(\zeta = -1\) cube face (\(i = 1, \ldots, 4\)) and the \(\zeta = +1\) face (\(i = 5, \ldots, 8\)).

References

• Hughes, T. J. R. (2000) The Finite Element Method — Linear Static and Dynamic Finite Element Analysis, Dover, §3.1 — trilinear shape functions.

• Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, §6.2.

Implementation: femorph_solver.elements.hex8.Hex8.

from __future__ import annotations

import matplotlib.pyplot as plt
import numpy as np

Corner natural-coordinate map (matches VTK_HEXAHEDRON ordering)

corners = 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,
)


def n_i(xi, eta, zeta, i):
    """Trilinear Lagrange shape function for corner i."""
    xi_i, eta_i, zeta_i = corners[i]
    return 0.125 * (1 + xi_i * xi) * (1 + eta_i * eta) * (1 + zeta_i * zeta)

Sample on a 41×41 grid on each ζ = ±1 face

n = 41
xi = eta = np.linspace(-1.0, 1.0, n)
XI, ETA = np.meshgrid(xi, eta, indexing="ij")

4×2 surface plot — one column per ζ-face, 4 rows for the 4 corners

fig, axes = plt.subplots(
    4,
    2,
    figsize=(10.0, 14.0),
    subplot_kw={"projection": "3d"},
    constrained_layout=True,
)
fig.suptitle("HEX8 shape functions $N_i(\\xi, \\eta, \\zeta)$ on cube faces", fontsize=13)

face_zeta = (-1.0, +1.0)
node_indices_per_face = ((0, 1, 2, 3), (4, 5, 6, 7))
for col, (zeta, idxs) in enumerate(zip(face_zeta, node_indices_per_face, strict=True)):
    for row, idx in enumerate(idxs):
        ax = axes[row, col]
        N = n_i(XI, ETA, zeta, idx)
        ax.plot_surface(XI, ETA, N, cmap="plasma", linewidth=0, antialiased=True)
        x_i, y_i, _ = corners[idx]
        ax.set_title(
            f"$N_{{{idx + 1}}}$ at corner $({int(x_i):+d}, {int(y_i):+d}, {int(zeta):+d})$",
            fontsize=10,
        )
        ax.set_xlabel(r"$\xi$")
        ax.set_ylabel(r"$\eta$")
        ax.set_zlabel(r"$N$")
        ax.set_xticks([-1, 0, 1])
        ax.set_yticks([-1, 0, 1])
        ax.set_zticks([0, 0.5, 1])
        ax.view_init(elev=22, azim=-55)

fig.show()

Verify partition of unity at every interior sample point

\(\sum_{i=1}^8 N_i(\xi, \eta, \zeta) = 1\) for any \((\xi, \eta, \zeta) \in \hat\Omega\). Sample on a 21×21×21 lattice and assert.

n3 = 21
xi3 = eta3 = zeta3 = np.linspace(-1.0, 1.0, n3)
XI3, ETA3, ZETA3 = np.meshgrid(xi3, eta3, zeta3, indexing="ij")
total = sum(n_i(XI3, ETA3, ZETA3, idx) for idx in range(8))
max_dev = float(np.max(np.abs(total - 1.0)))
assert max_dev < 1.0e-14, f"partition-of-unity violated: max |Σ N_i - 1| = {max_dev:.3e}"
print(f"max |Σ N_i - 1| over a 21^3 grid = {max_dev:.3e}")
print("OK — eight trilinear HEX8 shape functions form a partition of unity.")

max |Σ N_i - 1| over a 21^3 grid = 4.441e-16
OK — eight trilinear HEX8 shape functions form a partition of unity.

Verify Kronecker-delta interpolation at every corner

\(N_i(\boldsymbol{\xi}_j) = \delta_{ij}\): each shape function equals 1 at its own corner and 0 at all other corners — required so the nodal-displacement DOFs are literally the displacement values.

print()
for i in range(8):
    row = [n_i(*corners[j], i) for j in range(8)]
    expected = [1.0 if j == i else 0.0 for j in range(8)]
    np.testing.assert_allclose(row, expected, atol=1.0e-15)
print("OK — Kronecker-delta interpolation verified at every corner.")

OK — Kronecker-delta interpolation verified at every corner.