"""SOLID187 — 10-node quadratic tetrahedral solid.
Node ordering matches MAPDL and VTK_QUADRATIC_TETRA:
0-3 : corner nodes I, J, K, L
4 : mid-edge I-J
5 : mid-edge J-K
6 : mid-edge K-I
7 : mid-edge I-L
8 : mid-edge J-L
9 : mid-edge K-L
Reference tetrahedron with natural coordinates ``(ξ, η, ζ)``, volume
coordinates ``(L₁, L₂, L₃, L₄) = (1-ξ-η-ζ, ξ, η, ζ)``:
* **Corner** node ``i`` (Lᵢ = 1 at corner ``i``): Nᵢ = Lᵢ(2Lᵢ - 1)
* **Mid-edge** ``(i, j)``: N_{ij} = 4 Lᵢ Lⱼ
Both stiffness *and* mass use the 4-point Gauss rule on the reference
tet (degree-2 exact, weights ``1/24`` at the 4 symmetric Keast points).
This matches MAPDL exactly: a side-by-side comparison of the 30×30
element K and M against MAPDL's EMAT (after applying MAPDL's
kinematic constraint ``u_Q = (u_J+u_L)/2``, ``u_R = (u_K+u_L)/2`` on
the J-L and K-L midside slots) lands at ~7 × 10⁻¹³ relative on both
straight and skewed tets — i.e. bit-exact agreement with MAPDL's
internal 4-pt rule. At the assembled-global level the condensation is
automatic: eigenfrequencies of a plain-modal SOLID187 sector match
MAPDL's output to ~1 × 10⁻⁹ (down from ~6 × 10⁻⁴ when we used an
over-integrated 15-pt mass rule — the former "accuracy" was actually
drift away from MAPDL's lower-order formulation).
The 4-pt mass integrates a degree-4 integrand with a degree-2 rule, so
per-element M is rank 12 of 30 (PSD, not PD). That's fine: once two
elements share nodes the global M becomes PD, and lumped M is diagonal.
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 — 10-node quadratic Lagrange tetrahedron in
volume coordinates (``L₁, L₂, L₃, L₄``): Zienkiewicz, O.C. and
Taylor, R.L. *The Finite Element Method: Its Basis and
Fundamentals*, 7th ed., Butterworth-Heinemann, 2013, §8.8.2.
Also Cook, Malkus, Plesha, Witt, *Concepts and Applications of
Finite Element Analysis*, 4th ed., Wiley, 2002, §6.5 Table 6.5-1.
* **4-point Keast Gauss rule** (degree-2 exact on the unit
tetrahedron, weights ``1/24`` at the four symmetric points):
Keast, P., "Moderate-degree tetrahedral quadrature formulas,"
*Comput. Methods Appl. Mech. Engrg.* 55 (1986), pp. 339–348.
The 4-point rule integrates K exactly (degree-2 integrand from
``BᵀCB`` on quadratic shape functions) and M with rank
deficiency 12 of 30 that the global assembly absorbs.
* Kinematic midside constraints (``u_Q = ½(u_J + u_L)`` and
``u_R = ½(u_K + u_L)``) used for the MAPDL-EMAT condensation
comparison: standard degenerate-hex / degenerate-tet treatment in
Zienkiewicz, Taylor, and Zhu, *The Finite Element Method: Its
Basis and Fundamentals*, 7th ed., §8.8.1 discussion of collapsed
hex nodes.
* Voigt elastic matrix: Cook et al. §3.7.
MAPDL compatibility — specification source
------------------------------------------
* Ansys, Inc., *Ansys Mechanical APDL Element Reference*,
Release 2022R2, section "SOLID187 — 3-D 10-Node Tetrahedral
Structural Solid".
Short factual summary (paraphrased from that section):
* Topology: 10-node tetrahedron (4 corner + 6 mid-edge nodes);
3 translational DOFs per node; 30 DOFs per element.
* Integration rule: 4-point Gauss quadrature on the reference
tetrahedron for both K and M.
* Elastic material input: ``EX`` + ``PRXY``; density ``DENS``.
femorph-solver reimplements the kernel from Keast 1986 and the
standard volume-coordinate quadratic shape functions; the Ansys
Element Reference is consulted only as the compatibility
specification — what a MAPDL deck expects ``SOLID187`` 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
N_NODES = 10
N_DOF = 3 * N_NODES # 30
def _gauss_4pt() -> tuple[np.ndarray, np.ndarray]:
"""4-point Gauss rule on the reference tet (volume 1/6).
Exact for polynomials of degree 2 — sufficient for Ke on a
straight-edged tet with linear Jacobian (degree-2 integrand).
"""
a = (5.0 + 3.0 * np.sqrt(5.0)) / 20.0
b = (5.0 - np.sqrt(5.0)) / 20.0
pts = np.array(
[
[b, b, b],
[a, b, b],
[b, a, b],
[b, b, a],
],
dtype=np.float64,
)
wts = np.full(4, 1.0 / 24.0, dtype=np.float64)
return pts, wts
_GAUSS_POINTS, _GAUSS_WEIGHTS = _gauss_4pt()
def _shape_and_derivs(xi: float, eta: float, zeta: float) -> tuple[np.ndarray, np.ndarray]:
"""Return ``(N (10,), dN/dξ (10, 3))`` for the 10-node tet."""
L1 = 1.0 - xi - eta - zeta
L2 = xi
L3 = eta
L4 = zeta
N = np.array(
[
L1 * (2.0 * L1 - 1.0),
L2 * (2.0 * L2 - 1.0),
L3 * (2.0 * L3 - 1.0),
L4 * (2.0 * L4 - 1.0),
4.0 * L1 * L2,
4.0 * L2 * L3,
4.0 * L3 * L1,
4.0 * L1 * L4,
4.0 * L2 * L4,
4.0 * L3 * L4,
],
dtype=np.float64,
)
# Derivatives: ∂L₁/∂(ξ,η,ζ) = (-1,-1,-1); ∂L₂/∂ξ = 1, ∂L₃/∂η = 1, ∂L₄/∂ζ = 1
dN = np.empty((10, 3), dtype=np.float64)
# Corner nodes: Nᵢ = Lᵢ(2Lᵢ-1) → dNᵢ/dx = (4Lᵢ-1) · dLᵢ/dx
dN[0] = (4.0 * L1 - 1.0) * np.array([-1.0, -1.0, -1.0])
dN[1] = (4.0 * L2 - 1.0) * np.array([+1.0, 0.0, 0.0])
dN[2] = (4.0 * L3 - 1.0) * np.array([0.0, +1.0, 0.0])
dN[3] = (4.0 * L4 - 1.0) * np.array([0.0, 0.0, +1.0])
# Mid-edge nodes: N = 4 Lᵢ Lⱼ → dN/dx = 4(Lⱼ · dLᵢ/dx + Lᵢ · dLⱼ/dx)
dN[4] = 4.0 * (L2 * np.array([-1.0, -1.0, -1.0]) + L1 * np.array([+1.0, 0.0, 0.0]))
dN[5] = 4.0 * (L3 * np.array([+1.0, 0.0, 0.0]) + L2 * np.array([0.0, +1.0, 0.0]))
dN[6] = 4.0 * (L1 * np.array([0.0, +1.0, 0.0]) + L3 * np.array([-1.0, -1.0, -1.0]))
dN[7] = 4.0 * (L4 * np.array([-1.0, -1.0, -1.0]) + L1 * np.array([0.0, 0.0, +1.0]))
dN[8] = 4.0 * (L4 * np.array([+1.0, 0.0, 0.0]) + L2 * np.array([0.0, 0.0, +1.0]))
dN[9] = 4.0 * (L4 * np.array([0.0, +1.0, 0.0]) + L3 * np.array([0.0, 0.0, +1.0]))
return N, dN
def _b_matrix(dN_dx: np.ndarray) -> np.ndarray:
"""Return the 6×30 strain-displacement matrix from ``dN/dx (10, 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×30 shape-function matrix."""
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
[docs]
@register
class SOLID187(ElementBase):
name = "SOLID187"
n_nodes = N_NODES
n_dof_per_node = 3
vtk_cell_type = 24 # VTK_QUADRATIC_TETRA
[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"SOLID187 expects ({N_NODES}, 3) coords, got {coords.shape}")
E = float(material["EX"])
nu = float(material["PRXY"])
if _ext is not None:
return _ext.solid187_ke(coords, E, nu)
C = _elastic_matrix(E, nu)
K = np.zeros((N_DOF, N_DOF), 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"SOLID187: 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
[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"SOLID187 expects ({N_NODES}, 3) coords, got {coords.shape}")
rho = float(material["DENS"])
if _ext is not None:
return _ext.solid187_me(coords, rho, lumped)
M = np.zeros((N_DOF, N_DOF), 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
@staticmethod
def ke_batch(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
) -> np.ndarray | None:
if _ext is None:
return None
coords = np.ascontiguousarray(coords, dtype=np.float64)
return _ext.solid187_ke_batch(coords, float(material["EX"]), float(material["PRXY"]))
@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.solid187_me_batch(coords, float(material["DENS"]), bool(lumped))
[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, 10, 3)``; ``u_e``: ``(n_elem, 30)``.
Returns ``(n_elem, 10, 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.solid187_strain_batch(coords, u_e)