TET10 — 10-node quadratic tetrahedron ===================================== *MAPDL-compatible alias:* ``SOLID187``. **Kinematics.** 10-node quadratic tetrahedron (4 corners + 6 mid-edges); three translational DOFs per node (:math:`u_x, u_y, u_z`); 30 DOFs per element. **Stiffness.** :math:`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 :math:`M_e = \int \rho\, N^\top N\,\mathrm{d}V` with a 15-point Keast rule (the quadratic-shape-function product :math:`N^\top N` reaches polynomial degree 4 and needs the higher rule for rank). **MAPDL compatibility.** Reproduces SOLID187 at its default settings. .. contents:: On this page :local: :depth: 2 ---- Volume coordinates and shape functions -------------------------------------- The reference tetrahedron is conveniently described by *volume coordinates* (barycentric coordinates) :math:`L_1, L_2, L_3, L_4` satisfying :math:`\sum_i L_i = 1`. Switching between volume coordinates and a Cartesian natural triple is a linear map; we use :math:`(\xi, \eta, \zeta)` with .. math:: L_1 = 1 - \xi - \eta - \zeta, \quad L_2 = \xi, \quad L_3 = \eta, \quad L_4 = \zeta. Quadratic serendipity shape functions [Bathe2014]_ [ZTZ2013]_: * **Corner** :math:`i` (:math:`L_i = 1` at corner :math:`i`): .. math:: N_i = L_i\,(2 L_i - 1) * **Mid-edge** between corners :math:`i` and :math:`j`: .. math:: N_{ij} = 4\, L_i\, L_j Node ordering matches MAPDL and VTK_QUADRATIC_TETRA: corners :math:`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 :math:`B` are the usual 3-D Voigt quantities. The isotropic constitutive :math:`C` is identical to SOLID185'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): .. math:: \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 .. math:: (a, b, b, b),\ (b, a, b, b),\ (b, b, a, b),\ (b, b, b, a), where :math:`a = (5 + 3\sqrt{5})/20` and :math:`b = (5 - \sqrt{5})/20`. On a straight-edged tet with mid-edge nodes at the exact midpoints the Jacobian is constant, :math:`B^\top C B` reaches degree 2, and this 4-point rule integrates :math:`K_e` exactly. * **Mass** — the 15-point Keast rule (degree 5). :math:`N^\top N` reaches degree 4 (quadratic × quadratic), so the 4-point rule is rank-deficient; the 15-point rule restores full rank of :math:`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 :math:`(\xi, \eta, \zeta)` on a constant-Jacobian straight-edged tet, so :math:`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, but the default is adequate for the MAPDL-compatible path we target. Numpy walk-through ------------------ .. code-block:: python 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; :math:`K_e` has 6 rigid-body zero eigenvalues and 24 non-zero elastic modes on a unit regular tet. Formulation flag ---------------- SOLID187 has no user-facing formulation switch in the present implementation: there's just the one 4-point stiffness / 15-point mass path. KEYOPT variants (uniform reduced, enhanced-strain mixed-type, etc.) are documented in the MAPDL element reference but not yet exposed by femorph-solver. 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 CDB the global :math:`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 ------------- .. autoclass:: femorph_solver.elements.solid187.SOLID187 :members: References ---------- .. [Bathe2014] Bathe, K.-J. (2014). *Finite Element Procedures*, 2nd ed. Prentice-Hall / Klaus-Jürgen Bathe. https://www.klausjurgenbathe.com/fepbook/ .. [ZTZ2013] 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 .. [Keast1986] 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