HEX8 — 8-node hexahedral solid#

Kinematics. Trilinear isoparametric 8-node hexahedron; three translational DOFs per node (\(u_x, u_y, u_z\)); 24 DOFs per element.

Stiffness. Displacement-based continuum formulation \(K_e = \int_{\Omega_e} B^\top C B \, \mathrm{d}V\), integrated on the reference cube \([-1, 1]^3\) with 2×2×2 Gauss-Legendre.

Mass. Consistent \(M_e = \int_{\Omega_e} \rho N^\top N \, \mathrm{d}V\), same integration rule. Row-sum lumping available.

Formulations. Default "full" is Hughes volumetric B-bar on top of 2×2×2 Gauss; "enhanced_strain" is Simo–Rifai 9-parameter EAS with Pian–Tong Jacobian scaling; "plain_gauss" is plain 2×2×2 Gauss for cross-checks.

Kinematics and shape functions#

The formulation below is the standard isoparametric continuum brick documented in every FEM textbook [Bathe2014] [Hughes1987] [ZTZ2013] [Cook2001]; we give it here in a form that maps line-by-line to the implementation in femorph_solver.elements.hex8.

On the reference cube with natural coordinates \((\xi, \eta, \zeta) \in [-1, 1]^3\), trilinear shape functions

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

are associated with the eight corner nodes \((\xi_i, \eta_i, \zeta_i) \in \{\pm 1\}^3\). The isoparametric map takes the reference cube to the physical hex:

\[x(\xi, \eta, \zeta) = \sum_{i=1}^{8} N_i(\xi, \eta, \zeta)\, x_i, \qquad u(\xi, \eta, \zeta) = \sum_{i=1}^{8} N_i(\xi, \eta, \zeta)\, u_i.\]

The node-ordering convention matches VTK_HEXAHEDRON: corners of the bottom face (\(\zeta = -1\)) first in CCW order, then the top face.

Jacobian and strain-displacement matrix#

Derivatives of \(N_i\) with respect to physical coordinates are recovered via the Jacobian of the isoparametric map,

\[\frac{\partial N_i}{\partial x_k} = \bigl(J^{-1}\bigr)_{kj}\, \frac{\partial N_i}{\partial \xi_j}, \qquad J_{jk} = \frac{\partial x_k}{\partial \xi_j} = \sum_{i=1}^{8} \frac{\partial N_i}{\partial \xi_j}\, x_{i,k}.\]

Small-strain Voigt notation with engineering shear gives \(\varepsilon = [\varepsilon_{xx},\ \varepsilon_{yy},\ \varepsilon_{zz},\ \gamma_{xy},\ \gamma_{yz},\ \gamma_{xz}]^\top\) and the 6 × 24 strain-displacement matrix

\[\begin{split}B_i = \begin{bmatrix} \partial_x N_i & 0 & 0 \\ 0 & \partial_y N_i & 0 \\ 0 & 0 & \partial_z N_i \\ \partial_y N_i & \partial_x N_i & 0 \\ 0 & \partial_z N_i & \partial_y N_i \\ \partial_z N_i & 0 & \partial_x N_i \end{bmatrix}.\end{split}\]

Stacking \(B = [B_1 \,|\, B_2 \,|\, \dots \,|\, B_8]\) produces the full 6 × 24 block.

Isotropic constitutive matrix#

For linear isotropic elasticity with Lamé parameters \(\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}\) and shear modulus \(\mu = G = \frac{E}{2(1+\nu)}\),

\[\begin{split}C = \begin{bmatrix} \lambda + 2\mu & \lambda & \lambda & 0 & 0 & 0 \\ \lambda & \lambda + 2\mu & \lambda & 0 & 0 & 0 \\ \lambda & \lambda & \lambda + 2\mu & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \end{bmatrix}.\end{split}\]

Gauss-Legendre quadrature on \([-1, 1]^3\)#

The 2 × 2 × 2 tensor-product rule uses 8 points at \((\pm 1/\sqrt{3}, \pm 1/\sqrt{3}, \pm 1/\sqrt{3})\) with unit weights. The rule is exact for polynomial integrands up to degree 3 in each direction — enough to integrate the \(B^\top C B\) product exactly on undistorted hexes. Physical volume integration picks up the determinant of \(J\):

\[K_e = \sum_{g=1}^{8} B(\xi_g)^\top C \, B(\xi_g)\, \lvert J(\xi_g) \rvert\, w_g.\]

Numpy walk-through (plain displacement formulation)#

The following pure-numpy listing assembles \(K_e\) for a unit cube made of generic steel (E = 210 GPa, ν = 0.30). Running it on a 24-DOF element produces a symmetric matrix with six exact rigid-body zero eigenvalues and eighteen non-zero elastic modes.

import numpy as np

# 2x2x2 Gauss-Legendre on [-1, 1]^3.
g = 1.0 / np.sqrt(3.0)
gp = 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)])
w = np.ones(8)

# Natural coordinates of the 8 nodes (VTK_HEXAHEDRON order).
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)],
                dtype=float)

def shape_and_derivs(xi, eta, zeta):
    N  = 0.125 * (1 + xi_n[:, 0] * xi) \
               * (1 + xi_n[:, 1] * eta) \
               * (1 + xi_n[:, 2] * zeta)
    dN = np.empty((8, 3))
    dN[:, 0] = 0.125 * xi_n[:, 0] * (1 + xi_n[:, 1] * eta) * (1 + xi_n[:, 2] * zeta)
    dN[:, 1] = 0.125 * xi_n[:, 1] * (1 + xi_n[:, 0] * xi)  * (1 + xi_n[:, 2] * zeta)
    dN[:, 2] = 0.125 * xi_n[:, 2] * (1 + xi_n[:, 0] * xi)  * (1 + xi_n[:, 1] * eta)
    return N, dN

# Unit cube nodes.
X = np.array([(0,0,0), (1,0,0), (1,1,0), (0,1,0),
              (0,0,1), (1,0,1), (1,1,1), (0,1,1)], dtype=float)

E, nu = 2.1e11, 0.30
lam = E*nu / ((1+nu)*(1-2*nu))
mu  = E / (2*(1+nu))
C = np.diag([lam+2*mu, lam+2*mu, lam+2*mu, mu, mu, mu]).astype(float)
C[:3, :3] += lam * (1 - np.eye(3))    # off-diagonal λ coupling of ε_ii

K = np.zeros((24, 24))
for (xi, eta, zeta), wi in zip(gp, w):
    _, dNdxi = shape_and_derivs(xi, eta, zeta)
    J = dNdxi.T @ X
    detJ = np.linalg.det(J)
    dNdx = np.linalg.solve(J, dNdxi.T).T
    B = np.zeros((6, 24))
    for i in range(8):
        col = 3 * i
        B[0, col + 0] = dNdx[i, 0]
        B[1, col + 1] = dNdx[i, 1]
        B[2, col + 2] = dNdx[i, 2]
        B[3, col + 0] = dNdx[i, 1]
        B[3, col + 1] = dNdx[i, 0]
        B[4, col + 1] = dNdx[i, 2]
        B[4, col + 2] = dNdx[i, 1]
        B[5, col + 0] = dNdx[i, 2]
        B[5, col + 2] = dNdx[i, 0]
    K += B.T @ C @ B * detJ * wi

eig = np.sort(np.linalg.eigvalsh(K))
# eig[:6]  ≈ 0          # six rigid-body modes
# eig[6:]  > 0          # 18 elastic modes

Volumetric locking and Hughes B-bar#

The plain displacement formulation above suffers from volumetric locking as \(\nu \to 1/2\) or when the strain field has a strong volumetric component that the trilinear kinematics cannot accommodate. Hughes’s selective reduced integration [Hughes1980] projects the volumetric strain component onto the element-constant subspace:

\[\begin{split}\bar{B}_i = B_i + \tfrac{1}{3} \begin{bmatrix} \bar{b}_x - b_{i,x} & \bar{b}_y - b_{i,y} & \bar{b}_z - b_{i,z} \\ \bar{b}_x - b_{i,x} & \bar{b}_y - b_{i,y} & \bar{b}_z - b_{i,z} \\ \bar{b}_x - b_{i,x} & \bar{b}_y - b_{i,y} & \bar{b}_z - b_{i,z} \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix},\end{split}\]

where \(b_{i,k} = \partial_{x_k} N_i\) and \(\bar{b}_k = \frac{1}{V_e} \int_{\Omega_e} b_{i,k} \, \mathrm{d}V\) is the volume-averaged derivative. In practice the correction is implemented in closed form as

\[K_{\bar B} = K_\text{plain} + K_B \bigl(\tfrac{1}{V_e}\, S S^\top - H\bigr), \qquad K_B = \lambda + \tfrac{2}{3} G,\]

with \(S = \sum_g b(\xi_g)\lvert J(\xi_g)\rvert w_g\) and \(H = \sum_g b(\xi_g)b(\xi_g)^\top \lvert J(\xi_g)\rvert w_g\), where \(b(\xi_g)\) is the 24-vector of volumetric derivatives at Gauss point \(\xi_g\). This is femorph-solver’s default “full” formulation.

Enhanced assumed strain (Simo–Rifai)#

An alternative cure for volumetric and shear locking replaces the compatible strain \(\varepsilon = B u_e\) with an enhanced decomposition

\[\tilde{\varepsilon} = B u_e + G(\xi)\, \alpha,\]

where \(\alpha\) is a set of per-element strain parameters and \(G(\xi)\) is a 9-column enhancement matrix constructed from hierarchic bubble functions [SimoRifai1990]. Static condensation eliminates \(\alpha\) at the element level, producing the enhanced stiffness

\[\tilde{K}_e = K_e - K_{e\alpha} \bigl(K_{\alpha\alpha}\bigr)^{-1} K_{\alpha e},\]

with

\[\begin{split}K_e &= \int B^\top C B \,\mathrm{d}V, \\ K_{e\alpha} &= \int B^\top C \, \Gamma(\xi) \,\mathrm{d}V, \\ K_{\alpha\alpha} &= \int \Gamma(\xi)^\top C \, \Gamma(\xi) \,\mathrm{d}V,\end{split}\]

where \(\Gamma(\xi) = j_0 \, J_0^{-1}\, G(\xi)\, / \lvert J(\xi)\rvert\) is the enhancement in current-configuration coordinates under the Pian–Tong \(j_0 / j\) scaling.

Reach the EAS path via ELEMENTS.HEX8(integration="enhanced_strain") on assign(); the kernel keys the choice on material["_HEX8_INTEGRATION"] internally.

Mass matrix#

Consistent mass uses the same integration rule:

\[M_e = \int_{\Omega_e} \rho \, N^\top N \,\mathrm{d}V = \sum_{g=1}^{8} \rho \, N(\xi_g)^\top N(\xi_g) \, \lvert J(\xi_g)\rvert \, w_g.\]

Row-sum lumping \(M_e^{\text{lump}} = \text{diag}\bigl(\sum_j M_{e,ij}\bigr)\) is available via mass_matrix with lumped=True.

Integration flag#

Name

Description

"full" (default)

2×2×2 Gauss + Hughes volumetric B-bar. Cures volumetric locking at near-incompressible Poisson ratios.

"enhanced_strain"

Simo–Rifai EAS with 9 enhancement modes and Pian–Tong \(j_0 / j\) scaling. Cures both shear and volumetric locking; recovers >95 % of the Euler–Bernoulli tip deflection on a slender cantilever with one element through the thickness (vs <50 % for plain 2×2×2).

"plain_gauss"

Vanilla 2×2×2 Gauss with no B-bar / EAS treatment. Locks volumetrically at near-incompressible Poisson and shear- locks on bending-dominated geometries — useful only for cross-solver comparisons with codes that don’t ship B-bar.

Pick the formulation at assign() time via the ELEMENTS namespace:

import femorph_solver as fs
from femorph_solver import ELEMENTS

model = fs.Model.from_grid(grid)
model.assign(
    ELEMENTS.HEX8(integration="enhanced_strain"),  # or "full" (default, B-bar)
    {"EX": 2.1e11, "PRXY": 0.30},
)

Validation#

Rigid-body modes. On any undistorted hex, \(K_e\) should have exactly six zero eigenvalues (three translations + three rotations) and eighteen non-zero elastic modes. Verify directly:

import numpy as np
from femorph_solver.elements import get

HEX8 = get("HEX8")
X = np.array([(0,0,0), (1,0,0), (1,1,0), (0,1,0),
              (0,0,1), (1,0,1), (1,1,1), (0,1,1)], dtype=float)
K = HEX8.ke(X, {"EX": 2.1e11, "PRXY": 0.30}, np.empty(0))
eig = np.sort(np.linalg.eigvalsh(K))
assert np.allclose(eig[:6], 0, atol=1e-4 * eig[-1])   # 6 rigid-body modes
assert (eig[6:] > 0).all()                            # 18 elastic modes

Cantilever tip deflection. For a slender cantilever (\(L \gg h\)) with tip load \(P\), Euler–Bernoulli theory gives \(u_\text{tip} = P L^3 / (3 E I)\). A single hex through the thickness reaches this limit to within a few percent under the EAS formulation; the plain-displacement formulation severely under-predicts due to shear locking. The MacNeal–Harder benchmark set [MacNealHarder1985] ships the canonical variants of this problem; the SOLID185 gallery runs the slender-cantilever version end-to-end.

Cross-solver migration safety. femorph-solver’s HEX8 matches established commercial solvers (default B-bar formulation) on the bundled quarter-arc sector to machine precision — see the HEX8 cyclic-modal migration check.

API reference#

class femorph_solver.elements.hex8.Hex8[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 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).

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

  • coords ((n_nodes, 3) float64)

  • material (mapping with property names as keys ("EX", "PRXY", "DENS", …))

  • real (1-D array of real constants)

static ke_batch(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray | None[source]#

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.

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#

Return the element mass matrix in global DOF ordering.

References#

[Bathe2014]

Bathe, K.-J. (2014). Finite Element Procedures, 2nd ed. Prentice-Hall / Klaus-Jürgen Bathe. Chapters 4 (isoparametric elements) and 6 (continuum). https://www.klausjurgenbathe.com/fepbook/

[Hughes1987]

Hughes, T. J. R. (1987). The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Prentice-Hall; Dover reprint 2000 (ISBN 0-486-41181-8). Chapter 3 (isoparametric formulation), §4.5 (selective reduced integration). https://store.doverpublications.com/0486411818.html

[Hughes1980]

Hughes, T. J. R. (1980). “Generalization of selective integration procedures to anisotropic and nonlinear media.” International Journal for Numerical Methods in Engineering 15, 1413–1418. https://doi.org/10.1002/nme.1620150914

[SimoRifai1990]

Simo, J. C. and Rifai, M. S. (1990). “A class of mixed assumed strain methods and the method of incompatible modes.” International Journal for Numerical Methods in Engineering 29, 1595–1638. https://doi.org/10.1002/nme.1620290802

[MacNealHarder1985]

MacNeal, R. H. and Harder, R. L. (1985). “A proposed standard set of problems to test finite element accuracy.” Finite Elements in Analysis and Design 1 (1), 3–20. https://doi.org/10.1016/0168-874X(85)90003-4

[ZTZ2013]

Zienkiewicz, O. C., Taylor, R. L., and Zhu, J. Z. (2013). The Finite Element Method: Its Basis and Fundamentals, 7th ed. Butterworth-Heinemann (Elsevier). https://doi.org/10.1016/C2009-0-24909-9

[Cook2001]

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