TET10 reference geometry — corner + mid-edge nodes, Keast points

Renders the 10-node quadratic tetrahedron’s reference element \(\{(\xi_1, \xi_2, \xi_3) : \xi_i \ge 0, \sum_i \xi_i \le 1\}\):

• 4 corner nodes at the unit-tet vertices.

• 6 mid-edge nodes at the centre of every edge.

• Keast (1986) 4-point quadrature rule used for both \(K^e\) and \(M^e\) — degree-2 exact, weights \(1/24\).

Shape functions (volume coordinates \(L_i\)):

\[N^c_i = L_i (2 L_i - 1), \quad i = 1, \ldots, 4 \qquad\text{and}\qquad N^m_{ij} = 4 L_i L_j\]

for the four corners and six mid-edges respectively (Cook 2002 Table 6.5-1).

Keast 4-point rule on the unit tet:

\[(\xi_1, \xi_2, \xi_3) = \left(\tfrac{1 \pm 3\alpha}{4}, \tfrac{1 \mp \alpha}{4}, \tfrac{1 \mp \alpha}{4}\right), \quad \alpha = \tfrac{1}{\sqrt{5}}, \quad w_q = \tfrac{1}{24}.\]

References

• Keast, P. (1986) “Moderate-degree tetrahedral quadrature formulas,” Computer Methods in Applied Mechanics and Engineering 55, 339–348.

• Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, §9.2–§9.3.

• Hughes, T. J. R. (2000) The Finite Element Method, Dover, §3.8.

Implementation: femorph_solver.elements.tet10.Tet10.

from __future__ import annotations

import numpy as np
import pyvista as pv

Node layout — 4 corners + 6 mid-edges

corners = np.array(
    [
        [0.0, 0.0, 0.0],
        [1.0, 0.0, 0.0],
        [0.0, 1.0, 0.0],
        [0.0, 0.0, 1.0],
    ]
)
mid_edges = np.array(
    [
        # ordering matches VTK_QUADRATIC_TETRA: edges 0-1, 1-2, 2-0,
        # 0-3, 1-3, 2-3.
        0.5 * (corners[0] + corners[1]),
        0.5 * (corners[1] + corners[2]),
        0.5 * (corners[2] + corners[0]),
        0.5 * (corners[0] + corners[3]),
        0.5 * (corners[1] + corners[3]),
        0.5 * (corners[2] + corners[3]),
    ]
)
all_nodes = np.vstack([corners, mid_edges])

cells = np.hstack([[10], np.arange(10, dtype=np.int64)])
cell_types = np.array([pv.CellType.QUADRATIC_TETRA], dtype=np.uint8)
ref_tet = pv.UnstructuredGrid(cells, cell_types, all_nodes)

Keast 4-point integration points

alpha = 1.0 / np.sqrt(5.0)
plus = (1.0 + 3.0 * alpha) / 4.0
minus = (1.0 - alpha) / 4.0
keast = np.array(
    [
        [plus, minus, minus],
        [minus, plus, minus],
        [minus, minus, plus],
        [minus, minus, minus],
    ]
)

Render the figure

plotter = pv.Plotter(off_screen=True, window_size=(640, 480))
plotter.add_mesh(
    ref_tet,
    style="wireframe",
    color="black",
    line_width=2,
    label="reference TET10",
)
plotter.add_points(
    corners,
    render_points_as_spheres=True,
    point_size=18,
    color="black",
    label="corner nodes (4)",
)
plotter.add_points(
    mid_edges,
    render_points_as_spheres=True,
    point_size=12,
    color="#1f77b4",
    label="mid-edge nodes (6)",
)
plotter.add_points(
    keast,
    render_points_as_spheres=True,
    point_size=16,
    color="#d62728",
    label="Keast 4-pt (K + M)",
)
plotter.add_axes(line_width=4, color="black")
plotter.view_isometric()
plotter.camera.zoom(1.1)
plotter.add_legend(face=None, size=(0.22, 0.14), bcolor="white")
plotter.show()

Sanity: Keast weights sum to the unit-tet volume = 1/6

w_q = 1.0 / 24.0
sum_w = 4 * w_q
assert abs(sum_w - 1.0 / 6.0) < 1e-15
print(f"Keast 4-pt: Σ w_q = {sum_w:.6f}; unit-tet volume = {1 / 6:.6f}.  Match.")

Keast 4-pt: Σ w_q = 0.166667; unit-tet volume = 0.166667.  Match.

Verify the partition of unity at the Keast points

L = np.column_stack([1.0 - keast.sum(axis=1), keast[:, 0], keast[:, 1], keast[:, 2]])
# Corner shape functions L_i (2 L_i - 1).
N_corner = L * (2.0 * L - 1.0)
# Mid-edge shape functions 4 L_i L_j over the 6 edges.
edge_pairs = [(0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)]
N_mid = np.column_stack([4.0 * L[:, i] * L[:, j] for i, j in edge_pairs])
sums = N_corner.sum(axis=1) + N_mid.sum(axis=1)
print("partition-of-unity Σ N_i at each Keast point:")
for q, s in zip(keast, sums, strict=True):
    print(f"  ξ₁={q[0]:+.4f}  ξ₂={q[1]:+.4f}  ξ₃={q[2]:+.4f}  Σ={s:.15f}")
np.testing.assert_allclose(sums, 1.0, atol=1.0e-14)
print("OK — all 10 shape functions sum to 1 at every Keast point.")

partition-of-unity Σ N_i at each Keast point:
  ξ₁=+0.5854  ξ₂=+0.1382  ξ₃=+0.1382  Σ=1.000000000000000
  ξ₁=+0.1382  ξ₂=+0.5854  ξ₃=+0.1382  Σ=1.000000000000000
  ξ₁=+0.1382  ξ₂=+0.1382  ξ₃=+0.5854  Σ=1.000000000000000
  ξ₁=+0.1382  ξ₂=+0.1382  ξ₃=+0.1382  Σ=1.000000000000000
OK — all 10 shape functions sum to 1 at every Keast point.