"""Hex20 — 20-node quadratic hexahedral solid (serendipity).
Node ordering matches VTK_QUADRATIC_HEXAHEDRON: the 8 corner nodes
come first in Hex8 order, then the 4 bottom-face mid-edge nodes,
the 4 top-face mid-edge nodes, and finally the 4 vertical mid-edge
nodes connecting the bottom and top faces.
Shape functions on the reference cube ``ξ, η, ζ ∈ [-1, 1]``:
* **Corner** node with ``(ξᵢ, ηᵢ, ζᵢ) ∈ {±1}³``::
Nᵢ = (1/8)(1 + ξᵢ ξ)(1 + ηᵢ η)(1 + ζᵢ ζ)
· (ξᵢ ξ + ηᵢ η + ζᵢ ζ - 2)
* **Mid-edge** node with one natural coordinate zero (e.g. ``ξᵢ = 0``)::
Nᵢ = (1/4)(1 - ξ²)(1 + ηᵢ η)(1 + ζᵢ ζ)
Integration rule
----------------
``material["_HEX20_INTEGRATION"]`` selects the stiffness integration
rule:
* ``"reduced"`` (default) — 2×2×2 Gauss-Legendre (8 points). Single-
element Ke has 6 extra zero-energy hourglass modes; on a connected
mesh those modes almost always die under the assembly because a hex
whose midside is shared with a neighbour has no free hourglass
direction to excite. Make sure your mesh has no dangling midside
nodes before using the default; if in doubt switch to ``"full"``.
* ``"full"`` — 3×3×3 Gauss-Legendre (27 points). Ke is fully rank
(54 non-zero eigenvalues + 6 rigid-body zeros) on a single hex.
Mass integration uses the 14-point Irons rule by default; pass
``material["_HEX20_MASS"] = "consistent"`` for the textbook 3×3×3
consistent rule. The wedge / pyramid degenerate forms also use the
Irons 14-point rule (see ``_wedge15_pyr13.py``); the rule is
documented in the References section below.
Real constants
--------------
None required — geometry comes from node coordinates.
Material
--------
EX, PRXY : elastic moduli (both mandatory)
DENS : density (required for M_e)
References
----------
* Shape functions — 20-node serendipity hexahedron (corner nodes
weighted by ``(ξᵢξ + ηᵢη + ζᵢζ − 2)`` correction, mid-edge nodes
on zero-coordinate axes): Zienkiewicz, O.C. and Taylor, R.L.
*The Finite Element Method: Its Basis and Fundamentals*, 7th
ed., Butterworth-Heinemann, 2013, §8.7.3. Original derivation
in Ergatoudis, I., Irons, B.M. and Zienkiewicz, O.C.,
"Curved, isoparametric, 'quadrilateral' elements for finite
element analysis," *Int. J. Solids Struct.* 4 (1968), pp.
31–42, extended to 20-node hex in Zienkiewicz, O.C., *The
Finite Element Method in Engineering Science*, 2nd ed.,
McGraw-Hill, 1971.
* **2×2×2 reduced Gauss integration** (``"reduced"``, default):
Bedrosian, G., "Shape functions and integration formulas for
three-dimensional finite element analysis," *IJNME* 35 (1992),
pp. 95–108 — the canonical reference for the uniform-reduced
rule on the 20-node hex and its degenerate wedge / pyramid
forms. Hourglass-mode analysis of the reduced rule on a single
hex: Belytschko, T., Liu, W.K., Moran, B. and Elkhodary, K.,
*Nonlinear Finite Elements for Continua and Structures*, 2nd
ed., Wiley, 2014, §8.6.
* 3×3×3 full Gauss-Legendre rule (``"full"``, 27 points):
Zienkiewicz & Taylor §5.10 Table 5.3.
* **14-point Irons integration rule** (used for the consistent
mass and in the degenerate wedge / pyramid forms): Irons,
B.M., "Quadrature rules for brick based finite elements,"
*IJNME* 3 (1971), pp. 293–294. Originally derived for the
20-node hex with degenerate wedge / pyramid as special cases,
integrating polynomials through degree 5 exactly.
* Voigt elastic matrix (isotropic linear elasticity, engineering
shear convention): Cook, Malkus, Plesha, Witt *Concepts and
Applications of FEA*, 4th ed., Wiley, 2002, §3.7.
"""
from __future__ import annotations
import numpy as np
from femorph_solver.elements._base import ElementBase
from femorph_solver.elements._registry import register
from femorph_solver.elements.hex8 import _elastic_matrix
try:
from femorph_solver import _elements as _ext
except ImportError: # pragma: no cover
_ext = None
# Natural coordinates of the 20 nodes, in VTK_QUADRATIC_HEXAHEDRON
# order: 8 corners, then bottom mid-edges (I-J, J-K, K-L, L-I), then top
# mid-edges (M-N, N-O, O-P, P-M), then vertical mid-edges (I-M, J-N, K-O,
# L-P).
_XI_NODES = np.array(
[
# corners
[-1, -1, -1],
[+1, -1, -1],
[+1, +1, -1],
[-1, +1, -1],
[-1, -1, +1],
[+1, -1, +1],
[+1, +1, +1],
[-1, +1, +1],
# bottom-face mid-edges
[0, -1, -1],
[+1, 0, -1],
[0, +1, -1],
[-1, 0, -1],
# top-face mid-edges
[0, -1, +1],
[+1, 0, +1],
[0, +1, +1],
[-1, 0, +1],
# vertical mid-edges
[-1, -1, 0],
[+1, -1, 0],
[+1, +1, 0],
[-1, +1, 0],
],
dtype=np.float64,
)
N_NODES = 20
N_DOF = 3 * N_NODES # 60
def _gauss_3x3x3() -> tuple[np.ndarray, np.ndarray]:
"""3×3×3 Gauss-Legendre rule on [-1, 1]³ (27 points)."""
pts_1d = np.array([-np.sqrt(3.0 / 5.0), 0.0, +np.sqrt(3.0 / 5.0)], dtype=np.float64)
w_1d = np.array([5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0], dtype=np.float64)
pts = np.empty((27, 3), dtype=np.float64)
wts = np.empty(27, dtype=np.float64)
k = 0
for i in range(3):
for j in range(3):
for kk in range(3):
pts[k] = (pts_1d[i], pts_1d[j], pts_1d[kk])
wts[k] = w_1d[i] * w_1d[j] * w_1d[kk]
k += 1
return pts, wts
_GAUSS_POINTS, _GAUSS_WEIGHTS = _gauss_3x3x3()
def _gauss_2x2x2() -> tuple[np.ndarray, np.ndarray]:
"""2×2×2 Gauss-Legendre rule on [-1, 1]³ (8 points, uniform reduced)."""
g = 1.0 / np.sqrt(3.0)
pts = np.array(
[
[-g, -g, -g],
[+g, -g, -g],
[+g, +g, -g],
[-g, +g, -g],
[-g, -g, +g],
[+g, -g, +g],
[+g, +g, +g],
[-g, +g, +g],
],
dtype=np.float64,
)
wts = np.ones(8, dtype=np.float64)
return pts, wts
_GAUSS_POINTS_REDUCED, _GAUSS_WEIGHTS_REDUCED = _gauss_2x2x2()
def _irons_14pt() -> tuple[np.ndarray, np.ndarray]:
"""14-point Irons rule on the reference cube — Hex20 default mass rule.
This is **not** the 2×2×2 rule used for URI stiffness and **not**
the 3×3×3 consistent rule a textbook would prescribe; femorph-solver
uses Irons' (1971) lower-order rule by default because it drops the
full midside coupling and gives a denser-but-cheaper mass operator
that matches the natural symmetry of the serendipity hex. Pass
``material["_HEX20_MASS"] = "consistent"`` to opt into the textbook
3×3×3 rule.
Layout:
* 6 face-center points at ``(±a, 0, 0) / (0, ±a, 0) / (0, 0, ±a)``
with ``a = 0.7958224257542215`` and weight ``0.886426592797784``.
* 8 "corner" points at ``(±b, ±b, ±b)`` with
``b = 0.7587869106393281`` and weight ``0.3351800554016621``.
Reference: Irons, B. M. (1971) "Quadrature rules for brick based
finite elements", Int. J. Numer. Meth. Engng. 3 (2), 293–294.
"""
a = 0.7958224257542215
b = 0.7587869106393281
w_face = 0.8864265927977839
w_corner = 0.3351800554016621
face_pts = []
for coord_idx in range(3):
for sign in (-1.0, +1.0):
p = [0.0, 0.0, 0.0]
p[coord_idx] = sign * a
face_pts.append(p)
corner_pts = [
[sx * b, sy * b, sz * b] for sx in (-1.0, 1.0) for sy in (-1.0, 1.0) for sz in (-1.0, 1.0)
]
pts = np.array(face_pts + corner_pts, dtype=np.float64)
wts = np.concatenate(
[np.full(6, w_face, dtype=np.float64), np.full(8, w_corner, dtype=np.float64)]
)
return pts, wts
_IRONS_14_POINTS, _IRONS_14_WEIGHTS = _irons_14pt()
def _shape_and_derivs(xi: float, eta: float, zeta: float) -> tuple[np.ndarray, np.ndarray]:
"""Return ``(N (20,), dN/dξ (20, 3))`` for the 20-node serendipity hex."""
N = np.empty(N_NODES, dtype=np.float64)
dN = np.empty((N_NODES, 3), dtype=np.float64)
for i in range(N_NODES):
xi_i, eta_i, zeta_i = _XI_NODES[i]
abs_xi, abs_eta, abs_zeta = abs(xi_i), abs(eta_i), abs(zeta_i)
one_xi = 1.0 + xi_i * xi
one_eta = 1.0 + eta_i * eta
one_zeta = 1.0 + zeta_i * zeta
if abs_xi == 1 and abs_eta == 1 and abs_zeta == 1:
# Corner node.
combo = xi_i * xi + eta_i * eta + zeta_i * zeta - 2.0
N[i] = 0.125 * one_xi * one_eta * one_zeta * combo
# d/dξ: product rule across the four factors, two of which
# depend on ξ (one_xi and combo).
dN[i, 0] = 0.125 * one_eta * one_zeta * (xi_i * combo + one_xi * xi_i)
dN[i, 1] = 0.125 * one_xi * one_zeta * (eta_i * combo + one_eta * eta_i)
dN[i, 2] = 0.125 * one_xi * one_eta * (zeta_i * combo + one_zeta * zeta_i)
elif abs_xi == 0:
# Mid-edge on a ξ-directed edge.
q = 1.0 - xi * xi
N[i] = 0.25 * q * one_eta * one_zeta
dN[i, 0] = 0.25 * (-2.0 * xi) * one_eta * one_zeta
dN[i, 1] = 0.25 * q * eta_i * one_zeta
dN[i, 2] = 0.25 * q * one_eta * zeta_i
elif abs_eta == 0:
q = 1.0 - eta * eta
N[i] = 0.25 * q * one_xi * one_zeta
dN[i, 0] = 0.25 * q * xi_i * one_zeta
dN[i, 1] = 0.25 * (-2.0 * eta) * one_xi * one_zeta
dN[i, 2] = 0.25 * q * one_xi * zeta_i
else: # abs_zeta == 0
q = 1.0 - zeta * zeta
N[i] = 0.25 * q * one_xi * one_eta
dN[i, 0] = 0.25 * q * xi_i * one_eta
dN[i, 1] = 0.25 * q * one_xi * eta_i
dN[i, 2] = 0.25 * (-2.0 * zeta) * one_xi * one_eta
return N, dN
def _b_matrix(dN_dx: np.ndarray) -> np.ndarray:
"""Return the 6×60 strain-displacement matrix from ``dN/dx (20, 3)``."""
B = np.zeros((6, N_DOF), dtype=np.float64)
for i in range(N_NODES):
bx, by, bz = dN_dx[i]
col = 3 * i
B[0, col + 0] = bx
B[1, col + 1] = by
B[2, col + 2] = bz
B[3, col + 0] = by
B[3, col + 1] = bx
B[4, col + 1] = bz
B[4, col + 2] = by
B[5, col + 0] = bz
B[5, col + 2] = bx
return B
def _shape_matrix(N: np.ndarray) -> np.ndarray:
"""Return the 3×60 shape-function matrix ``Φ`` with three rows per DOF."""
Phi = np.zeros((3, N_DOF), dtype=np.float64)
for i in range(N_NODES):
Phi[0, 3 * i + 0] = N[i]
Phi[1, 3 * i + 1] = N[i]
Phi[2, 3 * i + 2] = N[i]
return Phi
# Shape-function derivatives at each Gauss point depend only on ``ξ``,
# not on the element coordinates — cache them at module load so the hot
# assembly path doesn't recompute them per element.
_DN_DXI_FULL = np.stack([_shape_and_derivs(*gp)[1] for gp in _GAUSS_POINTS], axis=0) # (27, 20, 3)
_DN_DXI_REDUCED = np.stack(
[_shape_and_derivs(*gp)[1] for gp in _GAUSS_POINTS_REDUCED], axis=0
) # (8, 20, 3)
def _ke_python(
coords: np.ndarray,
E: float,
nu: float,
gauss_points: np.ndarray,
gauss_weights: np.ndarray,
) -> np.ndarray:
"""Assemble Ke from an arbitrary Gauss rule (Python reference path)."""
C = _elastic_matrix(E, nu)
# Fast path: cached dN_dxi for the two standard rules.
if gauss_points is _GAUSS_POINTS:
dN_dxi_cache = _DN_DXI_FULL
elif gauss_points is _GAUSS_POINTS_REDUCED:
dN_dxi_cache = _DN_DXI_REDUCED
else:
dN_dxi_cache = None
K = np.zeros((N_DOF, N_DOF), dtype=np.float64)
for g, (gp, w) in enumerate(zip(gauss_points, gauss_weights)):
if dN_dxi_cache is not None:
dN_dxi = dN_dxi_cache[g]
else:
_, dN_dxi = _shape_and_derivs(*gp)
J = dN_dxi.T @ coords
det_J = float(np.linalg.det(J))
if det_J <= 0:
raise ValueError(f"HEX20: non-positive Jacobian {det_J:.3e} — check node ordering")
dN_dx = np.linalg.solve(J, dN_dxi.T).T
B = _b_matrix(dN_dx)
K += B.T @ C @ B * det_J * w
return K
def _ke_bbar(coords: np.ndarray, E: float, nu: float) -> np.ndarray:
"""Hughes B-bar stiffness for the 20-node hex (``"full"`` integration).
Same closed-form correction used in Hex8's B-bar path:
K_bbar = K_plain + K_B · ((1/V_e) S Sᵀ − H)
with K_B = λ + (2/3)G, S = Σ_gp b(gp)|J|w, H = Σ_gp b bᵀ|J|w, and
b(gp) = [∂N_i/∂x_d] is the 60-vector volumetric operator at gp.
The correction adds one O(N_DOF²) pass over the 3×3×3 Gauss points on
top of the plain-Gauss accumulation.
"""
lam = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
G = E / (2.0 * (1.0 + nu))
K_B = lam + (2.0 / 3.0) * G
C = _elastic_matrix(E, nu)
K = np.zeros((N_DOF, N_DOF), dtype=np.float64)
S = np.zeros(N_DOF, dtype=np.float64)
H = np.zeros((N_DOF, N_DOF), dtype=np.float64)
V_e = 0.0
for gp, w in zip(_GAUSS_POINTS, _GAUSS_WEIGHTS):
_, dN_dxi = _shape_and_derivs(*gp)
J = dN_dxi.T @ coords
det_J = float(np.linalg.det(J))
if det_J <= 0:
raise ValueError(
f"Hex20 (B-bar): non-positive Jacobian {det_J:.3e} — check node ordering"
)
dN_dx = np.linalg.solve(J, dN_dxi.T).T
B = _b_matrix(dN_dx)
b = np.empty(N_DOF, dtype=np.float64)
b[0::3] = dN_dx[:, 0]
b[1::3] = dN_dx[:, 1]
b[2::3] = dN_dx[:, 2]
dv = det_J * w
V_e += dv
S += b * dv
H += np.outer(b, b) * dv
K += B.T @ C @ B * dv
K += K_B * (np.outer(S, S) / V_e - H)
return K
[docs]
@register
class Hex20(ElementBase):
name = "HEX20"
n_nodes = N_NODES
n_dof_per_node = 3
vtk_cell_type = 25 # VTK_QUADRATIC_HEXAHEDRON
[docs]
@staticmethod
def ke(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
) -> np.ndarray:
coords = np.ascontiguousarray(coords, dtype=np.float64)
if coords.shape != (N_NODES, 3):
raise ValueError(f"Hex20 expects ({N_NODES}, 3) coords, got {coords.shape}")
E = float(material["EX"])
nu = float(material["PRXY"])
formulation = material.get("_HEX20_INTEGRATION", "reduced")
if formulation == "reduced":
# 2×2×2 Gauss — the femorph-solver default. Single-element Ke
# is rank-deficient (6 hourglass modes) but that's harmless
# in any connected mesh with shared midsides.
if _ext is not None and hasattr(_ext, "solid186_ke_reduced"):
return _ext.solid186_ke_reduced(coords, E, nu)
return _ke_python(coords, E, nu, _GAUSS_POINTS_REDUCED, _GAUSS_WEIGHTS_REDUCED)
if formulation == "full":
# 3×3×3 Gauss — slower than reduced; opt-in for locking-prone
# meshes where the default's hourglass modes can survive.
if _ext is not None:
return _ext.solid186_ke(coords, E, nu)
return _ke_python(coords, E, nu, _GAUSS_POINTS, _GAUSS_WEIGHTS)
raise ValueError(
f"HEX20: unknown _HEX20_INTEGRATION {formulation!r}; expected 'full' or 'reduced'"
)
[docs]
@staticmethod
def me(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
lumped: bool = False,
) -> np.ndarray:
coords = np.ascontiguousarray(coords, dtype=np.float64)
if coords.shape != (N_NODES, 3):
raise ValueError(f"Hex20 expects ({N_NODES}, 3) coords, got {coords.shape}")
rho = float(material["DENS"])
# Hex20 mass uses the 14-point Irons rule by default — drops the
# full midside coupling of the textbook 3×3×3 consistent rule for
# a denser-but-cheaper operator that matches the natural symmetry
# of the serendipity hex. Callers can opt into the textbook
# consistent rule via ``material["_HEX20_MASS"] = "consistent"``.
mass_rule = material.get("_HEX20_MASS", "irons14")
if mass_rule == "irons14":
pts, wts = _IRONS_14_POINTS, _IRONS_14_WEIGHTS
elif mass_rule == "consistent":
pts, wts = _GAUSS_POINTS, _GAUSS_WEIGHTS
else:
raise ValueError(
f"HEX20: unknown _HEX20_MASS {mass_rule!r}; expected 'irons14' or 'consistent'"
)
# C++ fast paths: prefer the kernel that matches ``mass_rule``.
if _ext is not None:
if mass_rule == "irons14" and hasattr(_ext, "solid186_me_irons14"):
return _ext.solid186_me_irons14(coords, rho, lumped)
if mass_rule == "consistent":
return _ext.solid186_me(coords, rho, lumped)
M = np.zeros((N_DOF, N_DOF), dtype=np.float64)
for gp, w in zip(pts, wts):
N, dN_dxi = _shape_and_derivs(*gp)
J = dN_dxi.T @ coords
det_J = float(np.linalg.det(J))
Phi = _shape_matrix(N)
M += rho * (Phi.T @ Phi) * det_J * w
if lumped:
M = np.diag(M.sum(axis=1))
return M
@staticmethod
def ke_batch(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
) -> np.ndarray | None:
if _ext is None:
return None
formulation = material.get("_HEX20_INTEGRATION", "reduced")
coords = np.ascontiguousarray(coords, dtype=np.float64)
E = float(material["EX"])
nu = float(material["PRXY"])
if formulation == "reduced" and hasattr(_ext, "solid186_ke_reduced_batch"):
return _ext.solid186_ke_reduced_batch(coords, E, nu)
if formulation == "full":
return _ext.solid186_ke_batch(coords, E, nu)
return None
@staticmethod
def me_batch(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
lumped: bool = False,
) -> np.ndarray | None:
if _ext is None:
return None
mass_rule = material.get("_HEX20_MASS", "irons14")
coords = np.ascontiguousarray(coords, dtype=np.float64)
rho = float(material["DENS"])
if mass_rule == "irons14" and hasattr(_ext, "solid186_me_irons14_batch"):
return _ext.solid186_me_irons14_batch(coords, rho, lumped)
if mass_rule == "consistent":
return _ext.solid186_me_batch(coords, rho, lumped)
# Unknown rule or C++ kernel missing — fall through to the Python loop.
return None
[docs]
@staticmethod
def eel_batch(
coords: np.ndarray,
u_e: np.ndarray,
material: dict[str, float] | None = None,
) -> np.ndarray | None:
"""Per-element elastic strain at every element node.
``coords``: ``(n_elem, 20, 3)``; ``u_e``: ``(n_elem, 60)``.
Returns ``(n_elem, 20, 6)`` — engineering-shear Voigt strain.
``material`` is accepted for signature uniformity with plane
kernels but is unused.
"""
del material # unused for 3D solids
if _ext is None:
return None
coords = np.ascontiguousarray(coords, dtype=np.float64)
u_e = np.ascontiguousarray(u_e, dtype=np.float64)
return _ext.solid186_strain_batch(coords, u_e)