"""SOLID186 — 20-node quadratic hexahedral solid (serendipity).
Node ordering matches MAPDL (and VTK_QUADRATIC_HEXAHEDRON): the 8 corner
nodes come first in SOLID185 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["_SOLID186_FORMULATION"]`` selects the stiffness integration
rule:
* ``"reduced"`` (default) — 2×2×2 Gauss-Legendre (8 points), matching
MAPDL's KEYOPT(2) = 0 uniform-reduced default. Single-element Ke has
6 extra zero-energy hourglass modes. On a connected mesh those modes
almost always die under the assembly — a hex whose midside is shared
with a neighbour has no free hourglass direction to excite — and the
resulting parity against MAPDL is ~3 × 10⁻⁵ on ptr.cdb vs. ~8 × 10⁻⁴
with full 3×3×3. 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), matching MAPDL's
KEYOPT(2) = 1 full-integration option. Ke is fully rank
(54 non-zero eigenvalues + 6 rigid-body zeros) on a single hex.
Mass integration is always 3×3×3 (not affected by the formulation flag).
Wedge / pyramid degenerate forms use the Irons 14-point rule on the
hex (see ``solid186_degenerate.py`` for the wedge / pyramid kernels);
the 14-point 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 MAPDL's KEYOPT(2) = 0
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.
MAPDL compatibility — specification source
------------------------------------------
* Ansys, Inc., *Ansys Mechanical APDL Element Reference*,
Release 2022R2, section "SOLID186 — 3-D 20-Node Structural
Solid".
Short factual summary (paraphrased from that section — facts are
not copyrightable; vendor text is not reproduced):
* Topology: 20-node hexahedron (8 corners + 12 mid-edge nodes);
3 translational DOFs per node; 60 DOFs per element. Degenerate
15-node wedge and 13-node pyramid forms dispatched to
``SOLID186W`` / ``SOLID186P`` kernels (see
``solid186_degenerate.py``).
* ``KEYOPT(2) = 0`` (default): "Uniform reduced integration" —
reproduced here under ``"reduced"``.
* ``KEYOPT(2) = 1``: "Full integration" (3×3×3 Gauss-Legendre) —
reproduced here under ``"full"``.
* Mass integration: 14-point Irons rule regardless of KEYOPT(2);
reproduced per the same Bedrosian 1992 + Irons 1971 references
cited in the Theory block above.
* Elastic material input: ``EX`` + ``PRXY``; density ``DENS``.
femorph-solver reimplements the kernel from the open references
above; the Ansys Element Reference is consulted only as the
compatibility specification — what a MAPDL deck expects
``SOLID186`` (and its degenerate forms) to do.
"""
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.solid185 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 MAPDL / 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 — MAPDL SOLID186 mass.
Verified bit-exact (~3 × 10⁻¹² rel Frobenius) against MAPDL's own
``Me`` in ptr.cdb's EMAT for SOLID186 KEYOPT(2)=0. This is **not**
the 2×2×2 rule used for URI stiffness and **not** the 3×3×3
consistent rule a textbook would prescribe: MAPDL specifically
drops the consistent mass's full midside coupling in favour of
Irons' (1971) lower-order 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"SOLID186: 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 (MAPDL KEYOPT(2)=1).
Same closed-form correction used in SOLID185'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"SOLID186 (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 SOLID186(ElementBase):
name = "SOLID186"
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"SOLID186 expects ({N_NODES}, 3) coords, got {coords.shape}")
E = float(material["EX"])
nu = float(material["PRXY"])
formulation = material.get("_SOLID186_FORMULATION", "reduced")
if formulation == "reduced":
# 2×2×2 Gauss, matching MAPDL KEYOPT(2)=0 (the MAPDL 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, matching MAPDL KEYOPT(2)=1. Slower and further
# from MAPDL's default; opt-in for locking-prone meshes.
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"SOLID186: unknown _SOLID186_FORMULATION {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"SOLID186 expects ({N_NODES}, 3) coords, got {coords.shape}")
rho = float(material["DENS"])
# MAPDL SOLID186 KEYOPT(2)=0 integrates Me with the 14-point Irons
# rule, not the 3×3×3 consistent rule a textbook would prescribe.
# Verified bit-exact (~3 × 10⁻¹² rel Frobenius) vs MAPDL's ptr.cdb
# EMAT on real distorted hexes. Callers can opt into the old
# 3×3×3 consistent rule via ``material["_SOLID186_MASS"] = "consistent"``.
mass_rule = material.get("_SOLID186_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"SOLID186: unknown _SOLID186_MASS {mass_rule!r}; "
f"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("_SOLID186_FORMULATION", "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("_SOLID186_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) -> 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.
"""
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)