HEX20 — 20-node quadratic hexahedral solid#
Kinematics. 20-node serendipity hexahedron; three translational DOFs per node (\(u_x, u_y, u_z\)); 60 DOFs per element. Also dispatches degenerate WEDGE15 (15-node) and PYR13 (13-node) shapes automatically when a foreign deck delivers collapsed-corner connectivity.
Stiffness. \(K_e = \int_{\Omega_e} B^\top C B\,\mathrm{d}V\)
integrated on the reference cube with 2×2×2 Gauss-Legendre
("reduced", default) or 3×3×3 ("full").
Mass. Consistent \(M_e = \int \rho\, N^\top N\,\mathrm{d}V\);
14-point Irons rule by default, 3×3×3 Gauss available via
material["_HEX20_MASS"] = "consistent".
Shape functions (serendipity 20-node hex)#
The 20-node brick is the quadratic serendipity extension of the trilinear HEX8 [ZTZ2013] [Bathe2014]. Eight corner nodes sit at \((\xi_i, \eta_i, \zeta_i) \in \{\pm 1\}^3\); twelve mid-edge nodes sit where exactly one natural coordinate is zero. There are no mid-face or body-centre nodes — that’s the hallmark of the serendipity family [Cook2001].
Corner node \(i\) (all three \(|\xi_i^c| = 1\)):
\[N_i(\xi, \eta, \zeta) = \tfrac{1}{8}\,(1 + \xi_i \xi)\,(1 + \eta_i \eta)\, (1 + \zeta_i \zeta)\,\bigl(\xi_i \xi + \eta_i \eta + \zeta_i \zeta - 2\bigr).\]Mid-edge node on a \(\xi\)-directed edge (\(\xi_i^m = 0\)):
\[N_i(\xi, \eta, \zeta) = \tfrac{1}{4}\,(1 - \xi^2)\,(1 + \eta_i \eta)\, (1 + \zeta_i \zeta),\]with the symmetric forms for \(\eta\)- and \(\zeta\)-directed mid-edges.
Jacobian, \(B\) matrix, and the isotropic constitutive \(C\) are the same six-component Voigt objects as for HEX8 — only the shape-function and derivative entries change.
3×3×3 Gauss quadrature#
Full integration uses 27 points at natural coordinates \(\xi_g \in \{-\sqrt{3/5},\, 0,\, +\sqrt{3/5}\}\) (and similarly for \(\eta\), \(\zeta\)) with weights \(w_g\) from the tensor product of the one-dimensional 3-point rule (\(5/9,\ 8/9,\ 5/9\)). The rule is exact for polynomials up to degree 5 in each direction, which integrates \(B^\top C B\) on an undistorted hex without approximation.
Reduced ("reduced", default) integration drops to the 2×2×2 rule
(same points as HEX8). On a single element this
produces six extra zero-energy modes (hourglass modes); in a patch of
elements the modes usually don’t connect and the global \(K\)
stays stable. No explicit hourglass stabilisation is applied yet —
on meshes with dangling midside nodes, switch to "full".
Numpy walk-through#
The following listing reuses the HEX8 pattern with the 20-node shape functions. On a unit cube with generic steel it produces the expected 6 rigid-body + 54 elastic eigenvalues:
import numpy as np
# Natural coordinates of the 20 nodes, in VTK_QUADRATIC_HEXAHEDRON
# order: corners (bottom 4, top 4), bottom-face midsides, top-face
# midsides, vertical midsides.
xi_n = 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),
( 0,-1,-1), (+1, 0,-1), ( 0,+1,-1), (-1, 0,-1),
( 0,-1,+1), (+1, 0,+1), ( 0,+1,+1), (-1, 0,+1),
(-1,-1, 0), (+1,-1, 0), (+1,+1, 0), (-1,+1, 0),
], dtype=float)
# 3x3x3 Gauss-Legendre.
p1d = np.array([-np.sqrt(3/5), 0.0, +np.sqrt(3/5)])
w1d = np.array([5/9, 8/9, 5/9])
gp = np.array([(a, b, c) for a in p1d for b in p1d for c in p1d])
w = np.array([a * b * c for a in w1d for b in w1d for c in w1d])
def shape_and_derivs(xi, eta, zeta):
N = np.empty(20)
dN = np.empty((20, 3))
for i, (xi_i, eta_i, zeta_i) in enumerate(xi_n):
one_xi = 1.0 + xi_i * xi
one_eta = 1.0 + eta_i * eta
one_zeta = 1.0 + zeta_i * zeta
if abs(xi_i) + abs(eta_i) + abs(zeta_i) == 3: # corner
combo = xi_i * xi + eta_i * eta + zeta_i * zeta - 2.0
N[i] = 0.125 * one_xi * one_eta * one_zeta * 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 xi_i == 0: # mid-edge on ξ
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 eta_i == 0: # mid-edge on η
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: # mid-edge on ζ
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
# Assemble K_e the same way as the HEX8 listing: loop over Gauss
# points, compute J = dN^T X, detJ, dN/dx = J^{-1} dN/dξ, build the
# 6x60 B, accumulate K += B^T C B detJ w.
The element passes the patch test to machine precision: a linear displacement field ramps exactly through every node of any distorted hex.
Degenerate shapes — wedge and pyramid#
Several foreign-deck formats represent quadratic wedge / pyramid elements as 20-slot hex connectivity with repeated node IDs at collapsed corners:
Wedge (15 nodes) — corner collapse
K = LandO = Pon adjacent vertical edges.Pyramid (13 nodes) — corner collapse
M = N = O = Pon all four top-face nodes (the apex).
Foreign-deck readers parse the degenerate connectivity and emit
VTK_QUADRATIC_WEDGE (15 pts) or VTK_QUADRATIC_PYRAMID (13 pts)
cells. femorph-solver’s
Wedge15 and
Pyr13 kernels
either run dedicated apex-singular shape functions ("bedrosian",
default) or restore the 20-slot layout and fold the 60 × 60 result
back down to 45 × 45 (wedge) or 39 × 39 (pyramid) via a constant
indicator matrix \(T\):
The hex-fold path is bit-for-bit identical to what the global assembler would compute if the element record had 20 node IDs with repeats; it is retained as a parity option for foreign decks that emit degenerate hexes natively.
Validated on a 5 166-element rotor-blade mesh (3 069 regular hex + 497 pyramid + 19 wedge + 1 581 tet): 12-mode modal agreement with established commercial solvers to 0.08 % relative on the same mesh (the residual is the plain-Gauss vs B-bar offset, not the degenerate dispatch).
Integration flag#
Name |
Description |
|---|---|
|
2×2×2 Gauss (8 points). 6 hourglass modes per element; no stabilisation yet — caller’s responsibility at mesh level (a connected mesh with shared midsides absorbs the hourglass modes automatically). |
|
3×3×3 Gauss (27 points). Fully ranked single-element \(K_e\) at ~3.4× the per-element flop cost. |
from femorph_solver import ELEMENTS
model.assign(
ELEMENTS.HEX20(integration="reduced"), # or "full" for 3×3×3 Gauss
{"EX": 2.1e11, "PRXY": 0.30},
)
Validation#
Rigid-body modes. 54 non-zero elastic modes + 6 rigid-body zeros
on an undistorted hex ("full"); 48 + 12 (6 rigid + 6 hourglass)
for "reduced".
Degenerate-shape correctness. WEDGE15 on a 15-node wedge and
PYR13 on a 13-node pyramid each have exactly 6 rigid-body zero
eigenvalues on the reference element.
Cross-solver migration safety. Plain modal on a quarter-arc blade sector (fixed root, 12 lowest modes) agrees with established commercial solvers to 0.08 % relative on the same mesh.
API reference#
- class femorph_solver.elements.hex20.Hex20[source]#
Bases:
ElementBase- static eel_batch(coords: ndarray, u_e: ndarray, material: dict[str, float] | None = None) ndarray | None[source]#
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.materialis accepted for signature uniformity with plane kernels but is unused.
- class femorph_solver.elements._wedge15_pyr13.Wedge15[source]#
Bases:
ElementBase15-node degenerate Hex20 wedge (K=L, O=P collapse).
Dedicated quadratic serendipity kernel with shape functions in area coords
(ξ₁, ξ₂, ζ), stiffness at 9-pt Gauss (3-pt triangle × 3-pt line), consistent mass at 21-pt Gauss (7-pt triangle × 3-pt line).
- class femorph_solver.elements._wedge15_pyr13.Pyr13[source]#
Bases:
ElementBase13-node degenerate Hex20 pyramid.
Default: dedicated Bedrosian apex-singular serendipity kernel integrated with a 2×2×2 Duffy collapsed-hex Gauss rule. The 2×2×2 rule (not 3×3×3) is the under-integrated quadrature that matches the foreign-deck convention; the higher-order
"consistent"rule is mathematically more accurate but ~13 % apart in Frobenius norm on a single element.The hex-fold wrapper (
T^T · K_20-node · T) is retained asmaterial["_PYR13_INTEGRATION"] = "hex_fold"for parity tests against foreign decks that emit hex-folded pyramids.
References#
Zienkiewicz, O. C., Taylor, R. L., and Zhu, J. Z. (2013). The Finite Element Method: Its Basis and Fundamentals, 7th ed. Butterworth-Heinemann. https://doi.org/10.1016/C2009-0-24909-9
Bathe, K.-J. (2014). Finite Element Procedures, 2nd ed. Prentice-Hall / Klaus-Jürgen Bathe. Ch. 5 (isoparametric elements). https://www.klausjurgenbathe.com/fepbook/
Cook, R. D., Malkus, D. S., Plesha, M. E., and Witt, R. J. (2001). Concepts and Applications of Finite Element Analysis, 4th ed. Wiley. ISBN 978-0-471-35605-9. https://www.wiley.com/en-us/Concepts+and+Applications+of+Finite+Element+Analysis%2C+4th+Edition-p-9780471356059