From 3a5a8213a11bcdd132660f7af6813bbb1e006515 Mon Sep 17 00:00:00 2001
From: jyviko <2501648+jyviko@users.noreply.github.com>
Date: Fri, 16 Oct 2026 12:46:30 -0400
Subject: [PATCH] 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/decks/src/abaqus.rs          |  4 +-
 crates/decks/src/lsdyna/elements.rs |  8 +++-
 crates/decks/src/nastran/mod.rs     |  3 +-
 crates/decks/tests/abaqus.rs        |  3 +-
 crates/engine/src/cmd/forms.rs      |  1 +
 crates/engine/src/cmd/mesh.rs       |  2 +-
 crates/engine/src/cmd/mirror.rs     |  2 +-
 crates/engine/src/view.rs           |  3 ++
 crates/fem/src/basis.rs             |  1 +
 crates/fem/src/element.rs           | 17 ++++++++-
 crates/fem/src/finite.rs            | 12 ++++--
 crates/fem/src/quadrature.rs        | 32 ++++++++++++++++
 crates/fem/src/shape.rs             |  7 ++++
 crates/fem/src/tests.rs             |  5 ++-
 crates/io/src/mesh.rs               |  3 +-
 crates/io/src/vtk.rs                |  9 +++--
 crates/mesh/src/mirror.rs           |  2 +
 crates/mesh/src/quality/jacobian.rs | 46 ++++++++++++++++++++--
 crates/mesh/src/quality/mod.rs      | 13 ++++++-
 crates/mesh/src/quality/shape.rs    | 11 ++++++
 crates/mesh/src/quality/tests.rs    | 28 ++++++++++++++
 crates/model/src/mesh.rs            | 19 ++++++++--
 crates/solver/tests/common/mod.rs   | 59 +++++++++++++++++++++++++++++
 crates/solver/tests/verification.rs | 50 ++++++++++++++++++++++++
 crates/ui-egui/src/theme.rs         |  5 +++
 docs/elements.md                    |  3 +-
 docs/fem.md                         |  3 +-
 docs/formats/fieldcraft.md          |  2 +-
 28 files changed, 325 insertions(+), 28 deletions(-)

diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs
index 0a3893b..aacd131 100644
--- a/crates/decks/src/abaqus.rs
+++ b/crates/decks/src/abaqus.rs
@@ -85,6 +85,7 @@
 //! | tet4 | C3D4 (H) | | | |
 //! | tet10 | C3D10 (H) | | | |
 //! | wedge6 | C3D6 (H) | | | |
+//! | pyramid5 | C3D5 (H) | | | |
 //! | hex8 | C3D8, C3D8R, C3D8I (H) | | | |
 //! | hex20 | C3D20, C3D20R (H) | | | |
 //! | tri3, quad4 | | S3 / S3R, S4 / S4R | | |
@@ -231,11 +232,12 @@ pub fn element_type(kind: ElementKind, section: &Section) -> std::result::Result
         (S::Solid {}, _, true, true) => return Err("incompatible modes with reduced integration".into()),
         (S::Solid {}, _, _, true) => return Err(format!("incompatible modes on {} elements (hex8 only)", kind.name())),
         (S::Solid {}, K::Tet4, ..) => "C3D4",
-        (S::Solid {}, K::Tet10 | K::Wedge6, true, _) => {
+        (S::Solid {}, K::Tet10 | K::Wedge6 | K::Pyramid5, true, _) => {
             return Err(format!("reduced integration of {} elements (Abaqus has no reduced-integration variant)", kind.name()));
         }
         (S::Solid {}, K::Tet10, ..) => "C3D10",
         (S::Solid {}, K::Wedge6, ..) => "C3D6",
+        (S::Solid {}, K::Pyramid5, ..) => "C3D5",
         (S::Solid {}, K::Hex8, r, _) => ["C3D8", "C3D8R"][usize::from(r)],
         (S::Solid {}, K::Hex20, r, _) => ["C3D20", "C3D20R"][usize::from(r)],
         (S::Shell { .. }, K::Tri3, r, _) => ["S3", "S3R"][usize::from(r)],
diff --git a/crates/decks/src/lsdyna/elements.rs b/crates/decks/src/lsdyna/elements.rs
index 7795ebf..84d984c 100644
--- a/crates/decks/src/lsdyna/elements.rs
+++ b/crates/decks/src/lsdyna/elements.rs
@@ -36,7 +36,7 @@ impl Keyword {
     pub(crate) fn of(kind: ElementKind) -> Keyword {
         match kind {
             ElementKind::Hex20 => Keyword::SolidH20,
-            ElementKind::Tet4 | ElementKind::Tet10 | ElementKind::Wedge6 | ElementKind::Hex8 => Keyword::Solid,
+            ElementKind::Tet4 | ElementKind::Tet10 | ElementKind::Wedge6 | ElementKind::Pyramid5 | ElementKind::Hex8 => Keyword::Solid,
             ElementKind::Tri3 | ElementKind::Tri6 | ElementKind::Quad4 | ElementKind::Quad8 => Keyword::Shell,
             ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => Keyword::Beam,
         }
@@ -64,6 +64,8 @@ pub(crate) fn order(kind: ElementKind) -> &'static [usize] {
         ElementKind::Tet4 => &[0, 1, 2, 3, 3, 3, 3, 3],
         ElementKind::Tet10 => &[0, 1, 2, 3, 4, 5, 6, 8, 9, 7],
         ElementKind::Wedge6 => &[1, 0, 3, 4, 2, 2, 5, 5],
+        // The degenerate hexahedron: the base, then the apex four times.
+        ElementKind::Pyramid5 => &[0, 1, 2, 3, 4, 4, 4, 4],
         ElementKind::Hex8 => &[0, 1, 2, 3, 4, 5, 6, 7],
         ElementKind::Hex20 => &[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19],
     }
@@ -95,7 +97,7 @@ pub(crate) struct Elform {
 ///
 /// Solids: hex8 2 (eight-point), reduced 1 (one-point) with hourglass control, incompatible
 /// modes −18 (eight-point enhanced strain with incompatible modes); wedges 2 (as degenerate
-/// hexahedra), reduced 15 (two-point pentahedron); tet4 10 (one point, exact for a linear
+/// hexahedra), reduced 15 (two-point pentahedron); pyramids 2 (as degenerate hexahedra); tet4 10 (one point, exact for a linear
 /// tetrahedron), hybrid 13 (nodal pressure); tet10 16; hex20 23. Shells: quad4 16 (fully
 /// integrated), reduced 2 (one-point) with hourglass control; tri3 4 (C0 triangle); quad8 23;
 /// tri6 24. Rods: 3 (truss). 2D (in the XY plane): plane stress 12, plane strain 13, both one-point
@@ -122,6 +124,8 @@ pub(crate) fn elform(kind: ElementKind, section: &SectionKind, f: Formulation, s
                 ElementKind::Hex8 => plain(2),
                 ElementKind::Wedge6 if reduced => plain(15),
                 ElementKind::Wedge6 => plain(2),
+                ElementKind::Pyramid5 if reduced => no("there is no reduced-integration form of"),
+                ElementKind::Pyramid5 => plain(2),
                 ElementKind::Tet4 => plain(10),
                 _ if reduced && kind.is_quadratic() => no("there is no reduced-integration form of"),
                 ElementKind::Tet10 => plain(16),
diff --git a/crates/decks/src/nastran/mod.rs b/crates/decks/src/nastran/mod.rs
index c910167..e2f4b36 100644
--- a/crates/decks/src/nastran/mod.rs
+++ b/crates/decks/src/nastran/mod.rs
@@ -6,7 +6,7 @@
 //! data (the entries). What the model holds becomes:
 //!
 //! - nodes and reference points: `GRID`;
-//! - solids: `CTETRA` (tet4, tet10), `CPENTA` (wedge6), `CHEXA` (hex8, hex20) with `PSOLID`;
+//! - solids: `CTETRA` (tet4, tet10), `CPENTA` (wedge6), `CPYRAM` (pyramid5), `CHEXA` (hex8, hex20) with `PSOLID`;
 //!   shells: `CTRIA3`, `CQUAD4`, `CTRIA6`, `CQUAD8` with `PSHELL`; beams: `CBEAM` with `PBEAML`;
 //!   rods: `CROD` with `PROD`. Nastran numbers a 20-node hexahedron's mid-side nodes bottom
 //!   edges, vertical edges, top edges, where the model has bottom, top, vertical
@@ -113,6 +113,7 @@ pub fn element_card(kind: ElementKind) -> &'static str {
         ElementKind::Quad8 => "CQUAD8",
         ElementKind::Tet4 | ElementKind::Tet10 => "CTETRA",
         ElementKind::Wedge6 => "CPENTA",
+        ElementKind::Pyramid5 => "CPYRAM",
         ElementKind::Hex8 | ElementKind::Hex20 => "CHEXA",
     }
 }
diff --git a/crates/decks/tests/abaqus.rs b/crates/decks/tests/abaqus.rs
index 74ebf37..0e1b87b 100644
--- a/crates/decks/tests/abaqus.rs
+++ b/crates/decks/tests/abaqus.rs
@@ -461,10 +461,11 @@ fn written_connectivity_has_a_positive_jacobian_in_abaqus_conventions() {
 #[test]
 fn model_faces_are_abaqus_faces() {
     // Each face's node positions (1-based, Abaqus element library, our transcription).
-    let table: [(ElementKind, &[&[usize]]); 3] = [
+    let table: [(ElementKind, &[&[usize]]); 4] = [
         (ElementKind::Hex8, &[&[1, 2, 3, 4], &[5, 8, 7, 6], &[1, 5, 6, 2], &[2, 6, 7, 3], &[3, 7, 8, 4], &[4, 8, 5, 1]]),
         (ElementKind::Tet4, &[&[1, 2, 3], &[1, 4, 2], &[2, 4, 3], &[3, 4, 1]]),
         (ElementKind::Wedge6, &[&[1, 2, 3], &[4, 6, 5], &[1, 4, 5, 2], &[2, 5, 6, 3], &[3, 6, 4, 1]]),
+        (ElementKind::Pyramid5, &[&[1, 2, 3, 4], &[1, 5, 2], &[2, 5, 3], &[3, 5, 4], &[4, 5, 1]]),
     ];
     for (kind, faces) in table {
         assert_eq!(kind.faces().len(), faces.len());
diff --git a/crates/engine/src/cmd/forms.rs b/crates/engine/src/cmd/forms.rs
index 142f1ae..5abbb28 100644
--- a/crates/engine/src/cmd/forms.rs
+++ b/crates/engine/src/cmd/forms.rs
@@ -1384,6 +1384,7 @@ pub const ELEMENT_ADD: &[Field] = &[
             ("tet4", "tet4"),
             ("tet10", "tet10"),
             ("wedge6", "wedge6"),
+            ("pyramid5", "pyramid5"),
             ("hex8", "hex8"),
             ("hex20", "hex20"),
         ],
diff --git a/crates/engine/src/cmd/mesh.rs b/crates/engine/src/cmd/mesh.rs
index 659434e..ef2adfa 100644
--- a/crates/engine/src/cmd/mesh.rs
+++ b/crates/engine/src/cmd/mesh.rs
@@ -40,7 +40,7 @@ pub const NODE_REMOVE: CommandSpec = CommandSpec {
 pub const ELEMENT_ADD: CommandSpec = CommandSpec {
     id: "mesh.element.add",
     label: "Add Element",
-    params: r#"{"label": 1, "kind": "rod2|beam2|tri3|tri6|quad4|quad8|tet4|tet10|wedge6|hex8|hex20", "nodes": [labels]} — adds or replaces"#,
+    params: r#"{"label": 1, "kind": "rod2|beam2|tri3|tri6|quad4|quad8|tet4|tet10|wedge6|pyramid5|hex8|hex20", "nodes": [labels]} — adds or replaces"#,
     form: Some(super::forms::ELEMENT_ADD),
     kind: CommandKind::Model,
     requires: &[Requirement::Nodes],
diff --git a/crates/engine/src/cmd/mirror.rs b/crates/engine/src/cmd/mirror.rs
index 67568ed..0cc5c2e 100644
--- a/crates/engine/src/cmd/mirror.rs
+++ b/crates/engine/src/cmd/mirror.rs
@@ -23,7 +23,7 @@ use crate::{CommandKind, CommandSpec, EngineError, Result, Session};
 pub const MIRROR: CommandSpec = CommandSpec {
     id: "mesh.mirror",
     label: "Mirror Mesh",
-    params: r#"{"plane": "Sym" (a datum plane) or {"point": [x, y, z], "normal": [nx, ny, nz]}, "bodies": ["Part/a"] (default: every body the mesh is tagged with), "as": {"Part/a": "Part/b"} (optional: the body each image is; default the design's body at the reflected place, else "<body>/mirror")} — replaces each image body's mesh (as mesh.tet.generate) with the reflection of the bodies' elements (tet4, tet10, wedge6, hex8, hex20) across the plane: each copy right-handed (its nodes in reflected order), the copies of the nodes on the plane merged with them so the mesh is conforming there, geometry tags carried to the image (its body, and on it the faces the reflected faces lie on, so sets by face name select both halves). Reports the elements and nodes added, the nodes shared on the plane and each image's body"#,
+    params: r#"{"plane": "Sym" (a datum plane) or {"point": [x, y, z], "normal": [nx, ny, nz]}, "bodies": ["Part/a"] (default: every body the mesh is tagged with), "as": {"Part/a": "Part/b"} (optional: the body each image is; default the design's body at the reflected place, else "<body>/mirror")} — replaces each image body's mesh (as mesh.tet.generate) with the reflection of the bodies' elements (tet4, tet10, wedge6, pyramid5, hex8, hex20) across the plane: each copy right-handed (its nodes in reflected order), the copies of the nodes on the plane merged with them so the mesh is conforming there, geometry tags carried to the image (its body, and on it the faces the reflected faces lie on, so sets by face name select both halves). Reports the elements and nodes added, the nodes shared on the plane and each image's body"#,
     form: Some(super::forms::MIRROR),
     kind: CommandKind::Model,
     requires: &[Requirement::Mesh],
diff --git a/crates/engine/src/view.rs b/crates/engine/src/view.rs
index bb41b15..24d711f 100644
--- a/crates/engine/src/view.rs
+++ b/crates/engine/src/view.rs
@@ -28,6 +28,7 @@ pub struct Palette {
     pub hex: Rgb,
     pub tet: Rgb,
     pub wedge: Rgb,
+    pub pyramid: Rgb,
     pub shell: Rgb,
     pub line_element: Rgb,
     pub mesh_edge: Rgb,
@@ -61,6 +62,7 @@ impl Palette {
         hex: Rgb(0x78, 0x9c, 0xc0),
         tet: Rgb(0x80, 0xb4, 0x8a),
         wedge: Rgb(0xc0, 0xaa, 0x78),
+        pyramid: Rgb(0xa8, 0x90, 0xc8),
         shell: Rgb(0xc0, 0x8c, 0x78),
         line_element: Rgb(0xc8, 0x8c, 0x1e),
         mesh_edge: Rgb(0x1e, 0x23, 0x2a),
@@ -100,6 +102,7 @@ impl Palette {
             ElementKind::Hex8 | ElementKind::Hex20 => self.hex,
             ElementKind::Tet4 | ElementKind::Tet10 => self.tet,
             ElementKind::Wedge6 => self.wedge,
+            ElementKind::Pyramid5 => self.pyramid,
             ElementKind::Tri3 | ElementKind::Tri6 | ElementKind::Quad4 | ElementKind::Quad8 => self.shell,
             ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => self.line_element,
         }
diff --git a/crates/fem/src/basis.rs b/crates/fem/src/basis.rs
index f9943f4..2c53480 100644
--- a/crates/fem/src/basis.rs
+++ b/crates/fem/src/basis.rs
@@ -207,6 +207,7 @@ impl Basis {
             ElementKind::Tet4 => Basis::Tet4,
             ElementKind::Tet10 => Basis::Tet10,
             ElementKind::Wedge6 => Basis::Wedge6,
+            ElementKind::Pyramid5 => Basis::Pyramid5,
             ElementKind::Hex8 => Basis::Hex8,
             ElementKind::Hex20 => Basis::Hex20,
         }
diff --git a/crates/fem/src/element.rs b/crates/fem/src/element.rs
index a34a819..b32fac1 100644
--- a/crates/fem/src/element.rs
+++ b/crates/fem/src/element.rs
@@ -47,6 +47,7 @@ use crate::{Error, Result};
 pub enum Family {
     Tetrahedron,
     Wedge,
+    Pyramid,
     Hexahedron,
     /// The 2D elements of plane and axisymmetric models ([`crate::plane`]).
     Triangle,
@@ -74,7 +75,7 @@ pub struct ElementType {
 }
 
 /// The element table.
-pub static ELEMENT_TYPES: [ElementType; 5] = [
+pub static ELEMENT_TYPES: [ElementType; 6] = [
     ElementType {
         kind: ElementKind::Tet4,
         family: Family::Tetrahedron,
@@ -108,6 +109,17 @@ pub static ELEMENT_TYPES: [ElementType; 5] = [
         mass: Rule::Wedge,
         incompatible_modes: false,
     },
+    ElementType {
+        kind: ElementKind::Pyramid5,
+        family: Family::Pyramid,
+        nodes: 5,
+        corners: 5,
+        shape: shape::pyramid5,
+        full: Rule::Pyramid(2),
+        reduced: None,
+        mass: Rule::Pyramid(3),
+        incompatible_modes: false,
+    },
     ElementType {
         kind: ElementKind::Hex8,
         family: Family::Hexahedron,
@@ -137,7 +149,7 @@ pub fn element_type(kind: ElementKind) -> Result<&'static ElementType> {
     ELEMENT_TYPES
         .iter()
         .find(|t| t.kind == kind)
-        .ok_or_else(|| Error::NotSupported(format!("{} elements (the solver has tet4, tet10, wedge6, hex8 and hex20 so far)", kind.name())))
+        .ok_or_else(|| Error::NotSupported(format!("{} elements (the solver has tet4, tet10, wedge6, pyramid5, hex8 and hex20 so far)", kind.name())))
 }
 
 /// How a solid element is formulated.
@@ -730,6 +742,7 @@ pub fn extrapolation(ty: &ElementType, rule: Rule) -> Result<Extrapolation> {
     let corner_basis = |xi: [f64; 3]| match ty.family {
         Family::Tetrahedron => shape::tet4(xi),
         Family::Wedge => shape::wedge6(xi),
+        Family::Pyramid => shape::pyramid5(xi),
         Family::Hexahedron => shape::hex8(xi),
         Family::Triangle => Basis::Tri3.eval(xi),
         Family::Quadrilateral => Basis::Quad4.eval(xi),
diff --git a/crates/fem/src/finite.rs b/crates/fem/src/finite.rs
index e340de2..7a2d272 100644
--- a/crates/fem/src/finite.rs
+++ b/crates/fem/src/finite.rs
@@ -833,9 +833,15 @@ mod tests {
         let x8: Vec<[f64; 3]> = x20[..8].to_vec();
         let x6 = vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.1], [0., 0., 1.], [1.1, 0., 1.], [0., 1., 1.2]];
         let x4 = vec![x10[0], x10[1], x10[2], x10[3]];
-        for (kind, x) in
-            [(ElementKind::Tet4, &x4), (ElementKind::Tet10, &x10), (ElementKind::Wedge6, &x6), (ElementKind::Hex8, &x8), (ElementKind::Hex20, &x20)]
-        {
+        let x5 = vec![[0., 0., 0.], [1., 0., 0.1], [1.1, 1., 0.], [0., 1., 0.], [0.5, 0.4, 0.9]];
+        for (kind, x) in [
+            (ElementKind::Tet4, &x4),
+            (ElementKind::Tet10, &x10),
+            (ElementKind::Wedge6, &x6),
+            (ElementKind::Pyramid5, &x5),
+            (ElementKind::Hex8, &x8),
+            (ElementKind::Hex20, &x20),
+        ] {
             for m in [&svk, &neo] {
                 let solid = Solid::new(kind, x, Formulation::Full).unwrap();
                 check_tangent(&solid, m, &displacement(x.len(), 0.08));
diff --git a/crates/fem/src/quadrature.rs b/crates/fem/src/quadrature.rs
index ff863fa..6d59917 100644
--- a/crates/fem/src/quadrature.rs
+++ b/crates/fem/src/quadrature.rs
@@ -66,6 +66,10 @@ pub enum Rule {
     /// Three triangle points times two Gauss points in ζ on the reference wedge: point k of each
     /// layer sits near corner k, the lower layer first.
     Wedge,
+    /// The collapsed (Duffy) product on the reference pyramid, n points per direction: Gauss–Legendre
+    /// in u, v ∈ [−1, 1] and w ∈ [0, 1], mapped by ξ = u(1 − w), η = v(1 − w), ζ = w with the
+    /// weight times (1 − w)², which makes the pyramid's rational functions polynomials in (u, v, w).
+    Pyramid(u8),
 }
 
 impl Rule {
@@ -125,6 +129,19 @@ impl Rule {
                 let tri = [[b, b], [a, b], [b, a]];
                 gauss_legendre(2).iter().flat_map(|&(z, wz)| tri.iter().map(move |t| Point { xi: [t[0], t[1], z], weight: wz / 6.0 })).collect()
             }
+            Rule::Pyramid(n) => {
+                let g = gauss_legendre(usize::from(n));
+                let mut out = Vec::with_capacity(g.len().pow(3));
+                for &(z, wz) in g {
+                    let (w, ww) = ((z + 1.0) / 2.0, wz / 2.0);
+                    for &(v, wv) in g {
+                        for &(u, wu) in g {
+                            out.push(Point { xi: [u * (1.0 - w), v * (1.0 - w), w], weight: wu * wv * ww * (1.0 - w) * (1.0 - w) });
+                        }
+                    }
+                }
+                out
+            }
         }
     }
 }
@@ -182,6 +199,21 @@ mod tests {
                 }
             }
         }
+        // Pyramid: ∫ ξ^a η^b ζ^c over the pyramid of base [−1, 1]² at ζ = 0 and apex (0, 0, 1)
+        // = ∫_0^1 ζ^c ∫ ξ^a ∫ η^b over the square of half-width 1 − ζ: for even a, b,
+        // 4/((a+1)(b+1)) ∫_0^1 ζ^c (1 − ζ)^(a+b+2) = 4/((a+1)(b+1)) · B(c + 1, a + b + 3).
+        let beta = |p: i32, q: i32| f(p - 1) * f(q - 1) / f(p + q - 1);
+        // Three points per direction are exact to degree 5 in w, i.e. a + b + c ≤ 3.
+        for a in 0..=3 {
+            for b in 0..=3 - a {
+                for c in 0..=3 - a - b {
+                    let exact = if a % 2 == 1 || b % 2 == 1 { 0.0 } else { 4.0 / f64::from((a + 1) * (b + 1)) * beta(c + 1, a + b + 3) };
+                    let got = integrate(Rule::Pyramid(3), a, b, c);
+                    assert!((got - exact).abs() < 1e-14, "pyramid ξ^{a} η^{b} ζ^{c}: {got} vs {exact}");
+                }
+            }
+        }
+        assert!((integrate(Rule::Pyramid(2), 0, 0, 0) - 4.0 / 3.0).abs() < 1e-15);
         // Triangle: ∫ ξ^a η^b = a! b! / (a + b + 2)!.
         for (rule, degree) in [(Rule::TriConical(3), 4), (Rule::TriConical(4), 6), (Rule::TriCentroid, 1), (Rule::Tri3, 2)] {
             for a in 0..=degree {
diff --git a/crates/fem/src/shape.rs b/crates/fem/src/shape.rs
index dfddb10..74c3a40 100644
--- a/crates/fem/src/shape.rs
+++ b/crates/fem/src/shape.rs
@@ -116,6 +116,13 @@ pub fn wedge6(xi: [f64; 3]) -> Shape {
     s
 }
 
+/// Five-node pyramid: base [−1, 1]² at ζ = 0, apex (0, 0, 1), the rational functions of
+/// [`crate::basis`] (Bedrosian 1992), conforming with tetrahedra on its triangles and with
+/// hexahedra on its base.
+pub fn pyramid5(xi: [f64; 3]) -> Shape {
+    crate::basis::Basis::Pyramid5.eval(xi)
+}
+
 /// Twenty-node serendipity hexahedron.
 pub fn hex20(xi: [f64; 3]) -> Shape {
     let mut s = Shape::new(20);
diff --git a/crates/fem/src/tests.rs b/crates/fem/src/tests.rs
index dd13b1d..d071a22 100644
--- a/crates/fem/src/tests.rs
+++ b/crates/fem/src/tests.rs
@@ -13,6 +13,7 @@ fn distorted(kind: ElementKind, curved: bool) -> Vec<[f64; 3]> {
     let corners: Vec<[f64; 3]> = match kind {
         ElementKind::Tet4 | ElementKind::Tet10 => vec![[0.1, -0.05, 0.0], [1.2, 0.1, 0.05], [0.2, 0.9, -0.1], [0.05, 0.15, 1.1]],
         ElementKind::Wedge6 => vec![[0.0, 0.0, 0.0], [2.1, 0.1, -0.1], [0.2, 1.3, 0.05], [0.1, -0.1, 0.9], [2.0, 0.05, 1.1], [0.15, 1.2, 1.0]],
+        ElementKind::Pyramid5 => vec![[0.0, 0.0, 0.0], [2.1, 0.1, -0.1], [2.3, 1.2, 0.05], [-0.1, 0.9, 0.1], [0.9, 0.6, 1.3]],
         _ => vec![
             [0.0, 0.0, 0.0],
             [2.1, 0.1, -0.1],
@@ -34,10 +35,11 @@ fn distorted(kind: ElementKind, curved: bool) -> Vec<[f64; 3]> {
     xyz
 }
 
-const CASES: [(ElementKind, Formulation); 10] = [
+const CASES: [(ElementKind, Formulation); 11] = [
     (ElementKind::Tet4, Formulation::Full),
     (ElementKind::Tet10, Formulation::Full),
     (ElementKind::Wedge6, Formulation::Full),
+    (ElementKind::Pyramid5, Formulation::Full),
     (ElementKind::Hex8, Formulation::Full),
     (ElementKind::Hex8, Formulation::IncompatibleModes),
     (ElementKind::Hex20, Formulation::Full),
@@ -156,6 +158,7 @@ fn extrapolation_reproduces_linear_fields_on_affine_elements() {
         let mut xyz: Vec<[f64; 3]> = match kind {
             ElementKind::Tet10 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]],
             ElementKind::Wedge6 => vec![[0., 0., 0.], [2., 0., 0.], [0.5, 1., 0.], [0., 0.2, 1.], [2., 0.2, 1.], [0.5, 1.2, 1.]],
+            ElementKind::Pyramid5 => vec![[0., 0., 0.], [2., 0., 0.], [2.5, 1., 0.], [0.5, 1., 0.], [1.2, 0.5, 1.]],
             _ => vec![[0., 0., 0.], [2., 0., 0.], [2.5, 1., 0.], [0.5, 1., 0.], [0., 0., 1.], [2., 0., 1.], [2.5, 1., 1.], [0.5, 1., 1.]],
         };
         let c = xyz.clone();
diff --git a/crates/io/src/mesh.rs b/crates/io/src/mesh.rs
index 602c69d..645f7ea 100644
--- a/crates/io/src/mesh.rs
+++ b/crates/io/src/mesh.rs
@@ -38,7 +38,7 @@ const PART: &str = "mesh.bin";
 const NO_TAG: u32 = u32::MAX;
 
 /// The element kind codes of the format.
-pub const KIND_CODES: [(ElementKind, u8); 12] = [
+pub const KIND_CODES: [(ElementKind, u8); 13] = [
     (ElementKind::Rod2, 1),
     (ElementKind::Beam2, 2),
     (ElementKind::Tri3, 3),
@@ -51,6 +51,7 @@ pub const KIND_CODES: [(ElementKind, u8); 12] = [
     (ElementKind::Hex8, 10),
     (ElementKind::Hex20, 11),
     (ElementKind::Rod3, 12),
+    (ElementKind::Pyramid5, 13),
 ];
 
 fn kind_code(kind: ElementKind) -> u8 {
diff --git a/crates/io/src/vtk.rs b/crates/io/src/vtk.rs
index eaf653a..f43b346 100644
--- a/crates/io/src/vtk.rs
+++ b/crates/io/src/vtk.rs
@@ -38,6 +38,7 @@ pub fn cell_type(kind: ElementKind) -> u8 {
         ElementKind::Tet4 => 10,                     // VTK_TETRA
         ElementKind::Tet10 => 24,                    // VTK_QUADRATIC_TETRA
         ElementKind::Wedge6 => 13,                   // VTK_WEDGE
+        ElementKind::Pyramid5 => 14,                 // VTK_PYRAMID
         ElementKind::Hex8 => 12,                     // VTK_HEXAHEDRON
         ElementKind::Hex20 => 25,                    // VTK_QUADRATIC_HEXAHEDRON
     }
@@ -347,6 +348,7 @@ mod tests {
                 [0., 0.5, 0.5],
             ],
             ElementKind::Wedge6 => &[[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [0., 1., 1.]],
+            ElementKind::Pyramid5 => &[[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0.5, 0.5, 1.]],
             ElementKind::Hex8 => &[[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [1., 1., 1.], [0., 1., 1.]],
             ElementKind::Hex20 => &[
                 [0., 0., 0.],
@@ -433,8 +435,8 @@ mod tests {
             assert_eq!(written, e.nodes, "{kind:?}");
         }
         let ids: Vec<u8> = ElementKind::ALL.iter().map(|&k| cell_type(k)).collect();
-        // rod2, rod3, beam2, tri3, tri6, quad4, quad8, tet4, tet10, wedge6, hex8, hex20
-        assert_eq!(ids, [3, 21, 3, 5, 22, 9, 23, 10, 24, 13, 12, 25]);
+        // rod2, rod3, beam2, tri3, tri6, quad4, quad8, tet4, tet10, wedge6, pyramid5, hex8, hex20
+        assert_eq!(ids, [3, 21, 3, 5, 22, 9, 23, 10, 24, 13, 14, 12, 25]);
     }
 
     /// The points of the one cell as written, in VTK's order.
@@ -459,7 +461,7 @@ mod tests {
 
     /// VTK's convention for each 3D cell, from its parametric coordinates: the base face's
     /// right-hand normal points to the rest of the cell (tetrahedron 0-1-2 to 3, hexahedron
-    /// 0-1-2-3 to 4-7, wedge 0-1-2 to 3-5). The signed measure is positive when the written cell
+    /// 0-1-2-3 to 4-7, wedge 0-1-2 to 3-5, pyramid 0-1-2-3 to 4). The signed measure is positive when the written cell
     /// follows it; VTK 9.7's vtkCellSizeFilter and vtkMeshQuality agree (a unit wedge written
     /// this way measures +0.5, reordered 0-2-1-3-5-4 it measures -0.5).
     #[test]
@@ -470,6 +472,7 @@ mod tests {
             (ElementKind::Hex8, 4, 4..8),
             (ElementKind::Hex20, 4, 4..8),
             (ElementKind::Wedge6, 3, 3..6),
+            (ElementKind::Pyramid5, 4, 4..5),
         ] {
             let p = cell_points(&write_vtu(&one(kind), &[], &[], Encoding::Ascii).unwrap());
             let b = &p[..base];
diff --git a/crates/mesh/src/mirror.rs b/crates/mesh/src/mirror.rs
index 52f5e8a..1d12526 100644
--- a/crates/mesh/src/mirror.rs
+++ b/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],
diff --git a/crates/mesh/src/quality/jacobian.rs b/crates/mesh/src/quality/jacobian.rs
index 815211a..2ea9dba 100644
--- a/crates/mesh/src/quality/jacobian.rs
+++ b/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])?) };
diff --git a/crates/mesh/src/quality/mod.rs b/crates/mesh/src/quality/mod.rs
index d86c49c..addc93b 100644
--- a/crates/mesh/src/quality/mod.rs
+++ b/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)?;
diff --git a/crates/mesh/src/quality/shape.rs b/crates/mesh/src/quality/shape.rs
index 5c92cc9..278492f 100644
--- a/crates/mesh/src/quality/shape.rs
+++ b/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);
diff --git a/crates/mesh/src/quality/tests.rs b/crates/mesh/src/quality/tests.rs
index 0bd50b8..113444a 100644
--- a/crates/mesh/src/quality/tests.rs
+++ b/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;
diff --git a/crates/model/src/mesh.rs b/crates/model/src/mesh.rs
index fa0b576..c708e96 100644
--- a/crates/model/src/mesh.rs
+++ b/crates/model/src/mesh.rs
@@ -17,6 +17,7 @@ pub struct Node {
 /// - rod3: the two ends, then the middle node;
 /// - tet10: mid-side nodes 5–10 on edges 1-2, 2-3, 3-1, 1-4, 2-4, 3-4;
 /// - wedge6: triangle 1-2-3, then 4-5-6 on the side its right-hand normal points to;
+/// - pyramid5: base 1-2-3-4, then the apex 5 on the side its right-hand normal points to;
 /// - hex8/hex20: face 1-2-3-4, then 5-6-7-8 above it; hex20's mid-side nodes 9–12 lie on the
 ///   edges of 1-2-3-4, 13–16 on those of 5-6-7-8, and 17–20 on edges 1-5, 2-6, 3-7, 4-8.
 ///
@@ -47,6 +48,9 @@ pub enum ElementKind {
     Tet10,
     /// Six-node wedge (pentahedron).
     Wedge6,
+    /// Five-node pyramid: a quadrilateral base and an apex; the transition between a hexahedron's
+    /// face and tetrahedra.
+    Pyramid5,
     /// Eight-node hexahedron.
     Hex8,
     /// Twenty-node hexahedron.
@@ -54,7 +58,7 @@ pub enum ElementKind {
 }
 
 impl ElementKind {
-    pub const ALL: [ElementKind; 12] = [
+    pub const ALL: [ElementKind; 13] = [
         ElementKind::Rod2,
         ElementKind::Rod3,
         ElementKind::Beam2,
@@ -65,6 +69,7 @@ impl ElementKind {
         ElementKind::Tet4,
         ElementKind::Tet10,
         ElementKind::Wedge6,
+        ElementKind::Pyramid5,
         ElementKind::Hex8,
         ElementKind::Hex20,
     ];
@@ -74,6 +79,7 @@ impl ElementKind {
             ElementKind::Rod2 | ElementKind::Beam2 => 2,
             ElementKind::Rod3 | ElementKind::Tri3 => 3,
             ElementKind::Quad4 | ElementKind::Tet4 => 4,
+            ElementKind::Pyramid5 => 5,
             ElementKind::Tri6 | ElementKind::Wedge6 => 6,
             ElementKind::Quad8 | ElementKind::Hex8 => 8,
             ElementKind::Tet10 => 10,
@@ -87,6 +93,7 @@ impl ElementKind {
             ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => 2,
             ElementKind::Tri3 | ElementKind::Tri6 => 3,
             ElementKind::Quad4 | ElementKind::Quad8 | ElementKind::Tet4 | ElementKind::Tet10 => 4,
+            ElementKind::Pyramid5 => 5,
             ElementKind::Wedge6 => 6,
             ElementKind::Hex8 | ElementKind::Hex20 => 8,
         }
@@ -102,7 +109,7 @@ impl ElementKind {
         match self {
             ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => 1,
             ElementKind::Tri3 | ElementKind::Tri6 | ElementKind::Quad4 | ElementKind::Quad8 => 2,
-            ElementKind::Tet4 | ElementKind::Tet10 | ElementKind::Wedge6 | ElementKind::Hex8 | ElementKind::Hex20 => 3,
+            ElementKind::Tet4 | ElementKind::Tet10 | ElementKind::Wedge6 | ElementKind::Pyramid5 | ElementKind::Hex8 | ElementKind::Hex20 => 3,
         }
     }
 
@@ -119,6 +126,7 @@ impl ElementKind {
             ElementKind::Tet4 => "tet4",
             ElementKind::Tet10 => "tet10",
             ElementKind::Wedge6 => "wedge6",
+            ElementKind::Pyramid5 => "pyramid5",
             ElementKind::Hex8 => "hex8",
             ElementKind::Hex20 => "hex20",
         }
@@ -135,7 +143,8 @@ impl ElementKind {
     ///   normal points out of the element. Hexahedra: 1 is the face of nodes 1-2-3-4, 2 the face
     ///   of 5-6-7-8, then the four sides along edges 1-2, 2-3, 3-4 and 4-1. Tetrahedra: 1 is
     ///   1-2-3, then the faces opposite node 3, node 1 and node 2. Wedges: the two triangles, then
-    ///   the sides along edges 1-2, 2-3 and 3-1.
+    ///   the sides along edges 1-2, 2-3 and 3-1. Pyramids: 1 is the base 1-2-3-4, then the
+    ///   triangles on edges 1-2, 2-3, 3-4 and 4-1 with the apex 5.
     /// - Shells: face 1 is the side the element's right-hand normal points to, face 2 the other.
     /// - Line elements have none.
     pub fn faces(self) -> &'static [&'static [usize]] {
@@ -148,6 +157,7 @@ impl ElementKind {
             ElementKind::Tet4 => &[&[0, 2, 1], &[0, 1, 3], &[1, 2, 3], &[0, 3, 2]],
             ElementKind::Tet10 => &[&[0, 2, 1, 6, 5, 4], &[0, 1, 3, 4, 8, 7], &[1, 2, 3, 5, 9, 8], &[0, 3, 2, 7, 9, 6]],
             ElementKind::Wedge6 => &[&[0, 2, 1], &[3, 4, 5], &[0, 1, 4, 3], &[1, 2, 5, 4], &[2, 0, 3, 5]],
+            ElementKind::Pyramid5 => &[&[0, 3, 2, 1], &[0, 1, 4], &[1, 2, 4], &[2, 3, 4], &[3, 0, 4]],
             ElementKind::Hex8 => &[&[0, 3, 2, 1], &[4, 5, 6, 7], &[0, 1, 5, 4], &[1, 2, 6, 5], &[2, 3, 7, 6], &[3, 0, 4, 7]],
             ElementKind::Hex20 => &[
                 &[0, 3, 2, 1, 11, 10, 9, 8],
@@ -174,6 +184,7 @@ impl ElementKind {
             ElementKind::Tet4 => &[&[0, 1], &[1, 2], &[2, 0], &[0, 3], &[1, 3], &[2, 3]],
             ElementKind::Tet10 => &[&[0, 1, 4], &[1, 2, 5], &[2, 0, 6], &[0, 3, 7], &[1, 3, 8], &[2, 3, 9]],
             ElementKind::Wedge6 => &[&[0, 1], &[1, 2], &[2, 0], &[3, 4], &[4, 5], &[5, 3], &[0, 3], &[1, 4], &[2, 5]],
+            ElementKind::Pyramid5 => &[&[0, 1], &[1, 2], &[2, 3], &[3, 0], &[0, 4], &[1, 4], &[2, 4], &[3, 4]],
             ElementKind::Hex8 => &[&[0, 1], &[1, 2], &[2, 3], &[3, 0], &[4, 5], &[5, 6], &[6, 7], &[7, 4], &[0, 4], &[1, 5], &[2, 6], &[3, 7]],
             ElementKind::Hex20 => &[
                 &[0, 1, 8],
@@ -539,6 +550,7 @@ mod tests {
             ElementKind::Quad4 | ElementKind::Quad8 => vec![[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.]],
             ElementKind::Tet4 | ElementKind::Tet10 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]],
             ElementKind::Wedge6 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [0., 1., 1.]],
+            ElementKind::Pyramid5 => vec![[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0.5, 0.5, 1.]],
             ElementKind::Hex8 | ElementKind::Hex20 => {
                 vec![[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [1., 1., 1.], [0., 1., 1.]]
             }
@@ -589,6 +601,7 @@ mod tests {
         assert_eq!(ElementKind::Hex8.faces().len(), 6);
         assert_eq!(ElementKind::Tet4.faces().len(), 4);
         assert_eq!(ElementKind::Wedge6.faces().len(), 5);
+        assert_eq!(ElementKind::Pyramid5.faces().len(), 5);
         assert_eq!(ElementKind::Hex20.edges().len(), 12);
         for kind in ElementKind::ALL {
             let mid_sides = kind.edges().iter().filter(|e| e.len() == 3).count();
diff --git a/crates/solver/tests/common/mod.rs b/crates/solver/tests/common/mod.rs
index 6715034..045c67f 100644
--- a/crates/solver/tests/common/mod.rs
+++ b/crates/solver/tests/common/mod.rs
@@ -97,6 +97,65 @@ pub fn grid(n: [usize; 3], kind: ElementKind, map: impl Fn([f64; 3]) -> [f64; 3]
     Grid { model, param, kind }
 }
 
+/// A hex–tet transition mesh on [0, 1]³ with `n` cells per direction: the cells with x below
+/// one half are hex8, every other cell is split from its centre into six pyramids over its faces,
+/// and each pyramid whose face is not shared with a hex cell into two tetrahedra (split along
+/// the diagonal through the face's lowest node label, so neighbours agree). The mesh is
+/// conforming: hexahedra meet pyramids on whole quadrilaterals, everything else on triangles.
+pub fn transition_grid(n: [usize; 3], map: impl Fn([f64; 3]) -> [f64; 3]) -> Grid {
+    let size = n.map(|k| 2 * k + 1);
+    let label = |i: [usize; 3]| (1 + i[0] + size[0] * (i[1] + size[1] * i[2])) as u32;
+    let param_of = |i: [usize; 3]| [0, 1, 2].map(|k| i[k] as f64 / (2 * n[k]) as f64);
+    let mut model = Model::default();
+    let mut param = BTreeMap::new();
+    let mut node = |model: &mut Model, i: [usize; 3]| -> u32 {
+        let l = label(i);
+        if let std::collections::btree_map::Entry::Vacant(slot) = param.entry(l) {
+            let p = param_of(i);
+            model.set_node(l, map(p)).unwrap();
+            slot.insert(p);
+        }
+        l
+    };
+    let is_hex = |c: [isize; 3]| (0..3).all(|k| c[k] >= 0 && (c[k] as usize) < n[k]) && 2 * c[0] < n[0] as isize;
+    const CORNERS: [[usize; 3]; 8] = [[0, 0, 0], [2, 0, 0], [2, 2, 0], [0, 2, 0], [0, 0, 2], [2, 0, 2], [2, 2, 2], [0, 2, 2]];
+    // Hex8's faces (outward) and the cell across each.
+    let across: [[isize; 3]; 6] = [[0, 0, -1], [0, 0, 1], [0, -1, 0], [1, 0, 0], [0, 1, 0], [-1, 0, 0]];
+    let mut next = 1u32;
+    for cz in 0..n[2] {
+        for cy in 0..n[1] {
+            for cx in 0..n[0] {
+                let base = [2 * cx, 2 * cy, 2 * cz];
+                let at = |d: [usize; 3]| [base[0] + d[0], base[1] + d[1], base[2] + d[2]];
+                let cell = [cx as isize, cy as isize, cz as isize];
+                let corners: Vec<u32> = CORNERS.iter().map(|c| node(&mut model, at(*c))).collect();
+                if is_hex(cell) {
+                    model.set_element(next, ElementKind::Hex8, corners).unwrap();
+                    next += 1;
+                    continue;
+                }
+                let apex = node(&mut model, at([1, 1, 1]));
+                for (face, d) in ElementKind::Hex8.faces().iter().zip(across) {
+                    // Reversed, the outward face's normal points into the cell, at the apex.
+                    let q: Vec<u32> = face.iter().rev().map(|&i| corners[i]).collect();
+                    if is_hex([cell[0] + d[0], cell[1] + d[1], cell[2] + d[2]]) {
+                        model.set_element(next, ElementKind::Pyramid5, vec![q[0], q[1], q[2], q[3], apex]).unwrap();
+                        next += 1;
+                        continue;
+                    }
+                    let lowest = (0..4).min_by_key(|&i| q[i]).unwrap();
+                    let (a, b, c, e) = if lowest % 2 == 0 { (q[0], q[1], q[2], q[3]) } else { (q[1], q[2], q[3], q[0]) };
+                    for tet in [[a, b, c, apex], [a, c, e, apex]] {
+                        model.set_element(next, ElementKind::Tet4, tet.to_vec()).unwrap();
+                        next += 1;
+                    }
+                }
+            }
+        }
+    }
+    Grid { model, param, kind: ElementKind::Pyramid5 }
+}
+
 pub fn near(a: f64, b: f64) -> bool {
     (a - b).abs() < 1e-9
 }
diff --git a/crates/solver/tests/verification.rs b/crates/solver/tests/verification.rs
index 80add0f..bafe856 100644
--- a/crates/solver/tests/verification.rs
+++ b/crates/solver/tests/verification.rs
@@ -77,6 +77,56 @@ fn patch_test_every_element() {
     }
 }
 
+/// The patch test and uniaxial tension on a conforming hex8–pyramid5–tet4 transition mesh
+/// ([`transition_grid`]), distorted inside: the linear field is reproduced exactly, so the three
+/// kinds are compatible across their shared faces.
+#[test]
+fn hex_pyramid_tet_transition_passes_the_patch_test() {
+    let a = [[1e-3, 2e-4, -3e-4], [5e-4, -2e-3, 1e-4], [-1e-4, 3e-4, 7e-4]];
+    let field = |x: [f64; 3]| [0, 1, 2].map(|i| (0..3).map(|j| a[i][j] * x[j]).sum::<f64>());
+    let map = |p: [f64; 3]| {
+        let b = (PI * p[0]).sin() * (PI * p[1]).sin() * (PI * p[2]).sin();
+        [p[0] + 0.08 * b, p[1] - 0.06 * b, p[2] + 0.07 * b]
+    };
+    let mut g = transition_grid([4, 2, 2], map);
+    let kinds: std::collections::BTreeSet<&str> = g.model.elements().values().map(|e| e.kind.name()).collect();
+    assert_eq!(kinds.into_iter().collect::<Vec<_>>(), ["hex8", "pyramid5", "tet4"]);
+    g.define(iso(1000.0, 0.25), None, full(), Procedure::Static {});
+    let on_boundary = |p: [f64; 3]| p.iter().any(|c| near(*c, 0.0) || near(*c, 1.0));
+    for (i, label) in g.nodes_where(on_boundary).into_iter().enumerate() {
+        let u = field(g.model.node(label).unwrap().xyz);
+        let set = format!("b{i}");
+        g.model.add_set(&set, fieldcraft_model::Set { nodes: vec![label], ..Default::default() }).unwrap();
+        g.bc(
+            &set,
+            BoundaryCondition::Displacement { target: Target::Set(set.clone()), values: [Some(u[0]), Some(u[1]), Some(u[2]), None, None, None] },
+        );
+    }
+    let r = static_solve(&g);
+    let mut worst: f64 = 0.0;
+    for (label, u) in r.nodes.iter().zip(&r.displacement) {
+        let want = field(g.model.node(*label).unwrap().xyz);
+        worst = worst.max((0..3).map(|k| (u[k] - want[k]).abs()).fold(0.0, f64::max));
+    }
+    eprintln!("transition patch: displacement error {worst:.1e}");
+    assert!(worst < 1e-12, "{worst}");
+
+    // Uniaxial tension with symmetry planes: σxx = p and ux = p x / E everywhere.
+    let (young, nu, p) = (2000.0, 0.3, 5.0);
+    let mut g = transition_grid([4, 2, 2], map);
+    g.define(iso(young, nu), None, full(), Procedure::Static {});
+    g.symmetry("x0", Symmetry::X, |q| near(q[0], 0.0));
+    g.symmetry("y0", Symmetry::Y, |q| near(q[1], 0.0));
+    g.symmetry("z0", Symmetry::Z, |q| near(q[2], 0.0));
+    g.face_set("end", |q| near(q[0], 1.0));
+    g.load("pull", Load::Pressure { set: "end".into(), magnitude: -p });
+    let r = static_solve(&g);
+    let worst = r.nodes.iter().zip(&r.displacement).map(|(l, u)| (u[0] - p * g.model.node(*l).unwrap().xyz[0] / young).abs()).fold(0.0, f64::max);
+    let s_err = r.stress.iter().map(|s| (s[0] - p).abs()).fold(0.0, f64::max);
+    eprintln!("transition tension: displacement error {worst:.1e}, stress error {s_err:.1e}");
+    assert!(worst < 1e-11 && s_err < 1e-9, "{worst} {s_err}");
+}
+
 /// Uniaxial tension by a consistent surface load on a distorted block with symmetry planes:
 /// σxx = p exactly, ux = p L / E, the lateral contraction −ν p y / E, and the reactions balance
 /// the load (the force form of the patch test, and a check of pressure loads).
diff --git a/crates/ui-egui/src/theme.rs b/crates/ui-egui/src/theme.rs
index f30833f..bc3894f 100644
--- a/crates/ui-egui/src/theme.rs
+++ b/crates/ui-egui/src/theme.rs
@@ -57,6 +57,7 @@ pub struct Tokens {
     pub hex: Color32,
     pub tet: Color32,
     pub wedge: Color32,
+    pub pyramid: Color32,
     pub shell: Color32,
     pub line_element: Color32,
     /// Element edges (the mesh wireframe).
@@ -144,6 +145,7 @@ impl Tokens {
             hex: Color32::from_rgb(0x78, 0x9c, 0xc0),
             tet: Color32::from_rgb(0x80, 0xb4, 0x8a),
             wedge: Color32::from_rgb(0xc0, 0xaa, 0x78),
+            pyramid: Color32::from_rgb(0xa8, 0x90, 0xc8),
             shell: Color32::from_rgb(0xc0, 0x8c, 0x78),
             line_element: Color32::from_rgb(0xc8, 0x8c, 0x1e),
             mesh_edge: Color32::from_rgb(0x1e, 0x23, 0x2a),
@@ -210,6 +212,7 @@ impl Tokens {
         t.hex = Color32::from_rgb(0x82, 0xa6, 0xc8);
         t.tet = Color32::from_rgb(0x8c, 0xc0, 0x96);
         t.wedge = Color32::from_rgb(0xc8, 0xb4, 0x82);
+        t.pyramid = Color32::from_rgb(0xb4, 0x9c, 0xd2);
         t.shell = Color32::from_rgb(0xc8, 0x96, 0x82);
         t.line_element = Color32::from_rgb(0xf0, 0xc0, 0x50);
         t.mesh_edge = Color32::from_rgb(0x14, 0x18, 0x1e);
@@ -265,6 +268,7 @@ impl Tokens {
             hex: rgb(self.hex),
             tet: rgb(self.tet),
             wedge: rgb(self.wedge),
+            pyramid: rgb(self.pyramid),
             shell: rgb(self.shell),
             line_element: rgb(self.line_element),
             mesh_edge: rgb(self.mesh_edge),
@@ -287,6 +291,7 @@ impl Tokens {
             ElementKind::Hex8 | ElementKind::Hex20 => self.hex,
             ElementKind::Tet4 | ElementKind::Tet10 => self.tet,
             ElementKind::Wedge6 => self.wedge,
+            ElementKind::Pyramid5 => self.pyramid,
             ElementKind::Tri3 | ElementKind::Tri6 | ElementKind::Quad4 | ElementKind::Quad8 => self.shell,
             ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => self.line_element,
         }
diff --git a/docs/elements.md b/docs/elements.md
index 19ec97a..36d4315 100644
--- a/docs/elements.md
+++ b/docs/elements.md
@@ -72,7 +72,8 @@ hex8/20/27 with first and second derivatives).
 | hex8 (full, incompatible modes), hex20 (full, reduced) | Structured and swept meshes | done |
 | wedge6 | Swept meshes from triangles, transitions | done |
 | hex8 one-point with hourglass control, nodal-pressure tet4 (explicit) | Explicit dynamics | done |
-| wedge15, pyramid5/13 | Quadratic transitions, tet–hex transition | planned (M5); wedge15 and pyramid5 shape functions done |
+| pyramid5 | Conforming tet–hex transition (a hexahedron's face to tetrahedra) | done (Bedrosian's rational functions, a collapsed 2 × 2 × 2 Gauss rule; C3D5, CPYRAM, LS-DYNA degenerate solid) |
+| wedge15, pyramid13 | Quadratic transitions | planned (M5); wedge15 shape functions done |
 | Hybrid (mixed u–p) solids | Nearly incompressible materials (rubber) | done for hex8 (constant pressure), hex20 (linear pressure) and tet10 (constant pressure), discontinuous pressure condensed, isotropic materials; tet4/wedge6 refused (they lock); continuous-pressure (Taylor–Hood) elements planned with an indefinite factorization |
 | Modified quadratic tet (tet10 for contact/plasticity) | Robust quadratic tets | planned |
 
diff --git a/docs/fem.md b/docs/fem.md
index 60a644d..15c8c05 100644
--- a/docs/fem.md
+++ b/docs/fem.md
@@ -76,7 +76,8 @@ rational pyramid, 2 for the quadratic ones).
 Evaluation at a natural point `xi: [f64; 3]` (unused coordinates are zero):
 
 - `eval(xi) -> Shape`: values `n` and first derivatives `dn[i] = ∂Nᵢ/∂(ξ, η, ζ)`. tet4, tet10,
-  wedge6, hex8, hex20 use the hand-written functions of `shape.rs` (the solid kernels' path).
+  wedge6, hex8, hex20 use the hand-written functions of `shape.rs` (the solid kernels' path);
+  the solid kernel's pyramid5 uses this module's.
 - `eval2(xi) -> Shape2`: adds `d2n[i][a][b] = ∂²Nᵢ/∂ξₐ∂ξ_b`, from second-order forward-mode
   automatic differentiation of the same formulas (`jet::Jet`): exact to round-off.
 - `gradients(&shape, xyz) -> Gradients`: physical `dx[i] = ∇ₓNᵢ = J⁻¹ ∇_ξNᵢ` and `det J`, with
diff --git a/docs/formats/fieldcraft.md b/docs/formats/fieldcraft.md
index 6e4202b..2d94c84 100644
--- a/docs/formats/fieldcraft.md
+++ b/docs/formats/fieldcraft.md
@@ -61,7 +61,7 @@ tagged nodes × { u32 node, u32 tag }
 ```
 
 Element kind codes: 1 rod2, 2 beam2, 3 tri3, 4 tri6, 5 quad4, 6 quad8, 7 tet4, 8 tet10, 9 wedge6,
-10 hex8, 11 hex20. Node order is the model's (Abaqus/VTK). A tag indexes the entity list, which
+10 hex8, 11 hex20, 12 rod3, 13 pyramid5. Node order is the model's (Abaqus/VTK). A tag indexes the entity list, which
 lists each geometry entity once. Coordinates keep every bit (signed zero included).
 
 ## Reading a file safely