"""Tet10 — 10-node quadratic tetrahedral solid.
Node ordering matches 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).
The 4-pt rule integrates K exactly (degree-2 integrand from ``BᵀCB``
on quadratic shape functions); the per-element M is degree-4 and
therefore rank 12 of 30 (PSD, not PD) under the same rule. That's
fine: once two elements share nodes the global M becomes PD, and
lumped M is diagonal regardless.
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.
* Voigt elastic matrix: Cook et al. §3.7.
"""
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.hex8 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 Tet10(ElementBase):
name = "TET10"
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"Tet10 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"TET10: 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"Tet10 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,
material: dict[str, float] | None = None,
) -> 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.
``material`` is accepted for signature uniformity with plane
kernels but is unused.
"""
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.solid187_strain_batch(coords, u_e)