//! Verification of the incompatible-modes quadrilateral (Wilson Q6 in
//! Taylor's QM6 form).
//!
//! # Why the element exists
//!
//! A bilinear `Quad4` interpolates displacement as `u = a + bx + cy + dxy`.
//! Put pure bending into it — the exact field `u = -kxy`, `v = k(x² + νy²)/2`
//! — and the `Quad4` can represent `u` but not `v`, because `v` is quadratic
//! in `x`. What it produces instead is a field whose shear strain
//! `γ = ∂u/∂y + ∂v/∂x` is non-zero everywhere except along the element
//! midlines, even though the exact bending field has `γ ≡ 0`. That spurious
//! shear absorbs strain energy the beam should not store, so the element is
//! far too stiff in bending — *shear locking*. The error grows with the square
//! of the element's length-to-height ratio, so it is worst on exactly the
//! coarse meshes one wants to use.
//!
//! Wilson's remedy is to add two internal *incompatible* modes,
//! `N₅ = 1 - ξ²` and `N₆ = 1 - η²`, each with an `x` and a `y` amplitude.
//! `N₅` supplies the missing `x²` term, so the element can now bend without
//! shearing. The four amplitudes are internal to the element and are removed
//! by static condensation before assembly, so the element still presents an
//! 8×8 stiffness to the mesh.
//!
//! # Why QM6 rather than Q6
//!
//! The added modes are *incompatible*: they do not vanish on the element
//! boundary, so displacement continuity between neighbours is broken and
//! convergence is no longer automatic. It has to be earned by passing the
//! patch test, which for the internal block requires
//!
//! ```text
//!     ∫ B_a dΩ = 0
//! ```
//!
//! Wilson's original Q6 evaluates `B_a` with the Jacobian at each quadrature
//! point. For a general quadrilateral `|J|` varies bilinearly over the
//! element, that integral is not zero, and Q6 fails the patch test on any
//! element that is not a parallelogram. Taylor, Beresford and Wilson's QM6
//! evaluates the incompatible-mode Jacobian **at the element centre** and
//! scales by `|J₀|/|J|`; the weight `|J|` then cancels and the integral
//! reduces to `|J₀| ∫ B_a⁰ dΩ`, whose integrand is linear in `ξ` and `η` and
//! therefore integrates to zero exactly on the symmetric 2×2 Gauss rule. That
//! is the whole content of the correction, and it is what the distorted-mesh
//! tests below check.
//!
//! # What each test rules out
//!
//! Every test here is written against a closed form or an algebraic invariant
//! that a plausible-but-wrong implementation violates:
//!
//! * A stiffness of zeros, or one that ignores the internal block entirely,
//!   is caught by [`qm6_is_strictly_softer_than_quad4`] — the difference
//!   `K_Q4 - K_QM6` must be positive semi-definite *and* non-zero.
//! * A sign error in the condensation (`Kuu + Kua Kaa⁻¹ Kau`) leaves the rigid
//!   body modes intact and the matrix symmetric, so only the semi-definiteness
//!   of that same difference, and the cantilever, catch it.
//! * Evaluating the incompatible block at the quadrature point instead of the
//!   centre — i.e. Q6 rather than QM6 — passes every rigid-body and symmetry
//!   check and still gives a good cantilever answer. Only the distorted-mesh
//!   patch tests catch it.

use nalgebra::{DMatrix, DVector, Vector3};
use rtx_fea::elements::{
    ElementMatrixComputer, FiniteElement, NaturalCoords, StandardFiniteElement,
};
use rtx_fea::mesh::ElementType;

/// Which quadrilateral formulation to build the element stiffness with.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum Formulation {
    /// The existing bilinear element, 2×2 Gauss.
    Quad4,
    /// Incompatible modes with static condensation, Taylor's QM6 form.
    Qm6,
}

impl Formulation {
    fn name(self) -> &'static str {
        match self {
            Formulation::Quad4 => "Quad4",
            Formulation::Qm6 => "QM6",
        }
    }
}

/// Build the 8×8 element stiffness for a four-node quadrilateral.
fn element_stiffness(coords: &[Vector3<f64>], e: f64, nu: f64, form: Formulation) -> DMatrix<f64> {
    let fe = StandardFiniteElement::new(ElementType::Quad4, coords.to_vec());
    let result = match form {
        Formulation::Quad4 => {
            ElementMatrixComputer::compute_stiffness_matrix(&fe, coords, e, nu, None)
        }
        Formulation::Qm6 => {
            ElementMatrixComputer::compute_stiffness_matrix_incompatible(&fe, coords, e, nu, None)
        }
    };
    result
        .unwrap_or_else(|err| panic!("{} stiffness failed: {err}", form.name()))
        .matrix
}

fn node(x: f64, y: f64) -> Vector3<f64> {
    Vector3::new(x, y, 0.0)
}

/// A quadrilateral that is neither a rectangle nor a parallelogram.
///
/// Both pairs of opposite edges are non-parallel, so `|J|` genuinely varies
/// over the element and the Q6/QM6 distinction has something to bite on.
fn distorted_element() -> Vec<Vector3<f64>> {
    vec![
        node(0.0, 0.0),
        node(2.0, 0.0),
        node(2.2, 1.0),
        node(0.6, 1.2),
    ]
}

/// Nodal values of the linear field `u = a₀ + a₁x + a₂y`, `v = b₀ + b₁x + b₂y`,
/// interleaved `[u₀, v₀, u₁, v₁, …]`.
fn linear_field(coords: &[Vector3<f64>], a: [f64; 3], b: [f64; 3]) -> DVector<f64> {
    let mut u = DVector::zeros(2 * coords.len());
    for (i, c) in coords.iter().enumerate() {
        u[2 * i] = a[0] + a[1] * c.x + a[2] * c.y;
        u[2 * i + 1] = b[0] + b[1] * c.x + b[2] * c.y;
    }
    u
}

// ---------------------------------------------------------------------------
// Premise check
// ---------------------------------------------------------------------------

/// The physical shape-function derivatives must be complete to first order on
/// a distorted element.
///
/// `Σ ∂Nᵢ/∂x = 0`, `Σ (∂Nᵢ/∂x) xᵢ = 1`, `Σ (∂Nᵢ/∂x) yᵢ = 0`, and the
/// transposes. These are the identities that make constant strain
/// representable, and every claim below rests on them. If the Jacobian
/// transform were transposed — a mistake that is invisible on rectangles,
/// where `J` is diagonal — the cross terms here would be non-zero and no patch
/// test could pass.
#[test]
fn physical_derivatives_are_first_order_complete_on_a_distorted_element() {
    let coords = distorted_element();
    let fe = StandardFiniteElement::new(ElementType::Quad4, coords.clone());

    for &(xi, eta) in &[(0.0, 0.0), (-0.6, 0.3), (0.577, -0.577)] {
        let nat = NaturalCoords::new_2d(xi, eta);
        let shape = fe.shape_functions(&nat).unwrap();
        let jac = fe.jacobian(&nat, &coords).unwrap();
        let d = jac.transform_derivatives(&shape.derivatives).unwrap();

        let mut sum_dx = 0.0;
        let mut sum_dy = 0.0;
        let mut dx_x = 0.0;
        let mut dx_y = 0.0;
        let mut dy_x = 0.0;
        let mut dy_y = 0.0;
        for i in 0..4 {
            sum_dx += d[(i, 0)];
            sum_dy += d[(i, 1)];
            dx_x += d[(i, 0)] * coords[i].x;
            dx_y += d[(i, 0)] * coords[i].y;
            dy_x += d[(i, 1)] * coords[i].x;
            dy_y += d[(i, 1)] * coords[i].y;
        }

        let tol = 1e-12;
        assert!(
            sum_dx.abs() < tol,
            "at ({xi}, {eta}): sum dN/dx = {sum_dx:e}"
        );
        assert!(
            sum_dy.abs() < tol,
            "at ({xi}, {eta}): sum dN/dy = {sum_dy:e}"
        );
        assert!(
            (dx_x - 1.0).abs() < tol,
            "at ({xi}, {eta}): sum (dN/dx) x = {dx_x}, expected 1"
        );
        assert!(
            (dy_y - 1.0).abs() < tol,
            "at ({xi}, {eta}): sum (dN/dy) y = {dy_y}, expected 1"
        );
        assert!(
            dx_y.abs() < tol,
            "at ({xi}, {eta}): sum (dN/dx) y = {dx_y:e}, expected 0"
        );
        assert!(
            dy_x.abs() < tol,
            "at ({xi}, {eta}): sum (dN/dy) x = {dy_x:e}, expected 0"
        );
    }
}

// ---------------------------------------------------------------------------
// 1. Patch test
// ---------------------------------------------------------------------------

/// Single-element patch test in its sharpest algebraic form.
///
/// Static condensation gives `K = Kuu - Kua Kaa⁻¹ Kau`, so for any nodal
/// vector `u`
///
/// ```text
///     K u = Kuu u - Kua Kaa⁻¹ (Kau u)
/// ```
///
/// and `Kuu` *is* the ordinary `Quad4` stiffness — same `B`, same quadrature.
/// So `K_QM6 u = K_Q4 u` holds **exactly** whenever `Kau u = 0`, and `Kau u`
/// is zero for a linear displacement field precisely when `∫ B_a dΩ = 0`,
/// which is the patch-test condition. Comparing the two force vectors
/// therefore tests the internal block directly, with no solver in the way.
///
/// Wilson's Q6 fails this on a non-parallelogram; QM6 passes it for any shape.
/// Three independent linear fields are used so the check covers all three
/// constant-strain states.
#[test]
fn qm6_patch_test_single_distorted_element() {
    let coords = distorted_element();
    let e = 1000.0;
    let nu = 0.3;

    let k_q4 = element_stiffness(&coords, e, nu, Formulation::Quad4);
    let k_qm6 = element_stiffness(&coords, e, nu, Formulation::Qm6);

    // Uniform x-stretch, uniform y-stretch, pure shear, and a general mix.
    let fields: [([f64; 3], [f64; 3]); 4] = [
        ([0.0, 1.0, 0.0], [0.0, 0.0, 0.0]),
        ([0.0, 0.0, 0.0], [0.0, 0.0, 1.0]),
        ([0.0, 0.0, 1.0], [0.0, 1.0, 0.0]),
        ([0.003, 0.011, -0.007], [-0.002, 0.005, 0.013]),
    ];

    for (a, b) in fields {
        let u = linear_field(&coords, a, b);
        let f_q4 = &k_q4 * &u;
        let f_qm6 = &k_qm6 * &u;
        let scale = f_q4.amax().max(1e-30);
        let diff = (&f_qm6 - &f_q4).amax();
        assert!(
            diff / scale < 1e-11,
            "linear field u={a:?} v={b:?}: QM6 nodal forces differ from the \
             constant-stress forces by {diff:e} (relative {:e}). The internal \
             modes were activated by a linear field, so the patch test fails — \
             this is what Q6 does on a non-parallelogram.",
            diff / scale
        );

        // And those forces really are a self-equilibrated constant-stress
        // state: zero resultant force and zero resultant moment.
        let fx: f64 = (0..4).map(|i| f_qm6[2 * i]).sum();
        let fy: f64 = (0..4).map(|i| f_qm6[2 * i + 1]).sum();
        let m: f64 = (0..4)
            .map(|i| coords[i].x * f_qm6[2 * i + 1] - coords[i].y * f_qm6[2 * i])
            .sum();
        assert!(fx.abs() / scale < 1e-12, "resultant Fx = {fx:e}");
        assert!(fy.abs() / scale < 1e-12, "resultant Fy = {fy:e}");
        assert!(m.abs() / scale < 1e-12, "resultant moment = {m:e}");
    }
}

/// The classical multi-element patch test: a patch of four distorted
/// quadrilaterals around one interior node.
///
/// The exact linear field is imposed on all eight boundary nodes and the two
/// interior degrees of freedom are solved for. A convergent element must
/// return the interior node to the exact field — no element in the patch may
/// generate any residual force under a constant-stress state. The interior
/// node sits well off centre so every one of the four elements is a distinct
/// non-parallelogram.
#[test]
fn qm6_patch_test_four_distorted_elements() {
    let model = patch_model();
    let e = 1000.0;
    let nu = 0.3;
    let a = [0.003, 0.011, -0.007];
    let b = [-0.002, 0.005, 0.013];

    let interior = model.nodes.len() - 1;
    let exact = linear_field(&model.nodes, a, b);

    for form in [Formulation::Quad4, Formulation::Qm6] {
        let k = model.assemble(e, nu, form);
        let f = DVector::zeros(k.nrows());
        let prescribed: Vec<(usize, f64)> = (0..interior)
            .flat_map(|n| [(2 * n, exact[2 * n]), (2 * n + 1, exact[2 * n + 1])])
            .collect();
        let u = solve_with_constraints(&k, &f, &prescribed);

        let du = (u[2 * interior] - exact[2 * interior]).abs();
        let dv = (u[2 * interior + 1] - exact[2 * interior + 1]).abs();
        let scale = exact.amax();
        assert!(
            du / scale < 1e-11 && dv / scale < 1e-11,
            "{}: interior node came out at ({:.15e}, {:.15e}), exact is \
             ({:.15e}, {:.15e}) — errors {du:e}, {dv:e}. The patch test fails.",
            form.name(),
            u[2 * interior],
            u[2 * interior + 1],
            exact[2 * interior],
            exact[2 * interior + 1],
        );
    }
}

// ---------------------------------------------------------------------------
// 2. Rigid body modes and the condensation itself
// ---------------------------------------------------------------------------

/// The condensed stiffness must keep exactly three zero eigenvalues in 2-D.
///
/// Two translations and one infinitesimal rotation — no more, no fewer. Fewer
/// means the condensation has locked out a rigid motion; more means it has
/// introduced a spurious zero-energy (hourglass) mode, which is the failure
/// mode of under-integration, the *other* common cure for locking.
#[test]
fn qm6_has_exactly_three_rigid_body_modes() {
    let coords = distorted_element();
    let e = 1000.0;
    let nu = 0.3;
    let k = element_stiffness(&coords, e, nu, Formulation::Qm6);

    assert_eq!(k.nrows(), 8, "condensed stiffness must stay 8x8");
    assert_eq!(k.ncols(), 8);

    let asym = (&k - k.transpose()).amax();
    assert!(
        asym / k.amax() < 1e-12,
        "condensed stiffness is not symmetric: max |K - K^T| = {asym:e}"
    );

    let mut eig: Vec<f64> = k.clone().symmetric_eigenvalues().iter().copied().collect();
    eig.sort_by(|p, q| p.partial_cmp(q).unwrap());
    let largest = eig[7];

    let zeros = eig.iter().filter(|v| v.abs() < 1e-10 * largest).count();
    assert_eq!(
        zeros,
        3,
        "expected exactly 3 rigid body modes, found {zeros}. Eigenvalues \
         (scaled by the largest, {largest:e}): {:?}",
        eig.iter().map(|v| v / largest).collect::<Vec<_>>()
    );
    assert!(
        eig[3] > 1e-6 * largest,
        "the first deformation mode has eigenvalue {:e}, only {:e} of the \
         largest — that is a spurious zero-energy mode",
        eig[3],
        eig[3] / largest
    );
    assert!(
        eig[0] > -1e-10 * largest,
        "stiffness is not positive semi-definite: smallest eigenvalue {:e}",
        eig[0]
    );
}

/// Rigid translation and infinitesimal rigid rotation must store no energy.
///
/// A zero matrix also passes the eigenvalue count above (it has eight zero
/// eigenvalues, not three, so in fact it does not — but a matrix that is
/// merely *scaled wrongly* would). This check pins the energy of the rigid
/// modes against the energy of a real deformation of the same amplitude, so it
/// is a ratio, not an absolute, and cannot be satisfied by a small matrix.
#[test]
fn qm6_stores_no_energy_in_rigid_body_motion() {
    let coords = distorted_element();
    let k = element_stiffness(&coords, 1000.0, 0.3, Formulation::Qm6);

    let translate_x = linear_field(&coords, [1.0, 0.0, 0.0], [0.0, 0.0, 0.0]);
    let translate_y = linear_field(&coords, [0.0, 0.0, 0.0], [1.0, 0.0, 0.0]);
    // u = -theta y, v = theta x with theta = 1.
    let rotate = linear_field(&coords, [0.0, 0.0, -1.0], [0.0, 1.0, 0.0]);
    // A genuine deformation of comparable amplitude, for the scale.
    let stretch = linear_field(&coords, [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]);

    let energy = |u: &DVector<f64>| (u.transpose() * &k * u)[(0, 0)];
    let reference = energy(&stretch);
    assert!(reference > 0.0, "reference stretch energy is {reference:e}");

    for (name, u) in [
        ("translation in x", &translate_x),
        ("translation in y", &translate_y),
        ("infinitesimal rotation", &rotate),
    ] {
        let w = energy(u);
        assert!(
            w.abs() / reference < 1e-12,
            "{name} stores energy {w:e}, which is {:e} of the stretch energy \
             {reference:e}",
            w.abs() / reference
        );
    }
}

/// Static condensation can only remove stiffness, never add it — and here it
/// must actually remove some.
///
/// The condensed energy is the minimum of the augmented energy over the
/// internal amplitudes, and `a = 0` is always admissible, so
/// `uᵀK_QM6 u ≤ uᵀK_Q4 u` for every `u`. Equivalently `K_Q4 - K_QM6 =
/// Kua Kaa⁻¹ Kau` is positive semi-definite. This is the invariant that pins
/// the *sign* and the *form* of the condensation: flipping the sign makes the
/// difference negative semi-definite, and dropping the internal block
/// altogether makes it exactly zero. Both are caught here and by nothing else
/// in this file except the cantilever.
#[test]
fn qm6_is_strictly_softer_than_quad4() {
    let coords = distorted_element();
    let e = 1000.0;
    let nu = 0.3;

    let k_q4 = element_stiffness(&coords, e, nu, Formulation::Quad4);
    let k_qm6 = element_stiffness(&coords, e, nu, Formulation::Qm6);
    let diff = &k_q4 - &k_qm6;

    let scale = k_q4.amax();
    let mut eig: Vec<f64> = diff
        .clone()
        .symmetric_eigenvalues()
        .iter()
        .copied()
        .collect();
    eig.sort_by(|p, q| p.partial_cmp(q).unwrap());

    assert!(
        eig[0] > -1e-10 * scale,
        "K_Q4 - K_QM6 has eigenvalue {:e} (scale {scale:e}); condensation made \
         the element STIFFER, so the sign of Kua Kaa^-1 Kau is wrong",
        eig[0]
    );
    assert!(
        eig[7] > 1e-3 * scale,
        "K_Q4 - K_QM6 is essentially zero (largest eigenvalue {:e} against a \
         stiffness scale of {scale:e}); the internal modes are not doing \
         anything",
        eig[7]
    );
    // Rank of the difference is at most 4 — that is the size of the internal
    // block. More than four non-negligible eigenvalues means the partition is
    // wrong.
    let nonzero = eig.iter().filter(|v| **v > 1e-10 * scale).count();
    assert!(
        nonzero <= 4,
        "K_Q4 - K_QM6 has rank {nonzero}, but the internal block has only 4 \
         degrees of freedom: {:?}",
        eig
    );
}

// ---------------------------------------------------------------------------
// 3. The point of the exercise: cantilever bending
// ---------------------------------------------------------------------------

/// Slender cantilever under an end load, against the Euler-Bernoulli closed
/// form.
///
/// ```text
///     delta = P L³ / (3 E I),   I = h³ / 12  (unit thickness)
/// ```
///
/// With `L = 10`, `h = 1`, `E = 1e7`, `P = 1`: `I = 1/12` and
/// `delta = 1000 / (3 · 1e7 / 12) = 4.0e-4`.
///
/// Poisson's ratio is set to zero deliberately. Beam theory has no transverse
/// contraction, so any non-zero `ν` puts a modelling difference between the
/// closed form and the two-dimensional answer that has nothing to do with the
/// element formulation, and would muddy exactly the comparison being made.
///
/// The remaining known difference is transverse shear, which Euler-Bernoulli
/// omits and plane-stress elasticity includes:
///
/// ```text
///     delta_shear = 6 P L / (5 G A) = 6 · 1 · 10 / (5 · 5e6 · 1) = 2.4e-6
/// ```
///
/// — 0.6% of the bending deflection at `L/h = 10`. So the correct 2-D answer
/// sits slightly *above* Euler-Bernoulli, and the clamped end (which restrains
/// the transverse contraction and the warping that the exact solution has) pulls
/// it back down a little. A converged answer within a couple of percent of
/// `4.0e-4` is what the physics allows; anything far below it is locking.
#[test]
fn cantilever_tip_deflection_quad4_versus_qm6() {
    const L: f64 = 10.0;
    const H: f64 = 1.0;
    const E: f64 = 1.0e7;
    const NU: f64 = 0.0;
    const P: f64 = 1.0;

    let inertia = H * H * H / 12.0;
    let euler_bernoulli = P * L * L * L / (3.0 * E * inertia);
    let shear = 6.0 * P * L / (5.0 * (E / 2.0) * H);
    let timoshenko = euler_bernoulli + shear;

    println!("\ncantilever L={L} h={H} E={E:e} nu={NU} P={P}");
    println!("  Euler-Bernoulli  delta = {euler_bernoulli:.6e}");
    println!(
        "  Timoshenko       delta = {timoshenko:.6e} (shear adds {:.2}%)",
        100.0 * shear / euler_bernoulli
    );

    // The coarse mesh the comparison is made on: ten square elements, one
    // through the depth. This is the regime where locking is at its worst and
    // where a practitioner would actually be.
    let coarse = (10usize, 1usize);
    let mut coarse_results = [0.0f64; 2];

    println!("\n  mesh     Quad4 tip      err        QM6 tip        err");
    for &(nx, ny) in &[coarse, (20, 2), (40, 4), (80, 8)] {
        let model = cantilever_model(L, H, nx, ny);
        let mut row = [0.0f64; 2];
        for (slot, form) in [Formulation::Quad4, Formulation::Qm6].iter().enumerate() {
            let d = model.cantilever_tip_deflection(E, NU, *form, P, L);
            row[slot] = d;
        }
        if (nx, ny) == coarse {
            coarse_results = row;
        }
        println!(
            "  {:>2}x{:<2}  {:.6e}  {:+7.2}%   {:.6e}  {:+7.2}%",
            nx,
            ny,
            row[0],
            100.0 * (row[0] - euler_bernoulli) / euler_bernoulli,
            row[1],
            100.0 * (row[1] - euler_bernoulli) / euler_bernoulli,
        );
    }

    let quad4 = coarse_results[0];
    let qm6 = coarse_results[1];
    let err_quad4 = (quad4 - euler_bernoulli) / euler_bernoulli;
    let err_qm6 = (qm6 - euler_bernoulli) / euler_bernoulli;

    println!(
        "\n  coarse mesh {}x{}: Quad4 {:.6e} ({:+.2}%), QM6 {:.6e} ({:+.2}%)",
        coarse.0,
        coarse.1,
        quad4,
        100.0 * err_quad4,
        qm6,
        100.0 * err_qm6
    );

    // Quad4 must be dramatically too stiff — that is the disease.
    assert!(
        err_quad4 < -0.25,
        "Quad4 tip deflection {quad4:.6e} is {:+.2}% off beam theory \
         {euler_bernoulli:.6e}; shear locking on a {}x{} mesh should be far \
         worse than that, so either the mesh or the load is not what this test \
         thinks it is",
        100.0 * err_quad4,
        coarse.0,
        coarse.1
    );

    // QM6 must land on beam theory. Tolerance set from the physics above:
    // Timoshenko shear adds 0.6%, the clamped end takes back a comparable
    // amount, so 5% is generous but still an order tighter than the locking
    // error it has to beat.
    assert!(
        err_qm6.abs() < 0.05,
        "QM6 tip deflection {qm6:.6e} is {:+.2}% off beam theory \
         {euler_bernoulli:.6e}; the incompatible modes are not curing the \
         locking",
        100.0 * err_qm6
    );

    assert!(
        err_qm6.abs() * 5.0 < err_quad4.abs(),
        "QM6 error {:+.2}% is not decisively better than Quad4's {:+.2}%",
        100.0 * err_qm6,
        100.0 * err_quad4
    );
}

// ---------------------------------------------------------------------------
// Test harness: assembly and a direct solve
// ---------------------------------------------------------------------------

/// A tiny quadrilateral mesh with a dense assembler, kept local to this test
/// so that the Quad4 and QM6 numbers differ only in the element formulation.
struct Model {
    nodes: Vec<Vector3<f64>>,
    elems: Vec<[usize; 4]>,
}

impl Model {
    /// Assemble the global stiffness with DOFs interleaved per node,
    /// `[u₀, v₀, u₁, v₁, …]`, matching the element DOF layout.
    fn assemble(&self, e: f64, nu: f64, form: Formulation) -> DMatrix<f64> {
        let n = self.nodes.len();
        let mut k = DMatrix::zeros(2 * n, 2 * n);
        for conn in &self.elems {
            let coords: Vec<Vector3<f64>> = conn.iter().map(|&i| self.nodes[i]).collect();
            let ke = element_stiffness(&coords, e, nu, form);
            for a in 0..4 {
                for b in 0..4 {
                    for i in 0..2 {
                        for j in 0..2 {
                            k[(2 * conn[a] + i, 2 * conn[b] + j)] += ke[(2 * a + i, 2 * b + j)];
                        }
                    }
                }
            }
        }
        k
    }

    /// Clamp `x = 0`, hang a total load `p` off the tip, return the mean
    /// vertical displacement of the tip edge.
    fn cantilever_tip_deflection(
        &self,
        e: f64,
        nu: f64,
        form: Formulation,
        p: f64,
        length: f64,
    ) -> f64 {
        let k = self.assemble(e, nu, form);
        let mut f = DVector::zeros(k.nrows());

        let tip: Vec<usize> = (0..self.nodes.len())
            .filter(|&i| (self.nodes[i].x - length).abs() < 1e-12)
            .collect();
        let share = -p / tip.len() as f64;
        for &i in &tip {
            f[2 * i + 1] += share;
        }

        let prescribed: Vec<(usize, f64)> = (0..self.nodes.len())
            .filter(|&i| self.nodes[i].x.abs() < 1e-12)
            .flat_map(|i| [(2 * i, 0.0), (2 * i + 1, 0.0)])
            .collect();

        let u = solve_with_constraints(&k, &f, &prescribed);
        let mean: f64 = tip.iter().map(|&i| u[2 * i + 1]).sum::<f64>() / tip.len() as f64;
        -mean
    }
}

/// `nx` by `ny` uniform mesh of the rectangle `[0, length] × [0, height]`,
/// nodes ordered counter-clockwise within each element.
fn cantilever_model(length: f64, height: f64, nx: usize, ny: usize) -> Model {
    let mut nodes = Vec::with_capacity((nx + 1) * (ny + 1));
    for j in 0..=ny {
        for i in 0..=nx {
            nodes.push(node(
                length * i as f64 / nx as f64,
                height * j as f64 / ny as f64,
            ));
        }
    }
    let idx = |i: usize, j: usize| j * (nx + 1) + i;
    let mut elems = Vec::with_capacity(nx * ny);
    for j in 0..ny {
        for i in 0..nx {
            elems.push([idx(i, j), idx(i + 1, j), idx(i + 1, j + 1), idx(i, j + 1)]);
        }
    }
    Model { nodes, elems }
}

/// Four distorted quadrilaterals filling `[0,2]²`, sharing one interior node
/// placed off centre. The interior node is the last entry so the boundary
/// nodes are `0..8`.
fn patch_model() -> Model {
    let nodes = vec![
        node(0.0, 0.0), // 0  corner
        node(1.0, 0.0), // 1  edge
        node(2.0, 0.0), // 2  corner
        node(2.0, 1.0), // 3  edge
        node(2.0, 2.0), // 4  corner
        node(1.0, 2.0), // 5  edge
        node(0.0, 2.0), // 6  corner
        node(0.0, 1.0), // 7  edge
        node(1.3, 0.9), // 8  interior, deliberately off centre
    ];
    let elems = vec![[0, 1, 8, 7], [1, 2, 3, 8], [8, 3, 4, 5], [7, 8, 5, 6]];
    Model { nodes, elems }
}

/// Solve `K u = f` with some degrees of freedom prescribed, by reduction.
fn solve_with_constraints(
    k: &DMatrix<f64>,
    f: &DVector<f64>,
    prescribed: &[(usize, f64)],
) -> DVector<f64> {
    let n = k.nrows();
    let mut u = DVector::zeros(n);
    let mut fixed = vec![false; n];
    for &(dof, value) in prescribed {
        fixed[dof] = true;
        u[dof] = value;
    }
    let free: Vec<usize> = (0..n).filter(|&d| !fixed[d]).collect();

    let mut kff = DMatrix::zeros(free.len(), free.len());
    let mut rhs = DVector::zeros(free.len());
    for (a, &da) in free.iter().enumerate() {
        let mut r = f[da];
        for d in 0..n {
            if fixed[d] {
                r -= k[(da, d)] * u[d];
            }
        }
        rhs[a] = r;
        for (b, &db) in free.iter().enumerate() {
            kff[(a, b)] = k[(da, db)];
        }
    }

    let sol = kff
        .lu()
        .solve(&rhs)
        .expect("reduced stiffness must be non-singular");
    for (a, &da) in free.iter().enumerate() {
        u[da] = sol[a];
    }
    u
}