diff --git a/crates/fem/src/dual.rs b/crates/fem/src/dual.rs
new file mode 100644
index 0000000..6090037
--- /dev/null
+++ b/crates/fem/src/dual.rs
@@ -0,0 +1,341 @@
+//! Forward-mode automatic differentiation in [`VARS`] variables, for kernels whose exact first
+//! and second derivatives are too long to write by hand (the mortar contact integrals of
+//! [`crate::mortar`], whose segmentation moves with the nodes).
+//!
+//! A kernel is written once over a [`Real`] scalar and run on `f64` for its value, on [`Dual`]
+//! for its gradient, or on [`Hyper`] for its gradient and Hessian. Arithmetic propagates the
+//! derivatives by the chain rule (truncated Taylor arithmetic: A. Griewank, A. Walther,
+//! *Evaluating Derivatives*, 2nd ed., SIAM 2008, ch. 3 and 13), so the derivatives are those of
+//! the code as it runs, to round-off: the branches a kernel takes (which polygon edges clip,
+//! which root of a quadratic) are decided on the values and differentiated where they are
+//! taken, which is the derivative of the piecewise smooth function wherever it is smooth. The
+//! automated linearisation of contact kernels this way is that of J. Lengiewicz, J. Korelc, S.
+//! Stupkiewicz, *Automation of finite element formulations for large deformation contact
+//! problems*, IJNME 85:1252–1279, 2011.
+//!
+//! [`jet::Jet`](crate::jet::Jet) is the three-variable sibling the shape-function library uses.
+
+use std::ops::{Add, Div, Mul, Neg, Sub};
+
+/// Independent variables of [`Dual`] and [`Hyper`]: the corners of two faces of up to four
+/// nodes, three coordinates each.
+pub const VARS: usize = 24;
+
+/// Entries of a packed symmetric [`VARS`] × [`VARS`] Hessian (its upper triangle by rows).
+pub const PACKED: usize = VARS * (VARS + 1) / 2;
+
+/// The position of (i, j) in a packed upper triangle (either order).
+pub fn packed(i: usize, j: usize) -> usize {
+    let (a, b) = if i <= j { (i, j) } else { (j, i) };
+    // Row a starts after rows 0..a, of VARS, VARS − 1, … entries.
+    a * VARS - a * a.saturating_sub(1) / 2 + (b - a)
+}
+
+/// A scalar a kernel can be written over: `f64`, [`Dual`] or [`Hyper`].
+pub trait Real:
+    Copy
+    + Add<Output = Self>
+    + Sub<Output = Self>
+    + Mul<Output = Self>
+    + Div<Output = Self>
+    + Neg<Output = Self>
+    + Add<f64, Output = Self>
+    + Sub<f64, Output = Self>
+    + Mul<f64, Output = Self>
+{
+    /// A constant.
+    fn cst(v: f64) -> Self;
+    /// The value.
+    fn val(&self) -> f64;
+    fn sqrt(self) -> Self;
+}
+
+impl Real for f64 {
+    fn cst(v: f64) -> f64 {
+        v
+    }
+    fn val(&self) -> f64 {
+        *self
+    }
+    fn sqrt(self) -> f64 {
+        f64::sqrt(self)
+    }
+}
+
+/// A value and its gradient.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Dual {
+    pub v: f64,
+    pub g: [f64; VARS],
+}
+
+impl Dual {
+    /// Variable `k` at value `v` (a constant when `k` is out of range).
+    pub fn var(v: f64, k: usize) -> Dual {
+        let mut d = Dual::cst(v);
+        if let Some(g) = d.g.get_mut(k) {
+            *g = 1.0;
+        }
+        d
+    }
+
+    fn chain(self, v: f64, d: f64) -> Dual {
+        Dual { v, g: self.g.map(|x| d * x) }
+    }
+}
+
+impl Real for Dual {
+    fn cst(v: f64) -> Dual {
+        Dual { v, g: [0.0; VARS] }
+    }
+    fn val(&self) -> f64 {
+        self.v
+    }
+    fn sqrt(self) -> Dual {
+        let s = self.v.sqrt();
+        self.chain(s, 0.5 / s)
+    }
+}
+
+impl Add for Dual {
+    type Output = Dual;
+    fn add(mut self, o: Dual) -> Dual {
+        self.v += o.v;
+        self.g.iter_mut().zip(o.g).for_each(|(a, b)| *a += b);
+        self
+    }
+}
+
+impl Sub for Dual {
+    type Output = Dual;
+    fn sub(mut self, o: Dual) -> Dual {
+        self.v -= o.v;
+        self.g.iter_mut().zip(o.g).for_each(|(a, b)| *a -= b);
+        self
+    }
+}
+
+impl Mul for Dual {
+    type Output = Dual;
+    fn mul(self, o: Dual) -> Dual {
+        let mut g = [0.0; VARS];
+        for k in 0..VARS {
+            g[k] = self.v * o.g[k] + o.v * self.g[k];
+        }
+        Dual { v: self.v * o.v, g }
+    }
+}
+
+impl Div for Dual {
+    type Output = Dual;
+    fn div(self, o: Dual) -> Dual {
+        let r = 1.0 / o.v;
+        let v = self.v * r;
+        let mut g = [0.0; VARS];
+        for k in 0..VARS {
+            g[k] = (self.g[k] - v * o.g[k]) * r;
+        }
+        Dual { v, g }
+    }
+}
+
+impl Neg for Dual {
+    type Output = Dual;
+    fn neg(self) -> Dual {
+        Dual { v: -self.v, g: self.g.map(|x| -x) }
+    }
+}
+
+impl Add<f64> for Dual {
+    type Output = Dual;
+    fn add(mut self, o: f64) -> Dual {
+        self.v += o;
+        self
+    }
+}
+
+impl Sub<f64> for Dual {
+    type Output = Dual;
+    fn sub(mut self, o: f64) -> Dual {
+        self.v -= o;
+        self
+    }
+}
+
+impl Mul<f64> for Dual {
+    type Output = Dual;
+    fn mul(self, o: f64) -> Dual {
+        Dual { v: self.v * o, g: self.g.map(|x| x * o) }
+    }
+}
+
+/// A value, its gradient and its Hessian (packed upper triangle, [`packed`]).
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Hyper {
+    pub v: f64,
+    pub g: [f64; VARS],
+    pub h: [f64; PACKED],
+}
+
+impl Hyper {
+    /// Variable `k` at value `v` (a constant when `k` is out of range).
+    pub fn var(v: f64, k: usize) -> Hyper {
+        let mut d = Hyper::cst(v);
+        if let Some(g) = d.g.get_mut(k) {
+            *g = 1.0;
+        }
+        d
+    }
+
+    /// ∂²/∂x_i∂x_j.
+    pub fn hessian(&self, i: usize, j: usize) -> f64 {
+        self.h.get(packed(i, j)).copied().unwrap_or(0.0)
+    }
+
+    /// f(self) from f's value `f0`, first derivative `f1` and second derivative `f2` at the
+    /// value: g = f1 g, H = f1 H + f2 g gᵀ.
+    fn chain(&self, f0: f64, f1: f64, f2: f64) -> Hyper {
+        let mut out = Hyper { v: f0, g: self.g.map(|x| f1 * x), h: [0.0; PACKED] };
+        let mut at = 0;
+        for i in 0..VARS {
+            let gi = f2 * self.g[i];
+            for j in i..VARS {
+                out.h[at] = f1 * self.h[at] + gi * self.g[j];
+                at += 1;
+            }
+        }
+        out
+    }
+}
+
+impl Real for Hyper {
+    fn cst(v: f64) -> Hyper {
+        Hyper { v, g: [0.0; VARS], h: [0.0; PACKED] }
+    }
+    fn val(&self) -> f64 {
+        self.v
+    }
+    fn sqrt(self) -> Hyper {
+        let s = self.v.sqrt();
+        self.chain(s, 0.5 / s, -0.25 / (s * s * s))
+    }
+}
+
+impl Add for Hyper {
+    type Output = Hyper;
+    fn add(mut self, o: Hyper) -> Hyper {
+        self.v += o.v;
+        self.g.iter_mut().zip(&o.g).for_each(|(a, b)| *a += b);
+        self.h.iter_mut().zip(&o.h).for_each(|(a, b)| *a += b);
+        self
+    }
+}
+
+impl Sub for Hyper {
+    type Output = Hyper;
+    fn sub(mut self, o: Hyper) -> Hyper {
+        self.v -= o.v;
+        self.g.iter_mut().zip(&o.g).for_each(|(a, b)| *a -= b);
+        self.h.iter_mut().zip(&o.h).for_each(|(a, b)| *a -= b);
+        self
+    }
+}
+
+impl Mul for Hyper {
+    type Output = Hyper;
+    fn mul(self, o: Hyper) -> Hyper {
+        let mut out = Hyper { v: self.v * o.v, g: [0.0; VARS], h: [0.0; PACKED] };
+        for k in 0..VARS {
+            out.g[k] = self.v * o.g[k] + o.v * self.g[k];
+        }
+        let mut at = 0;
+        for i in 0..VARS {
+            let (a, b) = (self.g[i], o.g[i]);
+            for j in i..VARS {
+                out.h[at] = self.v * o.h[at] + o.v * self.h[at] + a * o.g[j] + b * self.g[j];
+                at += 1;
+            }
+        }
+        out
+    }
+}
+
+impl Div for Hyper {
+    type Output = Hyper;
+    fn div(self, o: Hyper) -> Hyper {
+        let r = 1.0 / o.v;
+        self * o.chain(r, -r * r, 2.0 * r * r * r)
+    }
+}
+
+impl Neg for Hyper {
+    type Output = Hyper;
+    fn neg(self) -> Hyper {
+        Hyper { v: -self.v, g: self.g.map(|x| -x), h: self.h.map(|x| -x) }
+    }
+}
+
+impl Add<f64> for Hyper {
+    type Output = Hyper;
+    fn add(mut self, o: f64) -> Hyper {
+        self.v += o;
+        self
+    }
+}
+
+impl Sub<f64> for Hyper {
+    type Output = Hyper;
+    fn sub(mut self, o: f64) -> Hyper {
+        self.v -= o;
+        self
+    }
+}
+
+impl Mul<f64> for Hyper {
+    type Output = Hyper;
+    fn mul(self, o: f64) -> Hyper {
+        Hyper { v: self.v * o, g: self.g.map(|x| x * o), h: self.h.map(|x| x * o) }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn f<T: Real>(x: T, y: T) -> T {
+        (x * x * y + T::cst(1.0)).sqrt() / (y - x * 0.5) + x * 3.0
+    }
+
+    #[test]
+    fn packed_indices_cover_the_upper_triangle_once() {
+        let mut seen = vec![false; PACKED];
+        let mut at = 0;
+        for i in 0..VARS {
+            for j in i..VARS {
+                assert_eq!(packed(i, j), at);
+                assert_eq!(packed(j, i), at);
+                seen[at] = true;
+                at += 1;
+            }
+        }
+        assert!(seen.iter().all(|s| *s));
+    }
+
+    #[test]
+    fn derivatives_match_finite_differences() {
+        let (x, y) = (0.7, 2.3);
+        let h = 1e-5;
+        let fx = |x: f64, y: f64| f(x, y);
+        let gx = (fx(x + h, y) - fx(x - h, y)) / (2.0 * h);
+        let gy = (fx(x, y + h) - fx(x, y - h)) / (2.0 * h);
+        let hxy = (fx(x + h, y + h) - fx(x + h, y - h) - fx(x - h, y + h) + fx(x - h, y - h)) / (4.0 * h * h);
+        let hxx = (fx(x + h, y) - 2.0 * fx(x, y) + fx(x - h, y)) / (h * h);
+        let d = f(Dual::var(x, 3), Dual::var(y, 17));
+        assert!((d.v - fx(x, y)).abs() < 1e-14 && (d.g[3] - gx).abs() < 1e-8 && (d.g[17] - gy).abs() < 1e-8);
+        let e = f(Hyper::var(x, 3), Hyper::var(y, 17));
+        assert!((e.g[3] - gx).abs() < 1e-8 && (e.g[17] - gy).abs() < 1e-8);
+        assert!((e.hessian(3, 17) - hxy).abs() < 1e-4 && (e.hessian(17, 3) - hxy).abs() < 1e-4, "{} {hxy}", e.hessian(3, 17));
+        assert!((e.hessian(3, 3) - hxx).abs() < 1e-3);
+        assert_eq!(e.hessian(0, 5), 0.0);
+    }
+}
diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs
index 43b4e4d..c87ac37 100644
--- a/crates/fem/src/lib.rs
+++ b/crates/fem/src/lib.rs
@@ -30,6 +30,8 @@
 //! - [`face`]: consistent nodal forces of pressure on element faces;
 //! - [`contact`]: contact geometry: the closest point of an element face to a point and the
 //!   bucket grid of the global search;
+//! - [`mortar`]: mortar segmentation of a secondary face against a main face and the integration
+//!   points on the overlap, over [`dual`]'s automatic differentiation for the exact tangent;
 //! - [`heat`]: heat conduction (conductivity, capacity, sources, surface flux and convection,
 //!   heat flux); thermal strain loads and stresses are in [`element`];
 //! - [`uniform`]: uniform-strain (one-point) hexahedron and tetrahedron kinematics, hourglass
@@ -59,6 +61,7 @@ pub mod constitutive;
 pub mod contact;
 pub mod dense;
 pub mod dof;
+pub mod dual;
 pub mod elastoplastic;
 pub mod element;
 pub mod face;
@@ -68,6 +71,7 @@ pub mod heat;
 pub mod hybrid;
 pub mod jet;
 pub mod material;
+pub mod mortar;
 pub mod plane;
 pub mod plasticity;
 pub mod quadrature;
diff --git a/crates/fem/src/mortar.rs b/crates/fem/src/mortar.rs
new file mode 100644
index 0000000..34702a1
--- /dev/null
+++ b/crates/fem/src/mortar.rs
@@ -0,0 +1,483 @@
+//! Mortar segmentation and integration: the overlap of a secondary (non-mortar) face with a main
+//! (mortar) face, and the integration points on it where both faces' shape functions and
+//! positions are known.
+//!
+//! **Segmentation** (M. A. Puso, T. A. Laursen, *A mortar segment-to-segment contact method for
+//! large deformation solid mechanics*, CMAME 193:601–629, 2004, §3; A. Popp, M. W. Gee, W. A.
+//! Wall, *A finite deformation mortar contact formulation using a primal–dual active set
+//! strategy*, IJNME 79:1354–1391, 2009, §3.3). Both faces are projected along the secondary
+//! face's unit normal n onto its plane (through its centroid), and the main polygon is clipped by
+//! the secondary one (Sutherland–Hodgman; a quadrilateral whose projection is not convex is
+//! split into its two triangles first). Each overlap polygon is split into a fan of triangles
+//! integrated with the seven-point degree-5 rule (J. Radon 1948; A. H. Stroud, *Approximate
+//! Calculation of Multiple Integrals*, Prentice-Hall 1971, T2:5-1), which integrates the products
+//! of two shape functions that the tie and the patch test need exactly on triangles and
+//! parallelograms (on a general quadrilateral the inverse of its map is not a polynomial, and the
+//! rule is accurate to its degree).
+//!
+//! **Points.** At an integration point y of the plane, each face's coordinates are found by
+//! inverting its projected map: the linear one of a triangle, or the bilinear one of a
+//! quadrilateral, x(u, v) = a + b u + c v + d u v, whose inverse is the root in [0, 1] of the
+//! quadratic (b × d) u² + (b × c − h × d) u − h × c = 0 (h = y − a, × the scalar cross product of
+//! the plane). The projection of a bilinear patch onto a plane is the bilinear patch of its
+//! projected corners, so the main face's point found this way is the point of the main surface
+//! along n through y, warped or not. The [`Sample`] then holds both faces' shape functions, both
+//! points in space, the normal and the weight (the plane's area element times the rule's weight).
+//!
+//! The kernel is written over [`Real`], so the same code gives the values (`f64`), their
+//! gradient ([`crate::dual::Dual`]) or their gradient and Hessian ([`crate::dual::Hyper`]) with
+//! respect to the corner coordinates, variables `3 i + k` for secondary corner i and
+//! `12 + 3 j + k` for main corner j: the exact linearisation of the segmentation, the normal and
+//! the projection, which the consistent tangent of mortar contact needs (Popp et al. 2009, §4).
+
+use crate::dual::Real;
+
+/// A point in space over a scalar.
+pub type Point<T> = [T; 3];
+
+/// A point of the secondary face's plane.
+type Flat<T> = [T; 2];
+
+/// Radon's seven-point rule on a triangle, exact for degree 5: barycentric coordinates of the
+/// second and third corners and weights summing to one.
+fn rule() -> [([f64; 2], f64); 7] {
+    let s = 15f64.sqrt();
+    let (a, b) = ((6.0 - s) / 21.0, (9.0 + 2.0 * s) / 21.0);
+    let (c, d) = ((6.0 + s) / 21.0, (9.0 - 2.0 * s) / 21.0);
+    let (wa, wc) = ((155.0 - s) / 1200.0, (155.0 + s) / 1200.0);
+    [([1.0 / 3.0, 1.0 / 3.0], 0.225), ([a, a], wa), ([b, a], wa), ([a, b], wa), ([c, c], wc), ([d, c], wc), ([c, d], wc)]
+}
+
+/// One integration point of a segment.
+#[derive(Debug, Clone, Copy)]
+pub struct Sample<T> {
+    /// The rule's weight times the plane's area element (the secondary face's projected area).
+    pub weight: T,
+    /// The secondary face's shape functions at the point (0 past its corners).
+    pub secondary: [T; 4],
+    /// The main face's shape functions at the point.
+    pub main: [T; 4],
+    /// The point on the secondary face, and on the main face along the normal.
+    pub xs: Point<T>,
+    pub xm: Point<T>,
+    /// The secondary face's unit normal (outward from its element).
+    pub normal: Point<T>,
+}
+
+impl<T: Real> Sample<T> {
+    /// The gap along the secondary normal, (x_m − x_s) · n: positive where the faces are apart.
+    pub fn gap(&self) -> T {
+        dot(sub(self.xm, self.xs), self.normal)
+    }
+}
+
+fn sub<T: Real>(a: Point<T>, b: Point<T>) -> Point<T> {
+    [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
+}
+
+fn dot<T: Real>(a: Point<T>, b: Point<T>) -> T {
+    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
+}
+
+fn cross<T: Real>(a: Point<T>, b: Point<T>) -> Point<T> {
+    [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]]
+}
+
+fn scale<T: Real>(a: Point<T>, s: T) -> Point<T> {
+    [a[0] * s, a[1] * s, a[2] * s]
+}
+
+fn cross2<T: Real>(a: Flat<T>, b: Flat<T>) -> T {
+    a[0] * b[1] - a[1] * b[0]
+}
+
+fn sub2<T: Real>(a: Flat<T>, b: Flat<T>) -> Flat<T> {
+    [a[0] - b[0], a[1] - b[1]]
+}
+
+/// A face's area vector: its outward normal times its area (flat triangle; for a bilinear
+/// quadrilateral half the cross product of its diagonals, the area of its projection on the
+/// plane normal to it).
+pub fn area_vector<T: Real>(x: &[Point<T>]) -> Option<Point<T>> {
+    let half = |v: Point<T>| scale(v, T::cst(0.5));
+    match x {
+        [a, b, c] => Some(half(cross(sub(*b, *a), sub(*c, *a)))),
+        [a, b, c, d] => Some(half(cross(sub(*c, *a), sub(*d, *b)))),
+        _ => None,
+    }
+}
+
+/// A projected face as convex counter-clockwise pieces for the clipping: itself when it is
+/// convex (no artificial diagonal, whose crossings with the other face's edges would change the
+/// overlap's polygons, and the integration's round-off, as the faces move), else its two
+/// triangles.
+fn pieces<T: Real>(points: &[Flat<T>]) -> Vec<Vec<Flat<T>>> {
+    let mut whole = points.to_vec();
+    if signed_area(&whole).val() < 0.0 {
+        whole.reverse();
+    }
+    let n = whole.len();
+    let convex = (0..n).all(|i| {
+        let (a, b, c) = (whole[i], whole[(i + 1) % n], whole[(i + 2) % n]);
+        cross2(sub2(b, a), sub2(c, b)).val() > 0.0
+    });
+    if convex || n == 3 {
+        return vec![whole];
+    }
+    [[0, 1, 2], [0, 2, 3]].iter().map(|t| counter_clockwise(t, points)).collect()
+}
+
+fn signed_area<T: Real>(p: &[Flat<T>]) -> T {
+    let mut s = T::cst(0.0);
+    for (i, a) in p.iter().enumerate() {
+        let b = p.get((i + 1) % p.len()).copied().unwrap_or(*a);
+        s = s + cross2(*a, b);
+    }
+    s * 0.5
+}
+
+/// The corners `t` of `points`, counter-clockwise.
+fn counter_clockwise<T: Real>(t: &[usize; 3], points: &[Flat<T>]) -> Vec<Flat<T>> {
+    let mut out: Vec<Flat<T>> = t.iter().filter_map(|i| points.get(*i).copied()).collect();
+    if signed_area(&out).val() < 0.0 {
+        out.reverse();
+    }
+    out
+}
+
+/// A vertex this close to a clip edge, relative to the edge's length, is on it.
+const ON_EDGE: f64 = 1e-8;
+
+/// `subject` clipped by the convex counter-clockwise polygon `clip` (Sutherland–Hodgman), the
+/// sides decided on the values. A vertex within [`ON_EDGE`] of a clip edge's line counts as on
+/// it (kept, and never cut): the intersection of two nearly coincident edges is ill-conditioned,
+/// its position and above all its derivatives swamped by round-off, while the vertex itself is
+/// as good a corner of the overlap to that precision. Edges are cut only between a vertex
+/// clearly inside and one clearly outside, where the intersection is well conditioned.
+fn clip<T: Real>(subject: &[Flat<T>], clip: &[Flat<T>]) -> Vec<Flat<T>> {
+    let mut out: Vec<Flat<T>> = subject.to_vec();
+    for (i, a) in clip.iter().enumerate() {
+        let Some(b) = clip.get((i + 1) % clip.len()).copied() else { break };
+        let input = std::mem::take(&mut out);
+        let edge = sub2(b, *a);
+        let eps = ON_EDGE * (edge[0].val() * edge[0].val() + edge[1].val() * edge[1].val());
+        let side = |p: Flat<T>| cross2(edge, sub2(p, *a));
+        for (j, p) in input.iter().enumerate() {
+            let Some(q) = input.get((j + 1) % input.len()).copied() else { break };
+            let (sp, sq) = (side(*p), side(q));
+            if sp.val() >= -eps {
+                out.push(*p);
+            }
+            if (sp.val() > eps && sq.val() < -eps) || (sp.val() < -eps && sq.val() > eps) {
+                let t = sp / (sp - sq);
+                out.push([p[0] + t * (q[0] - p[0]), p[1] + t * (q[1] - p[1])]);
+            }
+        }
+        if out.is_empty() {
+            break;
+        }
+    }
+    out
+}
+
+/// The shape functions of a face (of 3 or 4 projected corners `p`) at plane point `y`: the
+/// inverse of its linear or bilinear map. `None` when it has no inverse there.
+fn shape_at<T: Real>(p: &[Flat<T>], y: Flat<T>) -> Option<[T; 4]> {
+    let zero = T::cst(0.0);
+    match p {
+        [p0, p1, p2] => {
+            let (e1, e2, h) = (sub2(*p1, *p0), sub2(*p2, *p0), sub2(y, *p0));
+            let det = cross2(e1, e2);
+            if !(det.val().abs() > 0.0) {
+                return None;
+            }
+            let v = cross2(h, e2) / det;
+            let w = cross2(e1, h) / det;
+            Some([T::cst(1.0) - v - w, v, w, zero])
+        }
+        [p0, p1, p2, p3] => {
+            let b = sub2(*p1, *p0);
+            let c = sub2(*p3, *p0);
+            let d = [p0[0] - p1[0] + p2[0] - p3[0], p0[1] - p1[1] + p2[1] - p3[1]];
+            let h = sub2(y, *p0);
+            let k2 = cross2(b, d);
+            let k1 = cross2(b, c) - cross2(h, d);
+            let k0 = -cross2(h, c);
+            let u = if k2.val().abs() <= 1e-12 * k1.val().abs() {
+                if !(k1.val().abs() > 0.0) {
+                    return None;
+                }
+                -k0 / k1
+            } else {
+                // The roots of k2 u² + k1 u + k0, computed without cancellation; the one in [0, 1]
+                // (the nearer when round-off puts both outside).
+                let disc = k1 * k1 - k2 * k0 * 4.0;
+                let root = if disc.val() > 0.0 { disc.sqrt() } else { T::cst(0.0) };
+                let q = if k1.val() >= 0.0 { -(k1 + root) * 0.5 } else { (root - k1) * 0.5 };
+                let mut candidates = [None, None];
+                if q.val().abs() > 0.0 {
+                    candidates[0] = Some(q / k2);
+                    candidates[1] = Some(k0 / q);
+                } else {
+                    candidates[0] = Some(-k1 / (k2 * 2.0));
+                }
+                let distance = |u: f64| {
+                    if u < 0.0 {
+                        -u
+                    } else if u > 1.0 {
+                        u - 1.0
+                    } else {
+                        0.0
+                    }
+                };
+                candidates.into_iter().flatten().filter(|u| u.val().is_finite()).min_by(|a, b| distance(a.val()).total_cmp(&distance(b.val())))?
+            };
+            // v from the better conditioned component of h − u b = v (c + d u).
+            let m = [c[0] + d[0] * u, c[1] + d[1] * u];
+            let k = if m[0].val().abs() >= m[1].val().abs() { 0 } else { 1 };
+            if !(m[k].val().abs() > 0.0) {
+                return None;
+            }
+            let v = (h[k] - b[k] * u) / m[k];
+            let (s, t) = (u * 2.0 - 1.0, v * 2.0 - 1.0);
+            let one = T::cst(1.0);
+            Some([(one - s) * (one - t) * 0.25, (one + s) * (one - t) * 0.25, (one + s) * (one + t) * 0.25, (one - s) * (one + t) * 0.25])
+        }
+        _ => None,
+    }
+}
+
+fn interpolate<T: Real>(x: &[Point<T>], n: &[T; 4]) -> Point<T> {
+    let mut out = [T::cst(0.0); 3];
+    for (xi, ni) in x.iter().zip(n) {
+        for k in 0..3 {
+            out[k] = out[k] + xi[k] * *ni;
+        }
+    }
+    out
+}
+
+/// Integrate over the overlap of secondary face `secondary` with main face `main` (3 or 4
+/// corners each, in their faces' node order): `f` is called at every integration point. `None`
+/// when a face has another number of corners or the secondary face has no normal (collapsed);
+/// faces that do not overlap give no point.
+pub fn integrate<T: Real>(secondary: &[Point<T>], main: &[Point<T>], mut f: impl FnMut(&Sample<T>)) -> Option<()> {
+    if !matches!(main.len(), 3 | 4) {
+        return None;
+    }
+    let area = area_vector(secondary)?;
+    let length = dot(area, area).sqrt();
+    if !(length.val() > 0.0 && length.val().is_finite()) {
+        return None;
+    }
+    let normal = scale(area, T::cst(1.0) / length);
+    let count = T::cst(secondary.len() as f64);
+    let mut center = [T::cst(0.0); 3];
+    for x in secondary {
+        for k in 0..3 {
+            center[k] = center[k] + x[k] / count;
+        }
+    }
+    // A basis of the plane, (e1, e2, n) right-handed, e1 towards the first corner.
+    let r = sub(*secondary.first()?, center);
+    let r = sub(r, scale(normal, dot(r, normal)));
+    let rl = dot(r, r).sqrt();
+    if !(rl.val() > 0.0) {
+        return None;
+    }
+    let e1 = scale(r, T::cst(1.0) / rl);
+    let e2 = cross(normal, e1);
+    let flat = |x: &Point<T>| -> Flat<T> {
+        let d = sub(*x, center);
+        [dot(d, e1), dot(d, e2)]
+    };
+    let own: Vec<Flat<T>> = secondary.iter().map(flat).collect();
+    let theirs: Vec<Flat<T>> = main.iter().map(flat).collect();
+    let rule = rule();
+    let (own_pieces, their_pieces) = (pieces(&own), pieces(&theirs));
+    for s in &own_pieces {
+        for m in &their_pieces {
+            let overlap = clip(m, s);
+            let Some(o) = overlap.first().copied() else { continue };
+            for w in overlap.windows(2).skip(1) {
+                let (p, q) = (w[0], w[1]);
+                let doubled = cross2(sub2(p, o), sub2(q, o));
+                if !(doubled.val() > 0.0) {
+                    continue;
+                }
+                for (bary, weight) in rule {
+                    let y = [0, 1].map(|k| o[k] + (p[k] - o[k]) * bary[0] + (q[k] - o[k]) * bary[1]);
+                    let (Some(ns), Some(nm)) = (shape_at(&own, y), shape_at(&theirs, y)) else { continue };
+                    f(&Sample {
+                        weight: doubled * (0.5 * weight),
+                        secondary: ns,
+                        main: nm,
+                        xs: interpolate(secondary, &ns),
+                        xm: interpolate(main, &nm),
+                        normal,
+                    });
+                }
+            }
+        }
+    }
+    Some(())
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::dual::{Dual, Hyper, VARS};
+
+    fn square(x0: f64, y0: f64, side: f64, z: f64, up: bool) -> Vec<[f64; 3]> {
+        let c = [[x0, y0, z], [x0 + side, y0, z], [x0 + side, y0 + side, z], [x0, y0 + side, z]];
+        if up { c.to_vec() } else { vec![c[0], c[3], c[2], c[1]] }
+    }
+
+    /// ∫ N_a over the overlap and ∫ N_a g, per secondary corner.
+    fn integrals(s: &[[f64; 3]], m: &[[f64; 3]]) -> ([f64; 4], [f64; 4], f64) {
+        let (mut d, mut g, mut area) = ([0.0; 4], [0.0; 4], 0.0);
+        integrate(s, m, |p| {
+            area += p.weight;
+            for a in 0..4 {
+                d[a] += p.weight * p.secondary[a];
+                g[a] += p.weight * p.secondary[a] * p.gap();
+            }
+        })
+        .unwrap();
+        (d, g, area)
+    }
+
+    #[test]
+    fn a_covered_face_integrates_its_shape_functions_and_gap() {
+        // Secondary face up (normal +z), main face facing it 0.1 above.
+        let s = square(0.0, 0.0, 1.0, 0.0, true);
+        let m = square(-0.5, -0.5, 2.0, 0.1, false);
+        let (d, g, area) = integrals(&s, &m);
+        assert!((area - 1.0).abs() < 1e-14);
+        for a in 0..4 {
+            assert!((d[a] - 0.25).abs() < 1e-14 && (g[a] - 0.025).abs() < 1e-14, "{d:?} {g:?}");
+        }
+        // Half covered: the overlap is the strip x > 0.5.
+        let (d, _, area) = integrals(&s, &square(0.5, -0.5, 2.0, 0.1, false));
+        assert!((area - 0.5).abs() < 1e-14);
+        assert!((d[0] - 0.0625).abs() < 1e-14 && (d[1] - 0.1875).abs() < 1e-14, "{d:?}");
+        // Apart: nothing.
+        let (_, _, area) = integrals(&s, &square(3.0, 0.0, 1.0, 0.1, false));
+        assert_eq!(area, 0.0);
+    }
+
+    /// The sums over a 3 × 3 grid of main faces (triangles in one corner) under secondary face
+    /// `s`, against a tensor-product Gauss rule on the face itself.
+    fn covered(s: &[[f64; 3]; 4]) -> ([f64; 4], [f64; 4], f64, f64) {
+        let mut d = [0.0; 4];
+        let mut area = 0.0;
+        for i in 0..3 {
+            for j in 0..3 {
+                let quad = square(-0.3 + 0.6 * i as f64, -0.3 + 0.6 * j as f64, 0.6, 0.2, false);
+                let faces: Vec<Vec<[f64; 3]>> =
+                    if i == 0 && j == 0 { vec![vec![quad[0], quad[1], quad[2]], vec![quad[0], quad[2], quad[3]]] } else { vec![quad] };
+                for f in faces {
+                    let (da, _, a) = integrals(s, &f);
+                    area += a;
+                    for k in 0..4 {
+                        d[k] += da[k];
+                    }
+                }
+            }
+        }
+        let mut reference = [0.0; 4];
+        let mut exact = 0.0;
+        let g = [-(0.6f64).sqrt(), 0.0, 0.6f64.sqrt()];
+        let w = [5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0];
+        for (a, wa) in g.iter().zip(w) {
+            for (b, wb) in g.iter().zip(w) {
+                let n = [(1.0 - a) * (1.0 - b), (1.0 + a) * (1.0 - b), (1.0 + a) * (1.0 + b), (1.0 - a) * (1.0 + b)].map(|v| v / 4.0);
+                let dxa: [f64; 3] =
+                    std::array::from_fn(|k| 0.25 * (-(1.0 - b) * s[0][k] + (1.0 - b) * s[1][k] + (1.0 + b) * s[2][k] - (1.0 + b) * s[3][k]));
+                let dxb: [f64; 3] =
+                    std::array::from_fn(|k| 0.25 * (-(1.0 - a) * s[0][k] - (1.0 + a) * s[1][k] + (1.0 + a) * s[2][k] + (1.0 - a) * s[3][k]));
+                let j = crate::frame::norm(crate::frame::cross(dxa, dxb));
+                exact += wa * wb * j;
+                for k in 0..4 {
+                    reference[k] += wa * wb * j * n[k];
+                }
+            }
+        }
+        (d, reference, area, exact)
+    }
+
+    #[test]
+    fn non_matching_main_faces_add_up_to_the_covered_integrals() {
+        // A parallelogram: its shape functions are polynomials of the plane, integrated exactly.
+        let (d, reference, area, exact) = covered(&[[0.0, 0.0, 0.0], [1.0, 0.2, 0.0], [1.1, 1.1, 0.0], [0.1, 0.9, 0.0]]);
+        assert!((area - exact).abs() < 1e-13, "{area} {exact}");
+        for k in 0..4 {
+            assert!((d[k] - reference[k]).abs() < 1e-13, "{d:?} {reference:?}");
+        }
+        // A general quadrilateral: its inverse map is not polynomial; the rule is accurate to
+        // its degree.
+        let (d, reference, area, exact) = covered(&[[0.0, 0.0, 0.0], [1.1, 0.1, 0.0], [1.0, 0.9, 0.0], [0.05, 1.2, 0.0]]);
+        assert!((area - exact).abs() < 1e-13, "{area} {exact}");
+        for k in 0..4 {
+            assert!((d[k] - reference[k]).abs() < 1e-5 * reference[k], "{d:?} {reference:?}");
+        }
+    }
+
+    /// ∫ N_a g over a warped, tilted pair, as a function of all 24 coordinates.
+    fn weighted_gap<T: Real>(x: &[T; VARS], a: usize) -> T {
+        let s: Vec<Point<T>> = (0..4).map(|i| [x[3 * i], x[3 * i + 1], x[3 * i + 2]]).collect();
+        let m: Vec<Point<T>> = (0..4).map(|j| [x[12 + 3 * j], x[13 + 3 * j], x[14 + 3 * j]]).collect();
+        let mut g = T::cst(0.0);
+        integrate(&s, &m, |p| g = g + p.weight * p.secondary[a] * p.gap()).unwrap();
+        g
+    }
+
+    #[test]
+    fn the_automatic_derivatives_are_those_of_the_integral() {
+        let base: [f64; VARS] = [
+            0.0, 0.0, 0.0, 1.0, 0.05, 0.02, 1.1, 0.9, -0.03, -0.05, 1.0, 0.01, // secondary
+            -0.2, -0.3, 0.12, -0.25, 1.15, 0.08, 1.2, 1.1, 0.1, 0.9, -0.2, 0.15, // main, facing down
+        ];
+        let a = 2;
+        let duals: [Dual; VARS] = std::array::from_fn(|k| Dual::var(base[k], k));
+        let hypers: [Hyper; VARS] = std::array::from_fn(|k| Hyper::var(base[k], k));
+        let d = weighted_gap(&duals, a);
+        let e = weighted_gap(&hypers, a);
+        let v = weighted_gap(&base, a);
+        assert!((d.v - v).abs() < 1e-14 && (e.v - v).abs() < 1e-14);
+        let h = 1e-6;
+        let shifted = |k: usize, by: f64| {
+            let mut x = base;
+            x[k] += by;
+            x
+        };
+        for k in 0..VARS {
+            let fd = (weighted_gap(&shifted(k, h), a) - weighted_gap(&shifted(k, -h), a)) / (2.0 * h);
+            assert!((d.g[k] - fd).abs() < 1e-7 && (e.g[k] - d.g[k]).abs() < 1e-13, "gradient {k}: {} {fd}", d.g[k]);
+            // The Hessian's row k against differences of the gradient.
+            let (plus, minus) = (shifted(k, h), shifted(k, -h));
+            let gp = weighted_gap(&std::array::from_fn::<Dual, VARS, _>(|i| Dual::var(plus[i], i)), a);
+            let gm = weighted_gap(&std::array::from_fn::<Dual, VARS, _>(|i| Dual::var(minus[i], i)), a);
+            for j in 0..VARS {
+                let fd = (gp.g[j] - gm.g[j]) / (2.0 * h);
+                assert!((e.hessian(k, j) - fd).abs() < 1e-6, "hessian ({k}, {j}): {} {fd}", e.hessian(k, j));
+            }
+        }
+        // A rigid translation changes nothing: the gradient sums to zero per direction.
+        for c in 0..3 {
+            let sum: f64 = (0..8).map(|n| d.g[3 * n + c]).sum();
+            assert!(sum.abs() < 1e-13, "{sum}");
+        }
+    }
+
+    #[test]
+    fn degenerate_faces_give_nothing() {
+        let s = square(0.0, 0.0, 1.0, 0.0, true);
+        let collapsed = [[0.0; 3], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]];
+        assert!(integrate(&collapsed, &s, |_| {}).is_none());
+        assert!(integrate(&s, &s[..2], |_| {}).is_none());
+        let mut count = 0;
+        integrate(&s, &[[f64::NAN, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], |_| count += 1);
+        assert_eq!(count, 0);
+    }
+}