TET10 — 10-node quadratic tetrahedron#
Kinematics. 10-node quadratic tetrahedron (4 corners + 6 mid-edges); three translational DOFs per node (\(u_x, u_y, u_z\)); 30 DOFs per element.
Stiffness. \(K_e = \int_{\Omega_e} B^\top C B\,\mathrm{d}V\) integrated on the reference tetrahedron with the 4-point Keast rule.
Mass. Consistent \(M_e = \int \rho\, N^\top N\,\mathrm{d}V\) with a 15-point Keast rule (the quadratic-shape-function product \(N^\top N\) reaches polynomial degree 4 and needs the higher rule for rank).
Volume coordinates and shape functions#
The reference tetrahedron is conveniently described by volume coordinates (barycentric coordinates) \(L_1, L_2, L_3, L_4\) satisfying \(\sum_i L_i = 1\). Switching between volume coordinates and a Cartesian natural triple is a linear map; we use \((\xi, \eta, \zeta)\) with
Quadratic serendipity shape functions [Bathe2014] [ZTZ2013]:
Corner \(i\) (\(L_i = 1\) at corner \(i\)):
\[N_i = L_i\,(2 L_i - 1)\]Mid-edge between corners \(i\) and \(j\):
\[N_{ij} = 4\, L_i\, L_j\]
Node ordering matches VTK_QUADRATIC_TETRA: corners \(I, J, K,
L\) (indices 0-3), then mid-edges I-J, J-K, K-I, I-L, J-L, K-L
(indices 4-9).
Jacobian and \(B\) are the usual 3-D Voigt quantities. The isotropic constitutive \(C\) is identical to HEX8’s.
Keast quadrature on a tetrahedron#
Tetrahedral quadrature rules are tabulated in Keast’s foundational 1986 paper [Keast1986]. Two rules are used:
Stiffness — the 4-point rule (exact for polynomials up to degree 2):
\[\int_{T} f \,\mathrm{d}V \approx \frac{V}{4} \sum_{g=1}^{4} f(\xi_g),\]with the four points at the tetrahedron-symmetric positions
\[(a, b, b, b),\ (b, a, b, b),\ (b, b, a, b),\ (b, b, b, a),\]where \(a = (5 + 3\sqrt{5})/20\) and \(b = (5 - \sqrt{5})/20\). On a straight-edged tet with mid-edge nodes at the exact midpoints the Jacobian is constant, \(B^\top C B\) reaches degree 2, and this 4-point rule integrates \(K_e\) exactly.
Mass — the 15-point Keast rule (degree 5). \(N^\top N\) reaches degree 4 (quadratic × quadratic), so the 4-point rule is rank-deficient; the 15-point rule restores full rank of \(M_e\) on a single element. The full rotor modal parity on a 5 166-element rotor confirmed we need degree-5 integration for the mass.
Why the 4-point rule suffices for stiffness#
The strain field from quadratic shape functions is linear in \((\xi, \eta, \zeta)\) on a constant-Jacobian straight-edged tet, so \(B^\top C B\) is quadratic. The 4-point Keast rule integrates quadratics exactly. On a curved-edge tet (mid-edges off-centre) the Jacobian varies and higher-degree rules can improve accuracy.
Numpy walk-through#
import numpy as np
# Keast 4-point stiffness rule on the reference tet.
a = (5 + 3 * np.sqrt(5)) / 20
b = (5 - np.sqrt(5)) / 20
gp = np.array([(a, b, b, b),
(b, a, b, b),
(b, b, a, b),
(b, b, b, a)]) # barycentric (L1, L2, L3, L4)
w = np.full(4, 1/24) # reference-tet volume is 1/6, weight = 1/6 * 1/4
# Volume coordinates at the 10 nodes.
Lvec = np.array([
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1],
[0.5, 0.5, 0, 0],
[0, 0.5, 0.5, 0],
[0.5, 0, 0.5, 0],
[0.5, 0, 0, 0.5],
[0, 0.5, 0, 0.5],
[0, 0, 0.5, 0.5],
])
def shape_and_derivs(L):
# Quadratic serendipity shape functions + dN/dL derivatives.
# Convert to (dN/dξ, dN/dη, dN/dζ) via dL/dξ = (-1,-1,-1), I_3 otherwise.
N = np.empty(10)
for i in range(4):
N[i] = L[i] * (2 * L[i] - 1)
edges = [(0,1),(1,2),(2,0),(0,3),(1,3),(2,3)]
for k, (i, j) in enumerate(edges):
N[4 + k] = 4 * L[i] * L[j]
# Derivatives w.r.t. natural coords (ξ, η, ζ) by chain rule —
# straightforward bookkeeping; see the full listing in the
# gallery example.
...
return N, dN
The one-cube patch test passes; \(K_e\) has 6 rigid-body zero eigenvalues and 24 non-zero elastic modes on a unit regular tet.
Formulation flag#
TET10 has no user-facing formulation switch in the present implementation: there is just the one 4-point stiffness / 15-point mass path. Reduced-integration / enhanced-strain variants are roadmap and will be added with their verification fixtures.
Validation#
Rigid-body modes. 6 zero eigenvalues + 24 elastic on an undistorted reference tet.
Patch test. A linear displacement field ramps exactly through every node of any distorted tet mesh.
Full-rotor mass rank. On the 5 166-element rotor-blade benchmark mesh the global \(M\) is positive-definite (strictly) — this was the diagnostic that caught the original 4-point rank deficiency in mass integration and motivated the upgrade to Keast-15.
API reference#
- class femorph_solver.elements.tet10.Tet10[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:(n_elem, 10, 3);u_e:(n_elem, 30). Returns(n_elem, 10, 6)— engineering-shear Voigt strain.materialis accepted for signature uniformity with plane kernels but is unused.
References#
Bathe, K.-J. (2014). Finite Element Procedures, 2nd ed. Prentice-Hall / Klaus-Jürgen Bathe. https://www.klausjurgenbathe.com/fepbook/
Zienkiewicz, O. C., Taylor, R. L., and Zhu, J. Z. (2013). The Finite Element Method: Its Basis and Fundamentals, 7th ed. Butterworth-Heinemann. https://doi.org/10.1016/C2009-0-24909-9
Keast, P. (1986). “Moderate-degree tetrahedral quadrature formulas.” Computer Methods in Applied Mechanics and Engineering 55 (3), 339–348. https://doi.org/10.1016/0045-7825(86)90059-9