"""Hex8 — 8-node linear hexahedral solid.
Tri-linear shape functions on the reference cube ``ξ, η, ζ ∈ [-1, 1]``
with 2×2×2 Gauss integration. For each node ``i = 0…7`` with natural
coordinates ``(ξᵢ, ηᵢ, ζᵢ) ∈ {±1}³``:
Nᵢ(ξ, η, ζ) = (1/8)(1 + ξᵢ ξ)(1 + ηᵢ η)(1 + ζᵢ ζ)
Node ordering matches VTK_HEXAHEDRON — the eight cube corners
traversed ``(-,-,-), (+,-,-), (+,+,-), (-,+,-), (-,-,+),
(+,-,+), (+,+,+), (-,+,+)``.
Stiffness / mass
----------------
K_e = ∫_V Bᵀ C B dV ≈ Σ_g Bᵀ(ξ_g) C B(ξ_g) det J(ξ_g)
M_e = ρ ∫_V Nᵀ N dV ≈ Σ_g ρ Nᵀ(ξ_g) N(ξ_g) det J(ξ_g)
with isotropic linear elasticity in Voigt form (engineering shears):
C = λ·[111,000] ⊕ 2G·[diag(1,1,1,½,½,½)]
where ``λ = Eν/((1+ν)(1-2ν))`` and ``G = E/(2(1+ν))``.
Real constants
--------------
None required — geometry comes from node coordinates.
Material
--------
EX, PRXY : elastic moduli (both mandatory)
DENS : density (required for M_e)
Lumped mass
-----------
Row-sum (HRZ) lumping distributes each row of the consistent mass onto
its diagonal. For 8-node hex this gives m_total / 8 per corner node —
the standard and physically correct lumped form.
Integration flag
----------------
``material["_HEX8_INTEGRATION"]`` selects the stiffness formulation:
* ``"full"`` (default) — 2×2×2 Gauss full integration with Hughes
volumetric B-bar. The dilatation part of ``B`` at each Gauss point
is replaced by its element-volume-weighted average so
near-incompressible regimes do not lock.
* ``"plain_gauss"`` — plain 2×2×2 Gauss without the volumetric
correction. Useful for cross-checks against codes that don't ship
B-bar (e.g. the scikit-fem cross-FEM verification fixture); not
recommended for production runs because it locks under
near-incompressibility.
* ``"enhanced_strain"`` — Simo-Rifai 9-parameter enhanced assumed
strain (EAS) with Pian-Tong Jacobian scaling. Suppresses shear
locking in bending-dominated regular hex meshes.
References
----------
* Shape functions — trilinear Lagrange on the 8-node hexahedron:
Zienkiewicz, O.C. and Taylor, R.L. *The Finite Element Method:
Its Basis and Fundamentals*, 7th ed., Butterworth-Heinemann,
2013, Ch. 6. Equivalent form in Cook, Malkus, Plesha, Witt,
*Concepts and Applications of Finite Element Analysis*, 4th ed.,
Wiley, 2002, §6.2.
* 2×2×2 Gauss quadrature rule (product rule, 3 integration
points → 8 product points): Zienkiewicz & Taylor §5.10 Table 5.3.
* **B-bar formulation** (``"full"``, default — volumetric
dilatation replaced by its element-volume-weighted mean to
avoid near-incompressible locking): Hughes, T.J.R. *The Finite
Element Method: Linear Static and Dynamic Finite Element
Analysis*, Dover, 2000 (reprint of Prentice-Hall 1987), §4.4
("B-bar method"). Original formulation in Hughes, T.J.R.,
"Generalization of selective integration procedures to
anisotropic and nonlinear media," *IJNME* 15 (1980), pp.
1413–1418.
* **Enhanced assumed strain** (``"enhanced_strain"`` — 9 internal
parameters, Pian-Tong Jacobian scaling): Simo, J.C. and Rifai,
M.S., "A class of mixed assumed strain methods and the method
of incompatible modes," *IJNME* 29 (1990), pp. 1595–1638.
* Row-sum / HRZ lumped mass: Hinton, E., Rock, T. and Zienkiewicz,
O.C., "A note on mass lumping and related processes in the finite
element method," *Earthquake Engng. Struct. Dyn.* 4 (1976), pp.
245–249. Zienkiewicz & Taylor §16.2.2.
"""
from __future__ import annotations
import numpy as np
from femorph_solver.elements._base import ElementBase
from femorph_solver.elements._registry import register
try: # Optional C++ kernel — matches the Python reference to ~1e-15.
from femorph_solver import _elements as _ext
except ImportError: # pragma: no cover - extension always built in CI
_ext = None
# Natural coordinates of the 8 corners, in VTK_HEXAHEDRON order.
_XI_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=np.float64,
)
_GP = 1.0 / np.sqrt(3.0)
_GAUSS_POINTS = np.array(
[
[-_GP, -_GP, -_GP],
[+_GP, -_GP, -_GP],
[+_GP, +_GP, -_GP],
[-_GP, +_GP, -_GP],
[-_GP, -_GP, +_GP],
[+_GP, -_GP, +_GP],
[+_GP, +_GP, +_GP],
[-_GP, +_GP, +_GP],
],
dtype=np.float64,
)
_GAUSS_WEIGHTS = np.ones(8, dtype=np.float64)
def _shape_and_derivs(xi: float, eta: float, zeta: float) -> tuple[np.ndarray, np.ndarray]:
"""Return ``(N (8,), dN/dξ (8, 3))`` for the tri-linear cube."""
s = _XI_SIGNS
one_p_xi = 1.0 + s[:, 0] * xi
one_p_eta = 1.0 + s[:, 1] * eta
one_p_zeta = 1.0 + s[:, 2] * zeta
N = 0.125 * one_p_xi * one_p_eta * one_p_zeta
dN = np.empty((8, 3), dtype=np.float64)
dN[:, 0] = 0.125 * s[:, 0] * one_p_eta * one_p_zeta
dN[:, 1] = 0.125 * s[:, 1] * one_p_xi * one_p_zeta
dN[:, 2] = 0.125 * s[:, 2] * one_p_xi * one_p_eta
return N, dN
def _b_matrix(dN_dx: np.ndarray) -> np.ndarray:
"""Return the 6×24 strain-displacement matrix from ``dN/dx (8, 3)``."""
B = np.zeros((6, 24), dtype=np.float64)
for i in range(8):
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×24 shape-function matrix ``Φ`` with three rows per DOF."""
Phi = np.zeros((3, 24), dtype=np.float64)
for i in range(8):
Phi[0, 3 * i + 0] = N[i]
Phi[1, 3 * i + 1] = N[i]
Phi[2, 3 * i + 2] = N[i]
return Phi
def _elastic_matrix(E: float, nu: float) -> np.ndarray:
"""Isotropic Voigt C (engineering shears)."""
lam = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
G = E / (2.0 * (1.0 + nu))
C = np.zeros((6, 6), dtype=np.float64)
for i in range(3):
for j in range(3):
C[i, j] = lam
C[i, i] += 2.0 * G
for i in (3, 4, 5):
C[i, i] = G
return C
def _voigt_strain_transform(A: np.ndarray) -> np.ndarray:
"""6×6 Voigt transform of the tensor map ``ε ↦ A ε Aᵀ`` (engineering shears).
With Voigt order ``[ε_xx, ε_yy, ε_zz, γ_xy, γ_yz, γ_xz]`` and ``A``
any 3×3 matrix, this returns ``T`` such that the Voigt vector of
``B = A ε Aᵀ`` equals ``T @ ε_Voigt``. Used to push the natural-
frame enhanced strain through the element-centroid Jacobian
``J_0⁻ᵀ`` into the physical frame for Pian-Tong EAS.
"""
a, b, c = A[0]
d, e, f = A[1]
g, h, i = A[2]
T = np.array(
[
[a * a, b * b, c * c, a * b, b * c, a * c],
[d * d, e * e, f * f, d * e, e * f, d * f],
[g * g, h * h, i * i, g * h, h * i, g * i],
[2 * a * d, 2 * b * e, 2 * c * f, a * e + b * d, b * f + c * e, a * f + c * d],
[2 * d * g, 2 * e * h, 2 * f * i, d * h + e * g, e * i + f * h, d * i + f * g],
[2 * a * g, 2 * b * h, 2 * c * i, a * h + b * g, b * i + c * h, a * i + c * g],
],
dtype=np.float64,
)
return T
def _eas_natural(xi: float, eta: float, zeta: float) -> np.ndarray:
"""Simo-Rifai 9-mode enhanced strain in natural coords (6×9 Voigt).
Three diagonal modes (ξ, η, ζ) for ε_xx, ε_yy, ε_zz plus two cross
modes each for γ_xy, γ_yz, γ_xz. Each column integrates to zero on
the biunit cube, which combined with the Pian-Tong ``j_0/j`` scaling
guarantees the constant-stress patch test.
"""
M = np.zeros((6, 9), dtype=np.float64)
M[0, 0] = xi
M[1, 1] = eta
M[2, 2] = zeta
M[3, 3] = xi
M[3, 4] = eta
M[4, 5] = eta
M[4, 6] = zeta
M[5, 7] = xi
M[5, 8] = zeta
return M
def _ke_plain_gauss(coords: np.ndarray, E: float, nu: float) -> np.ndarray:
"""Plain 2×2×2 Gauss stiffness — no B-bar, no EAS. Exposed so
femorph-solver can compare apples-to-apples with codes that lack
B-bar (scikit-fem, older academic kernels). Locks under
near-incompressibility on distorted meshes; not for production use.
"""
C = _elastic_matrix(E, nu)
K = np.zeros((24, 24), dtype=np.float64)
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"Hex8 (plain_gauss): 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 8-node hex (default ``"full"`` integration).
The volumetric part of ``B`` at each Gauss point is replaced by its
element-volume-weighted mean so near-incompressible regimes do not
lock. On a regular hex (unit cube) the dilatation is already
constant, so B̄ ≡ B and this reduces to plain 2×2×2 Gauss.
"""
C = _elastic_matrix(E, nu)
# First pass: gather B, the volumetric-operator row ``b`` (= trace
# contribution), and the Gauss-point volumes. ``b`` is the 1×24 row
# such that ``θ(gp) = b · u_e`` (θ = ∂u/∂x + ∂v/∂y + ∂w/∂z).
B_all: list[np.ndarray] = []
b_all: list[np.ndarray] = []
dv_all: list[float] = []
V_e = 0.0
b_vol_sum = np.zeros(24, dtype=np.float64)
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"Hex8 (B-bar): non-positive Jacobian {det_J:.3e} — check node ordering"
)
dN_dx = np.linalg.solve(J, dN_dxi.T).T # (8, 3)
B_all.append(_b_matrix(dN_dx))
b = np.empty(24, dtype=np.float64)
b[0::3] = dN_dx[:, 0]
b[1::3] = dN_dx[:, 1]
b[2::3] = dN_dx[:, 2]
b_all.append(b)
dv = det_J * w
dv_all.append(dv)
V_e += dv
b_vol_sum += b * dv
b_bar = b_vol_sum / V_e
# Second pass: form B̄(gp) = B(gp) + (1/3) · [1,1,1,0,0,0]ᵀ · (b̄ − b(gp))
# and accumulate Kₑ = Σ B̄ᵀ C B̄ · dv.
e3 = np.array([1.0, 1.0, 1.0, 0.0, 0.0, 0.0]) / 3.0
K = np.zeros((24, 24), dtype=np.float64)
for B, b, dv in zip(B_all, b_all, dv_all):
B_bar = B + np.outer(e3, b_bar - b)
K += B_bar.T @ C @ B_bar * dv
return K
def _ke_eas(coords: np.ndarray, E: float, nu: float) -> np.ndarray:
"""Simo-Rifai E9 enhanced-strain stiffness for 8-node hex.
Static condensation: ``K_e = K_uu − K_uα K_αα⁻¹ K_uαᵀ`` where the
enhanced-strain contribution is added as 9 internal α parameters
that are eliminated element-locally. See the module docstring for
the formulation flag that selects this path.
"""
C = _elastic_matrix(E, nu)
_, dN_dxi_0 = _shape_and_derivs(0.0, 0.0, 0.0)
J_0 = dN_dxi_0.T @ coords
det_J_0 = float(np.linalg.det(J_0))
if det_J_0 <= 0:
raise ValueError(f"Hex8 (EAS): non-positive centroid Jacobian {det_J_0:.3e}")
# Voigt transform of ``J_0⁻¹`` (not ``J_0⁻ᵀ``). The strain pullback
# from natural to physical coords is ε_phys = J₀⁻¹ ε_nat J₀⁻ᵀ, which
# equals J₀⁻ᵀ ε_nat J₀⁻¹ only when J₀ is symmetric — so the unsymmetric
# form is the correct one and is verified against analytical reference
# solutions on skewed / trapezoidal / curved hexes.
T_0 = _voigt_strain_transform(np.linalg.inv(J_0))
K_uu = np.zeros((24, 24), dtype=np.float64)
K_ua = np.zeros((24, 9), dtype=np.float64)
K_aa = np.zeros((9, 9), dtype=np.float64)
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"Hex8 (EAS): non-positive Jacobian {det_J:.3e} — check node ordering")
dN_dx = np.linalg.solve(J, dN_dxi.T).T
B = _b_matrix(dN_dx)
M_nat = _eas_natural(*gp)
E_phys = (det_J_0 / det_J) * T_0 @ M_nat
dv = det_J * w
CB = C @ B
K_uu += B.T @ CB * dv
K_ua += B.T @ (C @ E_phys) * dv
K_aa += E_phys.T @ (C @ E_phys) * dv
return K_uu - K_ua @ np.linalg.solve(K_aa, K_ua.T)
[docs]
@register
class Hex8(ElementBase):
name = "HEX8"
n_nodes = 8
n_dof_per_node = 3
vtk_cell_type = 12 # VTK_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 != (8, 3):
raise ValueError(f"Hex8 expects (8, 3) coords, got {coords.shape}")
E = float(material["EX"])
nu = float(material["PRXY"])
formulation = material.get("_HEX8_INTEGRATION", "full")
if formulation == "enhanced_strain":
return _ke_eas(coords, E, nu)
if formulation == "plain_gauss":
# Plain 2×2×2 Gauss — apples-to-apples with scikit-fem etc.
# Prefer the C++ kernel when available; Python fallback kept
# so the path works even without the extension built.
if _ext is not None:
return _ext.solid185_ke(coords, E, nu)
return _ke_plain_gauss(coords, E, nu)
if formulation != "full":
raise ValueError(
f"HEX8: unknown _HEX8_INTEGRATION {formulation!r}; "
"expected 'full', 'plain_gauss', or 'enhanced_strain'"
)
# Hughes B-bar — the default femorph-solver formulation.
if _ext is not None:
return _ext.solid185_ke_bbar(coords, E, nu)
return _ke_bbar(coords, E, nu)
[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 != (8, 3):
raise ValueError(f"Hex8 expects (8, 3) coords, got {coords.shape}")
rho = float(material["DENS"])
if _ext is not None:
return _ext.solid185_me(coords, rho, lumped)
M = np.zeros((24, 24), dtype=np.float64)
for gp, w in zip(_GAUSS_POINTS, _GAUSS_WEIGHTS):
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
[docs]
@staticmethod
def ke_batch(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
) -> np.ndarray | None:
"""Vectorized Ke over ``(n_elem, 8, 3)`` coords.
Routes to the C++ B-bar kernel on the default path, to plain
Gauss when the ``plain_gauss`` flag is set, and returns
``None`` for EAS (Python-only for now) so the assembler falls
back to per-element ``ke()`` dispatch.
"""
if _ext is None:
return None
formulation = material.get("_HEX8_INTEGRATION", "full")
if formulation == "enhanced_strain":
return None
coords = np.ascontiguousarray(coords, dtype=np.float64)
E = float(material["EX"])
nu = float(material["PRXY"])
if formulation == "plain_gauss":
return _ext.solid185_ke_batch(coords, E, nu)
if formulation != "full":
raise ValueError(
f"HEX8: unknown _HEX8_INTEGRATION {formulation!r}; "
"expected 'full', 'plain_gauss', or 'enhanced_strain'"
)
return _ext.solid185_ke_bbar_batch(coords, E, nu)
@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
coords = np.ascontiguousarray(coords, dtype=np.float64)
return _ext.solid185_me_batch(coords, float(material["DENS"]), lumped)
[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`` is ``(n_elem, 8, 3)`` and ``u_e`` is ``(n_elem, 24)`` —
the element DOF vector laid out as ``[ux0, uy0, uz0, ux1, …]``.
Returns ``(n_elem, 8, 6)`` with Voigt components
``[εxx, εyy, εzz, γxy, γyz, γxz]`` (engineering shears), or
``None`` if the C++ extension is unavailable.
``material`` is accepted for signature uniformity with plane
kernels but is unused by the 3D-solid path (full 6-component
Voigt strain is recovered directly from ``u_e``).
"""
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.solid185_strain_batch(coords, u_e)