Implicit solver: reference points, couplings, moments, wedge6 (M4) (#26)

* fem: wedge6 element (triangle × interval, 3 × 2 points), patch and cantilever tests

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

* solver: reference points, kinematic and distributing couplings, moments by constraint elimination

Coupled reference points get six degrees of freedom; kinematic couplings are rigid
links, distributing couplings least-squares interpolation constraints. Every
constraint row and prescribed value is eliminated master-slave (Abel-Shephard
pivot choice), so the reduced matrix stays SPD; no coupling gives the old
elimination bit for bit. Verified: rigid motion of a coupled block, clamping
through a point, distributing tension vs pressure, bending and tube torsion.

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

* engine: reference point results in solve.static and solve.probe

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

* docs: couplings, moments and wedge6 in the solver docs, literature and roadmap

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

* solver: simplify constraint elimination and reference point plumbing

Rows combine by one sort; T_f transposed once; K and M reduced in one pass on a
shared pattern; prescribed motion computed once; reference point loads and BCs go
through Resolved like node ones (Mesh::label_dofs); one dof-direction rule;
StepResult::locate for probes; dead Csr::restrict removed.

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

* solver: refuse oversized reductions before allocating them

A held distributing coupling point makes a dense (3N)^2 block; estimate it from
the constraint rows and refuse it against the factorization's memory budget
with a clear message, before anything is built. Other reductions with
dependents count their entries first (stopping at the budget) and reserve
fallibly.

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

---------

Co-authored-by: jyviko <2501648+jyviko@users.noreply.github.com>
Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

crates/fem/src/element.rs

@@ -7,7 +7,8 @@
 //!
 //! Formulations:
 //! - **full integration**: the rule that integrates the stiffness of an undistorted element
-//!   exactly (1 point for tet4, 4 for tet10, 2 × 2 × 2 for hex8, 3 × 3 × 3 for hex20);
+//!   exactly (1 point for tet4, 4 for tet10, 3 triangle points × 2 for wedge6, 2 × 2 × 2 for
+//!   hex8, 3 × 3 × 3 for hex20);
 //! - **reduced integration** for hex20 (2 × 2 × 2), which avoids locking at the cost of one
 //!   spurious mode that a mesh of more than one element restrains. Reduced hex8 needs hourglass
 //!   control and is not offered here yet (it arrives with the explicit solver, Flanagan &
@@ -43,6 +44,7 @@ use crate::{Error, Result};
 #[derive(Debug, Clone, Copy, PartialEq, Eq)]
 pub enum Family {
     Tetrahedron,
+    Wedge,
     Hexahedron,
 }
 
@@ -67,7 +69,7 @@ pub struct ElementType {
 }
 
 /// The element table.
-pub static ELEMENT_TYPES: [ElementType; 4] = [
+pub static ELEMENT_TYPES: [ElementType; 5] = [
     ElementType {
         kind: ElementKind::Tet4,
         family: Family::Tetrahedron,
@@ -90,6 +92,17 @@ pub static ELEMENT_TYPES: [ElementType; 4] = [
         mass: Rule::TetConical(4),
         incompatible_modes: false,
     },
+    ElementType {
+        kind: ElementKind::Wedge6,
+        family: Family::Wedge,
+        nodes: 6,
+        corners: 6,
+        shape: shape::wedge6,
+        full: Rule::Wedge,
+        reduced: None,
+        mass: Rule::Wedge,
+        incompatible_modes: false,
+    },
     ElementType {
         kind: ElementKind::Hex8,
         family: Family::Hexahedron,
@@ -119,7 +132,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, hex8 and hex20 so far)", kind.name())))
+        .ok_or_else(|| Error::NotSupported(format!("{} elements (the solver has tet4, tet10, wedge6, hex8 and hex20 so far)", kind.name())))
 }
 
 /// How a solid element is formulated.
@@ -565,6 +578,7 @@ pub fn extrapolation(ty: &ElementType, rule: Rule) -> Result<Extrapolation> {
     let c = ty.corners;
     let corner_basis = |xi: [f64; 3]| match ty.family {
         Family::Tetrahedron => shape::tet4(xi),
+        Family::Wedge => shape::wedge6(xi),
         Family::Hexahedron => shape::hex8(xi),
     };
     let mut corner_rows = vec![0.0; c * g];

crates/fem/src/lib.rs

@@ -4,7 +4,7 @@
 //! global state, so the implicit solver (`modelcraft-solver`) and the explicit solver that follows
 //! it use the same element technology:
 //! - [`quadrature`]: Gauss–Legendre and tetrahedron rules;
-//! - [`shape`]: shape functions of tet4, tet10, hex8, hex20 and of their faces;
+//! - [`shape`]: shape functions of tet4, tet10, wedge6, hex8, hex20;
 //! - [`material`]: isotropic and orthotropic elasticity matrices, von Mises stress;
 //! - [`element`]: the element table, stiffness (full, reduced, incompatible modes), consistent
 //!   and lumped mass, body forces, integration-point stresses and their extrapolation to nodes;

crates/fem/src/quadrature.rs

@@ -12,6 +12,9 @@
 //! - On the reference triangle {ξ, η ≥ 0, ξ + η ≤ 1}: the same collapsed product in two
 //!   dimensions (Jacobian 1 − u), exact for degree 2n − 2.
 //! - On the reference square [−1, 1]²: tensor Gauss–Legendre.
+//! - On the reference wedge (triangle × [−1, 1]): the symmetric three-point triangle rule
+//!   (points at area coordinates (⅔, ⅙, ⅙) and permutations, degree 2; Zienkiewicz, Taylor &
+//!   Zhu 2013, table 6.3) times two-point Gauss–Legendre in ζ (degree 3).
 
 /// A quadrature point: natural coordinates and weight.
 #[derive(Debug, Clone, Copy, PartialEq)]
@@ -54,6 +57,9 @@ pub enum Rule {
     Quad(u8),
     /// The collapsed product on the reference triangle, n points per direction (z = 0).
     TriConical(u8),
+    /// 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,
 }
 
 impl Rule {
@@ -103,6 +109,11 @@ impl Rule {
                     .flat_map(|&(u, wu)| g.iter().map(move |&(v, wv)| Point { xi: [u, v * (1.0 - u), 0.0], weight: wu * wv * (1.0 - u) }))
                     .collect()
             }
+            Rule::Wedge => {
+                let (a, b) = (2.0 / 3.0, 1.0 / 6.0);
+                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()
+            }
         }
     }
 }
@@ -149,8 +160,18 @@ mod tests {
                 assert!((got - exact_1d * 2.0).abs() < 1e-14, "quad {n} y^{a}");
             }
         }
-        // Triangle: ∫ ξ^a η^b = a! b! / (a + b + 2)!.
+        // Wedge: ∫ ξ^a η^b ζ^c = a! b! / (a + b + 2)! · ∫ ζ^c, for a + b ≤ 2 and c ≤ 3.
         let f = |n: i32| (1..=n).map(f64::from).product::<f64>();
+        for a in 0..=2 {
+            for b in 0..=2 - a {
+                for c in 0..=3 {
+                    let line = if c % 2 == 1 { 0.0 } else { 2.0 / f64::from(c + 1) };
+                    let got = integrate(Rule::Wedge, a, b, c);
+                    assert!((got - f(a) * f(b) / f(a + b + 2) * line).abs() < 1e-15, "wedge ξ^{a} η^{b} ζ^{c}");
+                }
+            }
+        }
+        // Triangle: ∫ ξ^a η^b = a! b! / (a + b + 2)!.
         for a in 0..=4 {
             for b in 0..=4 - a {
                 let got = integrate(Rule::TriConical(3), a, b, 0);

crates/fem/src/shape.rs

@@ -9,7 +9,11 @@
 //!   N = ⅛(1 + ξξᵢ)(1 + ηηᵢ)(1 + ζζᵢ)(ξξᵢ + ηηᵢ + ζζᵢ − 2), a mid-side node with ξᵢ = 0
 //!   N = ¼(1 − ξ²)(1 + ηηᵢ)(1 + ζζᵢ), and likewise), mid-side nodes 9–12 on the bottom edges,
 //!   13–16 on the top edges, 17–20 on the vertical edges (mid-side positions come from
-//!   [`modelcraft_model::ElementKind::edges`], the one place the node order is written down).
+//!   [`modelcraft_model::ElementKind::edges`], the one place the node order is written down);
+//! - the six-node wedge (pentahedron) as the product of the linear triangle {ξ, η ≥ 0,
+//!   ξ + η ≤ 1} (area coordinates L₁ = 1 − ξ − η, L₂ = ξ, L₃ = η) and the linear interval
+//!   ζ ∈ [−1, 1]: nodes 1–3 at ζ = −1 counter-clockwise seen from above, nodes 4–6 above them,
+//!   Nᵢ = Lᵢ (1 − ζ)/2 and Nᵢ₊₃ = Lᵢ (1 + ζ)/2.
 //!
 //! The functions are the textbook ones (O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, *The Finite
 //! Element Method: Its Basis and Fundamentals*, 7th ed., Butterworth-Heinemann 2013, ch. 6;
@@ -47,6 +51,9 @@ impl Shape {
 pub const HEX_CORNERS: [[f64; 3]; 8] =
     [[-1., -1., -1.], [1., -1., -1.], [1., 1., -1.], [-1., 1., -1.], [-1., -1., 1.], [1., -1., 1.], [1., 1., 1.], [-1., 1., 1.]];
 
+/// Natural coordinates of the wedge's corners.
+pub const WEDGE_CORNERS: [[f64; 3]; 6] = [[0., 0., -1.], [1., 0., -1.], [0., 1., -1.], [0., 0., 1.], [1., 0., 1.], [0., 1., 1.]];
+
 /// Natural coordinates of the tetrahedron's corners.
 pub const TET_CORNERS: [[f64; 3]; 4] = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]];
 
@@ -95,6 +102,20 @@ pub fn hex8(xi: [f64; 3]) -> Shape {
     s
 }
 
+/// Six-node wedge: the linear triangle times the linear interval.
+pub fn wedge6(xi: [f64; 3]) -> Shape {
+    let mut s = Shape::new(6);
+    let l = [1.0 - xi[0] - xi[1], xi[0], xi[1]];
+    let dl = [[-1.0, -1.0], [1.0, 0.0], [0.0, 1.0]];
+    for (h, f, df) in [(0, (1.0 - xi[2]) / 2.0, -0.5), (3, (1.0 + xi[2]) / 2.0, 0.5)] {
+        for i in 0..3 {
+            s.n[h + i] = l[i] * f;
+            s.dn[h + i] = [dl[i][0] * f, dl[i][1] * f, l[i] * df];
+        }
+    }
+    s
+}
+
 /// Twenty-node serendipity hexahedron.
 pub fn hex20(xi: [f64; 3]) -> Shape {
     let mut s = Shape::new(20);
@@ -147,6 +168,7 @@ mod tests {
                 v
             }
             "hex8" => HEX_CORNERS.to_vec(),
+            "wedge6" => WEDGE_CORNERS.to_vec(),
             _ => {
                 let mut v = HEX_CORNERS.to_vec();
                 for e in ElementKind::Hex20.edges() {
@@ -159,7 +181,7 @@ mod tests {
 
     #[test]
     fn kronecker_delta_partition_of_unity_and_derivatives() {
-        let kinds: [(&str, fn([f64; 3]) -> Shape); 4] = [("tet4", tet4), ("tet10", tet10), ("hex8", hex8), ("hex20", hex20)];
+        let kinds: [(&str, fn([f64; 3]) -> Shape); 5] = [("tet4", tet4), ("tet10", tet10), ("wedge6", wedge6), ("hex8", hex8), ("hex20", hex20)];
         for (name, f) in kinds {
             let nodes = nodes_of(name);
             for (i, p) in nodes.iter().enumerate() {

crates/fem/src/tests.rs

@@ -12,6 +12,7 @@ use crate::material::{isotropic, orthotropic, stress};
 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]],
         _ => vec![
             [0.0, 0.0, 0.0],
             [2.1, 0.1, -0.1],
@@ -33,9 +34,10 @@ fn distorted(kind: ElementKind, curved: bool) -> Vec<[f64; 3]> {
     xyz
 }
 
-const CASES: [(ElementKind, Formulation); 6] = [
+const CASES: [(ElementKind, Formulation); 7] = [
     (ElementKind::Tet4, Formulation::Full),
     (ElementKind::Tet10, Formulation::Full),
+    (ElementKind::Wedge6, Formulation::Full),
     (ElementKind::Hex8, Formulation::Full),
     (ElementKind::Hex8, Formulation::IncompatibleModes),
     (ElementKind::Hex20, Formulation::Full),
@@ -145,6 +147,7 @@ fn extrapolation_reproduces_linear_fields_on_affine_elements() {
     let affine = |kind: ElementKind| -> Vec<[f64; 3]> {
         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.]],
             _ => 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();
@@ -172,7 +175,8 @@ fn bad_input_is_an_error() {
     assert!(Solid::new(ElementKind::Hex8, &xyz[..7], Formulation::Full).is_err());
     assert!(Solid::new(ElementKind::Hex8, &xyz, Formulation::Reduced).is_err());
     assert!(Solid::new(ElementKind::Tet4, &xyz[..4], Formulation::IncompatibleModes).is_err());
-    assert!(Solid::new(ElementKind::Wedge6, &xyz[..6], Formulation::Full).is_err());
+    assert!(Solid::new(ElementKind::Quad4, &xyz[..4], Formulation::Full).is_err());
+    assert!(Solid::new(ElementKind::Wedge6, &xyz[..6], Formulation::IncompatibleModes).is_err());
     let mut nan = xyz.clone();
     nan[3][1] = f64::NAN;
     assert!(Solid::new(ElementKind::Hex8, &nan, Formulation::Full).is_err());