HEX20 reference geometry — corner + mid-edge nodes, Gauss points

Renders the 20-node serendipity hexahedron’s reference element \(\hat\Omega = [-1, 1]^3\):

• 8 corner nodes at \((\xi, \eta, \zeta) \in \{-1, +1\}^3\).

• 12 mid-edge nodes at the centre of every cube edge.

• 14-point Irons (1971) quadrature points used for the consistent-mass integral on this kernel.

• 2×2×2 reduced-integration Gauss points \((\xi, \eta, \zeta) \in \{\pm 1/\sqrt{3}\}^3\) used for the stiffness integral by default in femorph-solver (the femorph-solver default for HEX20).

Shape functions

The serendipity HEX20 basis (Hughes 2000 §3.7; Zienkiewicz & Taylor 2013 §6.6) splits the 20 nodal functions into corner and mid-edge groups:

\[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),\]

\[N^m_j(\xi, \eta, \zeta) = \tfrac{1}{4} (1 - \xi^{2})(1 + \eta_j\eta)(1 + \zeta_j\zeta) \quad \text{(for the 4 mid-edges with } \xi_j = 0\text{)},\]

with analogous expressions for the four mid-edges along each of \(\eta_j = 0\) and \(\zeta_j = 0\).

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 — original serendipity basis.

• Hughes, T. J. R. (2000) The Finite Element Method, Dover, §3.7.

• Irons, B. M. (1971) “Quadrature rules for brick based finite elements,” International Journal for Numerical Methods in Engineering 3 (2), 293–294 — 14-point degree-5 rule.

• 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 numpy as np
import pyvista as pv

Node layout

8 corners + 12 mid-edges = 20 nodes. We label corner indices 0..7 and mid-edge indices 8..19 to match the standard isoparametric ordering used by VTK_QUADRATIC_HEXAHEDRON.

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,
)
mid_edges = np.array(
    [
        # bottom (z = -1) ring: 8, 9, 10, 11
        [0, -1, -1],
        [+1, 0, -1],
        [0, +1, -1],
        [-1, 0, -1],
        # top (z = +1) ring: 12, 13, 14, 15
        [0, -1, +1],
        [+1, 0, +1],
        [0, +1, +1],
        [-1, 0, +1],
        # vertical mid-edges (z = 0 ring): 16, 17, 18, 19
        [-1, -1, 0],
        [+1, -1, 0],
        [+1, +1, 0],
        [-1, +1, 0],
    ],
    dtype=float,
)
all_nodes = np.vstack([corners, mid_edges])

cells = np.hstack([[20], np.arange(20, dtype=np.int64)])
cell_types = np.array([pv.CellType.QUADRATIC_HEXAHEDRON], dtype=np.uint8)
ref_cube = pv.UnstructuredGrid(cells, cell_types, all_nodes)

Quadrature points — the two rules HEX20 uses

# 2x2x2 Gauss-Legendre — stiffness (reduced).
g = 1.0 / np.sqrt(3.0)
gauss_8 = np.array(
    [(a, b, c) for a in (-g, +g) for b in (-g, +g) for c in (-g, +g)],
    dtype=float,
)

# Irons 14-point — consistent mass.  Six face-centre points at
# (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) and 8 corner-octant points at
# ±sqrt(19/30)·(1, 1, 1) with sign permutations.
b_ = float(np.sqrt(19.0 / 30.0))
faces = np.array(
    [[+1, 0, 0], [-1, 0, 0], [0, +1, 0], [0, -1, 0], [0, 0, +1], [0, 0, -1]],
    dtype=float,
)
octants = np.array(
    [[a * b_, b * b_, c * b_] for a in (-1, +1) for b in (-1, +1) for c in (-1, +1)],
    dtype=float,
)
irons_14 = np.vstack([faces, octants])

Render

plotter = pv.Plotter(off_screen=True, window_size=(640, 480))
plotter.add_mesh(
    ref_cube,
    style="wireframe",
    color="black",
    line_width=2,
    label="reference HEX20",
)
plotter.add_points(
    corners,
    render_points_as_spheres=True,
    point_size=18,
    color="black",
    label="corner nodes (8)",
)
plotter.add_points(
    mid_edges,
    render_points_as_spheres=True,
    point_size=12,
    color="#1f77b4",
    label="mid-edge nodes (12)",
)
plotter.add_points(
    gauss_8,
    render_points_as_spheres=True,
    point_size=14,
    color="#d62728",
    label="2×2×2 Gauss (stiffness)",
)
plotter.add_points(
    irons_14,
    render_points_as_spheres=True,
    point_size=10,
    color="#2ca02c",
    label="Irons 14-pt (mass)",
)
plotter.add_axes(line_width=4, color="black")
plotter.view_isometric()
plotter.camera.zoom(1.05)
plotter.add_legend(face=None, size=(0.24, 0.16), bcolor="white")
plotter.show()

Sanity: every Irons-14 point is inside \(\hat\Omega\)

assert (np.abs(irons_14) <= 1.0 + 1e-12).all(), "Irons-14 point left the cube"
print(f"Irons 14-pt rule: {len(irons_14)} points, all inside [-1, 1]^3.")
print("  6 face-centre points at ±1 along one axis")
print("  8 corner-octant points at ±sqrt(19/30)·(1,1,1) sign-perm")

Irons 14-pt rule: 14 points, all inside [-1, 1]^3.
  6 face-centre points at ±1 along one axis
  8 corner-octant points at ±sqrt(19/30)·(1,1,1) sign-perm