fem: pyramid5 element (solver, quality, decks, VTK)

A five-node pyramid, the conforming transition between a hexahedron's face
and tetrahedra: Bedrosian's rational shape functions with a collapsed (Duffy)
Gauss rule, a certified Jacobian on the base square with exact volume, a
scaled Jacobian, C3D5 / CPYRAM / LS-DYNA degenerate solid, VTK_PYRAMID and
the binary mesh code 13. A distorted hex8-pyramid5-tet4 transition mesh
passes the patch test and uniaxial tension exactly.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo

crates/mesh/src/mirror.rs

@@ -23,6 +23,8 @@ pub fn reflected_order(kind: ElementKind) -> Option<&'static [usize]> {
         ElementKind::Tet10 => &[0, 2, 1, 3, 6, 5, 4, 7, 9, 8],
         // Each triangle's last two corners swap.
         ElementKind::Wedge6 => &[0, 2, 1, 3, 5, 4],
+        // The base's corners in reverse cyclic order; the apex stays.
+        ElementKind::Pyramid5 => &[0, 3, 2, 1, 4],
         // Each quadrilateral's corners in reverse cyclic order from its first.
         ElementKind::Hex8 => &[0, 3, 2, 1, 4, 7, 6, 5],
         ElementKind::Hex20 => &[0, 3, 2, 1, 4, 7, 6, 5, 11, 10, 9, 8, 15, 14, 13, 12, 16, 19, 18, 17],

crates/mesh/src/quality/jacobian.rs

@@ -19,6 +19,13 @@
 //!   lies on the side of face 1-2-3 that its right-hand normal points to.
 //! - wedge6: triangle 1-2-3, then 4-5-6 above it; positive when 4-5-6 lies on the side of 1-2-3's
 //!   right-hand normal.
+//! - pyramid5: base 1-2-3-4, then the apex 5 on the side of its right-hand normal. Its map is
+//!   rational in the pyramid's coordinates (ξ, η, ζ), but its Jacobian determinant is not: with
+//!   the base's bilinear map `B(s, t)` on the unit square (ξ = (2s − 1)(1 − ζ), likewise η),
+//!   the columns are `∂x/∂ξ = ½ B_s`, `∂x/∂η = ½ B_t` and `∂x/∂ζ = x₅ − B + (s − ½) B_s +
+//!   (t − ½) B_t`, all bilinear in (s, t) and constant along each line to the apex, so the
+//!   determinant is a polynomial of degree (3, 3) on the square and the element's volume is
+//!   4/3 of its integral there.
 //! - hex8/hex20: face 1-2-3-4, then 5-6-7-8 above it; mid-side nodes 9–12 on the edges of 1-2-3-4,
 //!   13–16 on those of 5-6-7-8, 17–20 on edges 1-5, 2-6, 3-7, 4-8.
 
@@ -213,12 +220,43 @@ fn geometry(kind: ElementKind, xyz: &[[f64; 3]]) -> Result<[Bernstein; 3]> {
     Ok(g)
 }
 
+/// The pyramid's Jacobian columns (see the module docs) as bilinear polynomials on the unit
+/// square of its base.
+fn pyramid_columns(xyz: &[[f64; 3]]) -> Result<[[Bernstein; 3]; 3]> {
+    let node = |i: usize| xyz.get(i).copied().ok_or_else(|| Error::Invalid("pyramid5 needs 5 nodes".into()));
+    let [p0, p1, p2, p3, apex] = [node(0)?, node(1)?, node(2)?, node(3)?, node(4)?];
+    let square = [Factor { dim: 1, deg: 1 }, Factor { dim: 1, deg: 1 }];
+    let zero = Bernstein::zero(&square)?;
+    let mut columns = [[zero.clone(), zero.clone(), zero.clone()], [zero.clone(), zero.clone(), zero.clone()], [zero.clone(), zero.clone(), zero]];
+    // A bilinear polynomial's coefficients are its values at the square's corners.
+    for (s, t) in [(0usize, 0usize), (1, 0), (1, 1), (0, 1)] {
+        let (sf, tf) = (s as f64, t as f64);
+        let base = |k: usize| p0[k] * (1.0 - sf) * (1.0 - tf) + p1[k] * sf * (1.0 - tf) + p2[k] * sf * tf + p3[k] * (1.0 - sf) * tf;
+        let b_s = |k: usize| (p1[k] - p0[k]) * (1.0 - tf) + (p2[k] - p3[k]) * tf;
+        let b_t = |k: usize| (p3[k] - p0[k]) * (1.0 - sf) + (p2[k] - p1[k]) * sf;
+        for (k, top) in apex.iter().enumerate() {
+            let values = [0.5 * b_s(k), 0.5 * b_t(k), top - base(k) + (sf - 0.5) * b_s(k) + (tf - 0.5) * b_t(k)];
+            for (column, v) in columns.iter_mut().zip(values) {
+                if let Some(c) = column.get_mut(k) {
+                    c.set([s, t, 0], v);
+                }
+            }
+        }
+    }
+    Ok(columns)
+}
+
 /// The Jacobian of a solid element whose node coordinates are already checked.
 pub(crate) fn jacobian(kind: ElementKind, xyz: &[[f64; 3]]) -> Result<Jacobian> {
-    let (_, dirs) = domain(kind)?;
-    let g = geometry(kind, xyz)?;
-    let column = |(f, k): (usize, usize)| -> Result<[Bernstein; 3]> { Ok([g[0].derivative(f, k)?, g[1].derivative(f, k)?, g[2].derivative(f, k)?]) };
-    let columns = [column(dirs[0])?, column(dirs[1])?, column(dirs[2])?];
+    let columns = if kind == ElementKind::Pyramid5 {
+        pyramid_columns(xyz)?
+    } else {
+        let (_, dirs) = domain(kind)?;
+        let g = geometry(kind, xyz)?;
+        let column =
+            |(f, k): (usize, usize)| -> Result<[Bernstein; 3]> { Ok([g[0].derivative(f, k)?, g[1].derivative(f, k)?, g[2].derivative(f, k)?]) };
+        [column(dirs[0])?, column(dirs[1])?, column(dirs[2])?]
+    };
     let [a, b, c] = &columns;
     // det [a b c] = a · (b × c)
     let cross = |i: usize, j: usize| -> Result<Bernstein> { b[i].mul(&c[j])?.axpy(-1.0, &b[j].mul(&c[i])?) };

crates/mesh/src/quality/mod.rs

@@ -197,7 +197,13 @@ fn checked(kind: ElementKind, xyz: &[[f64; 3]], p: &[[f64; 3]]) -> Result<(Optio
     }
     let jacobian::Jacobian { columns, det } = jacobian::jacobian(kind, p)?;
     let floor = jacobian::admissible_floor(&columns);
-    Ok((Some(columns), jacobian::certify(det, floor)?))
+    let mut check = jacobian::certify(det, floor)?;
+    if kind == ElementKind::Pyramid5 {
+        // The determinant lives on the base's unit square; the pyramid's reference volume is 4/3
+        // of it (jacobian.rs).
+        check.volume *= 4.0 / 3.0;
+    }
+    Ok((Some(columns), check))
 }
 
 /// The exact orientation of a tetrahedron (positive for node 4 on the side of face 1-2-3's
@@ -249,6 +255,11 @@ pub fn element_quality(kind: ElementKind, xyz: &[[f64; 3]]) -> Result<ElementQua
             let edges = [(0, 1), (1, 2), (2, 0), (3, 4), (4, 5), (5, 3), (0, 3), (1, 4), (2, 5)];
             vec![(Metric::ScaledJacobian, shape::wedge_scaled_jacobian(c)), (Metric::EdgeRatio, shape::edge_ratio(&c, &edges))]
         }
+        ElementKind::Pyramid5 => {
+            let c: [[f64; 3]; 5] = corners(&p)?;
+            let edges = [(0, 1), (1, 2), (2, 3), (3, 0), (0, 4), (1, 4), (2, 4), (3, 4)];
+            vec![(Metric::ScaledJacobian, shape::pyramid_scaled_jacobian(c)), (Metric::EdgeRatio, shape::edge_ratio(&c, &edges))]
+        }
         ElementKind::Hex8 | ElementKind::Hex20 => vec![(Metric::EdgeRatio, shape::edge_ratio(&corners::<8>(&p)?, &jacobian::HEX_EDGES))],
     };
     let (columns, check) = checked(kind, xyz, &p)?;

crates/mesh/src/quality/shape.rs

@@ -184,6 +184,17 @@ pub(crate) fn wedge_scaled_jacobian(p: [P; 6]) -> f64 {
         .fold(f64::INFINITY, |m, &x| m.min(x))
 }
 
+/// Pyramid scaled Jacobian: the smallest determinant of the unit edge vectors at the four base
+/// corners (to the next and previous base corners and to the apex, right-handed), divided by its
+/// value √2/2 on the pyramid whose edges are all equal, so that pyramid scores 1.
+pub(crate) fn pyramid_scaled_jacobian(p: [P; 5]) -> f64 {
+    let worst = [(0, 1, 3), (1, 2, 0), (2, 3, 1), (3, 0, 2)]
+        .map(|(o, a, b)| corner_det(sub(p[a], p[o]), sub(p[b], p[o]), sub(p[4], p[o])))
+        .iter()
+        .fold(f64::INFINITY, |m, &x| m.min(x));
+    (worst / std::f64::consts::FRAC_1_SQRT_2).clamp(-1.0, 1.0)
+}
+
 /// Scaled Jacobian and Knupp shape of one Jacobian matrix `[a b c]` (hexahedron corners and centre).
 pub(crate) fn jacobian_metrics(a: P, b: P, c: P) -> (f64, f64) {
     let d = det(a, b, c);

crates/mesh/src/quality/tests.rs

@@ -206,6 +206,34 @@ fn wedges() {
     assert_eq!(q.validity(), Some(Validity::Invalid));
 }
 
+#[test]
+fn pyramids() {
+    let h = 0.5_f64.sqrt();
+    let regular = [[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0.5, 0.5, h]];
+    let q = element_quality(ElementKind::Pyramid5, &regular).unwrap();
+    assert!(close(q.size, h / 3.0), "{}", q.size);
+    assert_eq!(q.validity(), Some(Validity::Valid));
+    assert!(close(metric(&q, Metric::ScaledJacobian), 1.0));
+    assert!(close(metric(&q, Metric::EdgeRatio), 1.0));
+
+    // A sheared apex keeps the volume (base × height / 3) and stays valid while it is over the base.
+    let sheared = [[0., 0., 0.], [2., 0., 0.], [2., 1., 0.], [0., 1., 0.], [1.7, 0.2, 1.0]];
+    let q = element_quality(ElementKind::Pyramid5, &sheared).unwrap();
+    assert!(close(q.size, 2.0 / 3.0), "{}", q.size);
+    assert_eq!(q.validity(), Some(Validity::Valid));
+
+    // Over a planar base the map is valid wherever the apex is above the base's plane, even past
+    // an edge (an oblique pyramid). A bow-tie base folds it; an apex below turns it inside out.
+    let oblique = [[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [1.6, 0.5, 0.3]];
+    assert_eq!(element_quality(ElementKind::Pyramid5, &oblique).unwrap().validity(), Some(Validity::Valid));
+    let bow_tie = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [1., 1., 0.], [0.5, 0.5, h]];
+    assert_eq!(element_quality(ElementKind::Pyramid5, &bow_tie).unwrap().validity(), Some(Validity::Invalid));
+    let below: Vec<[f64; 3]> = regular.iter().map(|p| [p[0], p[1], -p[2]]).collect();
+    let q = element_quality(ElementKind::Pyramid5, &below).unwrap();
+    assert!(close(q.size, -h / 3.0));
+    assert_eq!(q.validity(), Some(Validity::Invalid));
+}
+
 #[test]
 fn shells() {
     let h = 3.0_f64.sqrt() / 2.0;