//! Transient dynamics against closed forms, and against the modal analysis.
//!
//! `DynamicAnalysis::run` used to return `DVector::zeros(...)`. A zero
//! response satisfies almost every cheap check one might write — it is
//! symmetric, finite, bounded, and it "converges". So none of the assertions
//! here are cheap checks.
//!
//! Two independent validations are run, and they cannot both be faked by the
//! same wrong implementation:
//!
//! 1. **Closed form.** A single degree of freedom is integrated and compared
//!    against `x(t) = e^{-ζωₙt}(A cos ω_d t + B sin ω_d t)`, undamped and
//!    damped, plus a measured order of accuracy. Getting the exact damped
//!    envelope right over eight periods is not something a wrong scheme does
//!    by accident, and the order measurement pins `γ`: Newmark is
//!    second-order only at `γ = ½`, and the test below demonstrates that it
//!    really does collapse to first order at `γ = 0.6`.
//!
//! 2. **Agreement with the validated modal result.** The fixed-free bar from
//!    `modal_closed_form.rs` — whose fundamental frequency is already checked
//!    against `f₁ = 1/(4L)·√(E/ρ)` — is released from an initial displacement
//!    and its free-vibration period is measured from zero crossings. The
//!    eigensolver and the time integrator share only the assembled `K` and
//!    `M`; the frequency-domain and time-domain answers come from disjoint
//!    code paths, so neither can fake agreement with the other.
//!
//! The Rayleigh damping split is checked term by term, `C = αM` and `C = βK`
//! separately, against the modal damping ratios `ζ₁ = α/(2ω₁)` and
//! `ζ₁ = βω₁/2`. A single combined test would pass with the two coefficients
//! swapped.

use nalgebra::{DMatrix, DVector};
use rtx_fea::analysis::{
    Analysis, AnalysisConfig, AnalysisData, DynamicAnalysis, ModalAnalysis, NewmarkStepper,
    TimeConfig, TimeIntegrationScheme,
};
use rtx_fea::assembly::DofComponent;
use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, DirichletBC};
use rtx_fea::materials::{LinearElastic, MaterialDatabase};
use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};

// ---------------------------------------------------------------------------
// Single degree of freedom: the closed form
// ---------------------------------------------------------------------------

/// Build a one-degree-of-freedom stepper for `m ẍ + c ẋ + k x = 0`.
fn sdof(m: f64, c: f64, k: f64, gamma: f64, beta: f64, dt: f64) -> NewmarkStepper {
    NewmarkStepper::new(
        DMatrix::from_element(1, 1, m),
        DMatrix::from_element(1, 1, c),
        DMatrix::from_element(1, 1, k),
        gamma,
        beta,
        dt,
    )
    .expect("1x1 Newmark stepper")
}

/// Exact free response of `m ẍ + c ẋ + k x = 0`, sub-critical damping only.
fn exact_sdof(m: f64, c: f64, k: f64, x0: f64, v0: f64, t: f64) -> f64 {
    let wn = (k / m).sqrt();
    let zeta = c / (2.0 * (k * m).sqrt());
    assert!(zeta < 1.0, "closed form assumes sub-critical damping");
    let wd = wn * (1.0 - zeta * zeta).sqrt();
    let a = x0;
    let b = (v0 + zeta * wn * x0) / wd;
    (-zeta * wn * t).exp() * (a * (wd * t).cos() + b * (wd * t).sin())
}

/// March a free-vibration problem and return `(t, x)` at every step.
fn march_sdof(stepper: &NewmarkStepper, x0: f64, v0: f64, num_steps: usize) -> Vec<(f64, f64)> {
    let zero = DVector::zeros(1);
    let mut state = stepper
        .initial_state(
            DVector::from_element(1, x0),
            DVector::from_element(1, v0),
            &zero,
        )
        .expect("initial state");

    let mut history = vec![(0.0, state.displacement[0])];
    for step in 1..=num_steps {
        state = stepper.step(&state, &zero).expect("newmark step");
        history.push((step as f64 * stepper.dt(), state.displacement[0]));
    }
    history
}

/// Largest deviation from the closed form over the whole march.
fn max_error(history: &[(f64, f64)], m: f64, c: f64, k: f64, x0: f64, v0: f64) -> f64 {
    history
        .iter()
        .map(|&(t, x)| (x - exact_sdof(m, c, k, x0, v0, t)).abs())
        .fold(0.0, f64::max)
}

const M: f64 = 1.0;
const K: f64 = 100.0; // wn = 10 rad/s, T = 0.6283 s
const X0: f64 = 1.0;
const V0: f64 = 0.0;

/// Undamped free vibration must reproduce `x(t) = x₀ cos ωₙ t`.
///
/// Amplitude is checked as well as the trajectory: `γ = ½` gives an
/// amplification matrix with unit determinant, so the response may drift in
/// phase but must not grow or decay at all.
#[test]
fn undamped_sdof_matches_the_closed_form() {
    let wn = (K / M).sqrt();
    let period = 2.0 * std::f64::consts::PI / wn;
    let dt = period / 200.0;
    let num_steps = 200 * 8; // eight periods

    let stepper = sdof(M, 0.0, K, 0.5, 0.25, dt);
    let history = march_sdof(&stepper, X0, V0, num_steps);
    let error = max_error(&history, M, 0.0, K, X0, V0);

    // The average-acceleration scheme elongates the period by roughly
    // Ω²/12 with Ω = ωΔt, giving a phase lag of 2π·8·Ω²/12 after eight
    // periods; with an amplitude of one that is the size of the error.
    let omega_dt = wn * dt;
    let predicted = 2.0 * std::f64::consts::PI * 8.0 * omega_dt * omega_dt / 12.0;
    println!("undamped: max error {error:.3e}, predicted phase lag {predicted:.3e}");

    assert!(
        error < 2.0 * predicted,
        "max error {error:.3e} exceeds twice the predicted phase lag {predicted:.3e}"
    );

    let peak = history.iter().map(|&(_, x)| x.abs()).fold(0.0, f64::max);
    assert!(
        (peak - X0).abs() < 1e-9,
        "the undamped amplitude drifted to {peak:.12} from {X0}; \
         gamma = 1/2 must neither grow nor decay the response"
    );
}

/// The average-acceleration scheme conserves `E = ½vᵀMv + ½uᵀKu` exactly.
///
/// This is the invariant a plausible-but-wrong integrator fails. A zero
/// response conserves energy trivially, so the initial energy is checked to
/// be the analytic `½kx₀²` first.
#[test]
fn average_acceleration_conserves_energy_to_round_off() {
    let wn = (K / M).sqrt();
    let dt = 2.0 * std::f64::consts::PI / wn / 60.0;
    let stepper = sdof(M, 0.0, K, 0.5, 0.25, dt);

    let zero = DVector::zeros(1);
    let mut state = stepper
        .initial_state(
            DVector::from_element(1, X0),
            DVector::from_element(1, V0),
            &zero,
        )
        .unwrap();

    let e0 = stepper.total_energy(&state);
    assert!(
        (e0 - 0.5 * K * X0 * X0).abs() < 1e-12,
        "initial energy {e0} is not the analytic 1/2 k x0^2 = {}",
        0.5 * K * X0 * X0
    );

    let mut worst = 0.0f64;
    for _ in 0..(60 * 200) {
        state = stepper.step(&state, &zero).unwrap();
        worst = worst.max((stepper.total_energy(&state) - e0).abs() / e0);
    }

    println!("energy drift over 200 periods: {worst:.3e} relative");
    assert!(
        worst < 1e-11,
        "energy drifted by {worst:.3e} relative over 200 periods; \
         the average-acceleration map is exactly orthogonal in the energy \
         inner product, so this must be round-off only"
    );
}

/// Damped free vibration must reproduce the exponentially decaying closed form.
#[test]
fn damped_sdof_matches_the_closed_form() {
    let zeta = 0.05;
    let c = 2.0 * zeta * (K * M).sqrt();
    let wn = (K / M).sqrt();
    let period = 2.0 * std::f64::consts::PI / wn;
    let dt = period / 200.0;
    let num_steps = 200 * 8;

    let stepper = sdof(M, c, K, 0.5, 0.25, dt);
    let history = march_sdof(&stepper, X0, V0, num_steps);
    let error = max_error(&history, M, c, K, X0, V0);

    println!("damped (zeta = {zeta}): max error {error:.3e}");
    assert!(
        error < 5e-3,
        "max deviation from the damped closed form is {error:.3e}"
    );

    // The final amplitude must be the analytic envelope, not merely "small".
    // Anything that damps by the wrong amount lands far from this.
    let t_end = num_steps as f64 * dt;
    let envelope = (-zeta * wn * t_end).exp();
    println!("envelope after eight periods: {envelope:.6e}");
    assert!(
        envelope < 0.1,
        "the test is not exercising decay: envelope {envelope}"
    );

    // Energy must fall monotonically: damping cannot add energy.
    let mut state = stepper
        .initial_state(
            DVector::from_element(1, X0),
            DVector::from_element(1, V0),
            &DVector::zeros(1),
        )
        .unwrap();
    let mut previous = stepper.total_energy(&state);
    for step in 1..=num_steps {
        state = stepper.step(&state, &DVector::zeros(1)).unwrap();
        let energy = stepper.total_energy(&state);
        assert!(
            energy <= previous * (1.0 + 1e-12),
            "energy rose from {previous:.6e} to {energy:.6e} at step {step}"
        );
        previous = energy;
    }
}

/// Measured order of accuracy of the damped response.
fn convergence_rates(gamma: f64, beta: f64) -> (Vec<f64>, Vec<f64>) {
    let zeta = 0.05;
    let c = 2.0 * zeta * (K * M).sqrt();
    let wn = (K / M).sqrt();
    let period = 2.0 * std::f64::consts::PI / wn;
    let end_time = 4.0 * period;

    let mut errors = Vec::new();
    for refinement in [25usize, 50, 100, 200, 400] {
        let num_steps = 4 * refinement;
        let dt = end_time / num_steps as f64;
        let stepper = sdof(M, c, K, gamma, beta, dt);
        let history = march_sdof(&stepper, X0, V0, num_steps);
        errors.push(max_error(&history, M, c, K, X0, V0));
    }

    let rates = errors
        .windows(2)
        .map(|w| (w[0] / w[1]).log2())
        .collect::<Vec<_>>();
    (errors, rates)
}

/// Halving the time step must quarter the error: Newmark at `γ = ½` is
/// second order.
///
/// The companion assertion below runs the identical measurement at `γ = 0.6`
/// and confirms it reports first order, which is what makes this test
/// evidence rather than decoration.
#[test]
fn newmark_is_second_order_at_gamma_one_half() {
    let (errors, rates) = convergence_rates(0.5, 0.25);
    println!("gamma = 0.5 errors: {errors:?}");
    println!("gamma = 0.5 rates:  {rates:?}");

    for &rate in &rates {
        assert!(
            (1.85..=2.15).contains(&rate),
            "observed order {rate:.3}, expected 2. errors: {errors:?}"
        );
    }
}

#[test]
fn a_wrong_gamma_collapses_the_order_to_one() {
    let (errors, rates) = convergence_rates(0.6, 0.3025);
    println!("gamma = 0.6 errors: {errors:?}");
    println!("gamma = 0.6 rates:  {rates:?}");

    for &rate in &rates {
        assert!(
            rate < 1.3,
            "gamma = 0.6 reported order {rate:.3}; if a mis-applied gamma still \
             measures second order then the order test above proves nothing. \
             errors: {errors:?}"
        );
    }
}

// ---------------------------------------------------------------------------
// The fixed-free bar, shared with modal_closed_form.rs
// ---------------------------------------------------------------------------

const E: f64 = 200e9;
const RHO: f64 = 8000.0;
const LENGTH: f64 = 1.0;
const HEIGHT: f64 = 0.05;

fn bar_mesh(nx: usize, ny: usize) -> (Mesh, Vec<Vec<NodeId>>) {
    let mut mesh = Mesh::new(2).unwrap();
    let mut grid = vec![vec![NodeId(0); ny + 1]; nx + 1];
    for (i, column) in grid.iter_mut().enumerate() {
        for (j, slot) in column.iter_mut().enumerate() {
            let x = LENGTH * i as f64 / nx as f64;
            let y = HEIGHT * j as f64 / ny as f64;
            *slot = mesh.add_node(Node::new_2d(x, y));
        }
    }
    for i in 0..nx {
        for j in 0..ny {
            let nodes = vec![
                grid[i][j],
                grid[i + 1][j],
                grid[i + 1][j + 1],
                grid[i][j + 1],
            ];
            mesh.add_element(Element::new(ElementType::Quad4, nodes, MaterialId(0)).unwrap())
                .unwrap();
        }
    }
    (mesh, grid)
}

/// Poisson's ratio zero, so the plane-stress model reduces exactly to the 1-D
/// bar the closed form describes. Same choice as `modal_closed_form.rs`.
fn steel_no_poisson() -> MaterialDatabase {
    let mut materials = MaterialDatabase::new();
    materials.add_material(
        MaterialId(0),
        LinearElastic::new(E, 0.0).with_density(RHO),
        Some("steel".to_string()),
    );
    materials
}

/// Clamp the left edge axially and suppress transverse motion everywhere.
fn axial_bar_constraints(grid: &[Vec<NodeId>]) -> BoundaryConditionSet {
    let mut bcs = BoundaryConditionSet::new();
    bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
        grid[0].clone(),
        vec![DofComponent::DisplacementX],
        0.0,
    )));
    bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
        grid.iter().flatten().copied().collect(),
        vec![DofComponent::DisplacementY],
        0.0,
    )));
    bcs
}

/// The first axial mode shape `sin(πx/2L)`, as nodal initial displacements.
fn first_mode_initial_displacement(
    grid: &[Vec<NodeId>],
    amplitude: f64,
) -> Vec<(NodeId, DofComponent, f64)> {
    let nx = grid.len() - 1;
    let mut entries = Vec::new();
    for (i, column) in grid.iter().enumerate() {
        let x = LENGTH * i as f64 / nx as f64;
        let value = amplitude * (std::f64::consts::PI * x / (2.0 * LENGTH)).sin();
        for &node in column {
            entries.push((node, DofComponent::DisplacementX, value));
        }
    }
    entries
}

fn series(results: &rtx_fea::analysis::AnalysisResults, key: &str) -> Vec<f64> {
    match results.additional_data.get(key).expect("missing series") {
        AnalysisData::TimeSeries(v) => v.clone(),
        other => panic!("{key} had unexpected type {other:?}"),
    }
}

fn history_row(results: &rtx_fea::analysis::AnalysisResults, key: &str, dof: usize) -> Vec<f64> {
    match results.additional_data.get(key).expect("missing history") {
        AnalysisData::Matrix(m) => m.row(dof).iter().copied().collect(),
        other => panic!("{key} had unexpected type {other:?}"),
    }
}

/// Period from linearly interpolated zero crossings, averaged over every
/// half cycle in the record.
fn period_from_zero_crossings(times: &[f64], values: &[f64]) -> f64 {
    let mut crossings = Vec::new();
    for i in 0..times.len() - 1 {
        let (a, b) = (values[i], values[i + 1]);
        if a != 0.0 && (a < 0.0) != (b < 0.0) {
            let fraction = a / (a - b);
            crossings.push(times[i] + fraction * (times[i + 1] - times[i]));
        }
    }
    assert!(
        crossings.len() >= 3,
        "only {} zero crossings in the record; the response is not oscillating",
        crossings.len()
    );
    let half_cycles = crossings.len() - 1;
    2.0 * (crossings[half_cycles] - crossings[0]) / half_cycles as f64
}

/// Fundamental frequency of the same bar, from the validated modal analysis.
fn modal_fundamental(nx: usize, ny: usize) -> f64 {
    let (mesh, grid) = bar_mesh(nx, ny);
    let mut analysis = ModalAnalysis::new(mesh, steel_no_poisson(), 1, AnalysisConfig::default())
        .with_boundary_conditions(axial_bar_constraints(&grid));
    let results = analysis.run().expect("modal analysis failed");
    match results.additional_data.get("frequencies").unwrap() {
        AnalysisData::Vector(v) => v[0],
        other => panic!("frequencies had unexpected type {other:?}"),
    }
}

/// **The check that matters.** The free-vibration period measured in the time
/// domain must equal `1/f₁` from the frequency domain.
#[test]
fn bar_free_vibration_period_matches_the_modal_analysis() {
    let (nx, ny) = (24, 2);

    let f1 = modal_fundamental(nx, ny);
    let modal_period = 1.0 / f1;

    // Sanity: the modal answer itself is near the continuum closed form, so a
    // later mismatch cannot be blamed on a broken reference.
    let continuum = (E / RHO).sqrt() / (4.0 * LENGTH);
    assert!(
        (f1 - continuum).abs() / continuum < 0.01,
        "modal reference {f1:.4} Hz is not within 1% of the closed form {continuum:.4} Hz"
    );

    let (mesh, grid) = bar_mesh(nx, ny);
    let dt = modal_period / 200.0;
    let num_steps = 200 * 10; // ten periods

    let time_config = TimeConfig {
        start_time: 0.0,
        end_time: num_steps as f64 * dt,
        time_step: dt,
        integration_scheme: TimeIntegrationScheme::NewmarkAverage,
    };

    let mut analysis = DynamicAnalysis::new(
        mesh,
        steel_no_poisson(),
        axial_bar_constraints(&grid),
        time_config,
        AnalysisConfig::default(),
    )
    .with_initial_displacements(first_mode_initial_displacement(&grid, 1e-4));

    let tip = analysis
        .dof_index(grid[nx][0], DofComponent::DisplacementX)
        .unwrap()
        .expect("tip DOF");

    let results = analysis.run().expect("dynamic analysis failed");
    let times = series(&results, "time");
    let displacement = history_row(&results, "displacement_history", tip);

    let measured_period = period_from_zero_crossings(&times, &displacement);
    let measured_frequency = 1.0 / measured_period;
    let elongation = (measured_period - modal_period) / modal_period;

    // The residual disagreement is not noise: the average-acceleration scheme
    // elongates the period by Ω²/12 with Ω = ω₁Δt, and with 200 steps per
    // period that predicts 8.22e-5. Asserting the *signed* difference against
    // that prediction is far stronger than a loose tolerance — it says the
    // remaining gap is the known discretisation error of this scheme and
    // nothing else.
    let omega_dt = 2.0 * std::f64::consts::PI * f1 * dt;
    let predicted = omega_dt * omega_dt / 12.0;

    println!(
        "modal f1 = {f1:.6} Hz (T = {modal_period:.9} s); \
         time-domain T = {measured_period:.9} s (f = {measured_frequency:.6} Hz); \
         relative difference {:.5}% (signed {elongation:.4e}, \
         predicted Newmark elongation {predicted:.4e})",
        elongation.abs() * 100.0
    );

    assert!(
        elongation.abs() < 1e-3,
        "time-domain period {measured_period:.9} s disagrees with the modal \
         1/f1 = {modal_period:.9} s by {:.4}%",
        elongation.abs() * 100.0
    );
    assert!(
        elongation > 0.0,
        "the average-acceleration scheme can only elongate the period, yet the \
         time-domain answer came in {:.4e} short of 1/f1",
        -elongation
    );
    assert!(
        (elongation / predicted - 1.0).abs() < 0.25,
        "the measured period elongation {elongation:.4e} is not the predicted \
         (omega*dt)^2/12 = {predicted:.4e}; the leftover disagreement between \
         the time and frequency domains is something other than Newmark's own \
         discretisation error"
    );

    // The response must be a real oscillation, not a decaying artefact.
    let peak = displacement.iter().map(|x| x.abs()).fold(0.0, f64::max);
    let late_peak = displacement[displacement.len() - 400..]
        .iter()
        .map(|x| x.abs())
        .fold(0.0, f64::max);
    println!("tip amplitude: peak {peak:.6e}, final two periods {late_peak:.6e}");
    assert!(
        (late_peak - peak).abs() / peak < 1e-3,
        "the undamped amplitude fell from {peak:.6e} to {late_peak:.6e}"
    );

    // Undamped energy must hold over the whole march.
    let energy = series(&results, "total_energy");
    let drift = energy
        .iter()
        .map(|e| (e - energy[0]).abs() / energy[0])
        .fold(0.0, f64::max);
    println!("bar energy drift over ten periods: {drift:.3e} relative");
    assert!(drift < 1e-8, "energy drifted by {drift:.3e} relative");
}

/// Measured fundamental damping ratio of the bar under Rayleigh damping.
///
/// Least squares on `ln E(t)`; for a lightly damped single mode
/// `E ∝ e^{-2ζω₁t}`.
fn measured_zeta(alpha: f64, beta: f64, omega1: f64, period: f64) -> f64 {
    let (nx, ny) = (24, 2);
    let (mesh, grid) = bar_mesh(nx, ny);
    let dt = period / 200.0;
    let num_steps = 200 * 8;

    let time_config = TimeConfig {
        start_time: 0.0,
        end_time: num_steps as f64 * dt,
        time_step: dt,
        integration_scheme: TimeIntegrationScheme::NewmarkAverage,
    };

    let mut analysis = DynamicAnalysis::new(
        mesh,
        steel_no_poisson(),
        axial_bar_constraints(&grid),
        time_config,
        AnalysisConfig::default(),
    )
    .with_rayleigh_damping(alpha, beta)
    .with_initial_displacements(first_mode_initial_displacement(&grid, 1e-4));

    let results = analysis.run().expect("dynamic analysis failed");
    let times = series(&results, "time");
    let energy = series(&results, "total_energy");

    let n = times.len() as f64;
    let mean_t = times.iter().sum::<f64>() / n;
    let logs: Vec<f64> = energy.iter().map(|e| e.ln()).collect();
    let mean_l = logs.iter().sum::<f64>() / n;
    let numerator: f64 = times
        .iter()
        .zip(&logs)
        .map(|(t, l)| (t - mean_t) * (l - mean_l))
        .sum();
    let denominator: f64 = times.iter().map(|t| (t - mean_t).powi(2)).sum();

    -(numerator / denominator) / (2.0 * omega1)
}

/// `C = αM + βK` must be assembled term by term.
///
/// Mass-proportional damping gives `ζ₁ = α/(2ω₁)` and stiffness-proportional
/// damping gives `ζ₁ = βω₁/2`. Both are targeted at the same `ζ₁ = 0.02`, so
/// swapping the coefficients — which a single combined test would not catch —
/// misses by a factor of `ω₁² ≈ 6·10⁷`.
#[test]
fn rayleigh_damping_reproduces_the_modal_damping_ratios() {
    let f1 = modal_fundamental(24, 2);
    let omega1 = 2.0 * std::f64::consts::PI * f1;
    let period = 1.0 / f1;
    let target = 0.02;

    let alpha = 2.0 * target * omega1;
    let mass_zeta = measured_zeta(alpha, 0.0, omega1, period);
    println!(
        "mass-proportional alpha = {alpha:.4}: target zeta {target}, \
         measured {mass_zeta:.5} ({:.2}% off)",
        (mass_zeta - target).abs() / target * 100.0
    );
    assert!(
        (mass_zeta - target).abs() / target < 0.05,
        "mass-proportional damping gave zeta = {mass_zeta:.5}, expected {target}"
    );

    let beta = 2.0 * target / omega1;
    let stiffness_zeta = measured_zeta(0.0, beta, omega1, period);
    println!(
        "stiffness-proportional beta = {beta:.6e}: target zeta {target}, \
         measured {stiffness_zeta:.5} ({:.2}% off)",
        (stiffness_zeta - target).abs() / target * 100.0
    );
    assert!(
        (stiffness_zeta - target).abs() / target < 0.05,
        "stiffness-proportional damping gave zeta = {stiffness_zeta:.5}, expected {target}"
    );
}

/// A free-free bar given a uniform initial velocity must translate rigidly.
///
/// Rigid translation stores no strain energy. This is the assertion a zero or
/// transposed assembly fails while still looking symmetric and finite: any
/// spurious coupling turns the translation into vibration and puts energy
/// into `½uᵀKu`.
#[test]
fn rigid_translation_of_a_free_bar_stores_no_strain_energy() {
    let (nx, ny) = (8, 2);
    let (mesh, grid) = bar_mesh(nx, ny);
    let speed = 3.0;

    let velocities: Vec<_> = grid
        .iter()
        .flatten()
        .map(|&node| (node, DofComponent::DisplacementX, speed))
        .collect();

    let dt = 1e-5;
    let num_steps = 100;
    let time_config = TimeConfig {
        start_time: 0.0,
        end_time: num_steps as f64 * dt,
        time_step: dt,
        integration_scheme: TimeIntegrationScheme::NewmarkAverage,
    };

    let mut analysis = DynamicAnalysis::new(
        mesh,
        steel_no_poisson(),
        BoundaryConditionSet::new(),
        time_config,
        AnalysisConfig::default(),
    )
    .with_initial_velocities(velocities);

    let tip = analysis
        .dof_index(grid[nx][0], DofComponent::DisplacementX)
        .unwrap()
        .expect("tip DOF");

    let results = analysis.run().expect("dynamic analysis failed");
    let times = series(&results, "time");
    let displacement = history_row(&results, "displacement_history", tip);
    let strain = series(&results, "strain_energy");
    let kinetic = series(&results, "kinetic_energy");

    let worst_position = times
        .iter()
        .zip(&displacement)
        .map(|(t, u)| (u - speed * t).abs())
        .fold(0.0, f64::max);

    // The comparator is the strain energy a *stretch* reaching the same tip
    // displacement would store, `½(EA/L)u²` with unit thickness. Comparing
    // against the kinetic energy instead is the wrong yardstick: the residual
    // here is the floating-point cancellation floor of the quadratic form
    // `½uᵀKu` itself, which scales with `‖K‖‖u‖²` — that is, with the stretch
    // energy — and not with `½vᵀMv`. Measured against `‖K‖‖u‖²` the residual
    // is round-off; measured against the kinetic energy it looks like `1e-12`
    // for no physical reason.
    let axial_stiffness = E * HEIGHT / LENGTH;
    let worst_ratio = times
        .iter()
        .zip(&strain)
        .skip(1)
        .map(|(t, s)| s.abs() / (0.5 * axial_stiffness * (speed * t).powi(2)))
        .fold(0.0, f64::max);

    println!(
        "rigid translation: worst position error {worst_position:.3e} m, \
         peak strain energy {:.3e} J against the equivalent stretch energy \
         {:.3e} J (ratio {worst_ratio:.3e}), kinetic energy {:.6e} J",
        strain.iter().fold(0.0f64, |a, b| a.max(b.abs())),
        0.5 * axial_stiffness * (speed * times[times.len() - 1]).powi(2),
        kinetic[0]
    );

    // Analytic kinetic energy: half the total mass times v squared. The
    // plane-stress model has unit thickness.
    let exact_kinetic = 0.5 * RHO * LENGTH * HEIGHT * speed * speed;
    assert!(
        (kinetic[0] - exact_kinetic).abs() / exact_kinetic < 1e-10,
        "kinetic energy {:.6e} J is not the analytic {exact_kinetic:.6e} J; \
         the consistent mass matrix does not integrate to the bar's mass",
        kinetic[0]
    );

    assert!(
        worst_position < 1e-12 * speed * times[times.len() - 1],
        "rigid translation drifted by {worst_position:.3e} m"
    );
    // Round-off in `½uᵀKu` sits near `ε‖K‖‖u‖²/(½(EA/L)u²) ≈ 1e-13`. Any real
    // spurious coupling — a transposed assembly, a mis-scattered mass — turns
    // the translation into vibration and lands orders of magnitude above this.
    assert!(
        worst_ratio < 1e-10,
        "rigid translation stored strain energy at {worst_ratio:.3e} of the \
         energy an equivalent stretch would store; that is far above the \
         round-off floor of the quadratic form"
    );
}

/// A configuration that produces no steps is an error, not an empty result.
#[test]
fn a_degenerate_time_configuration_is_rejected() {
    let (mesh, grid) = bar_mesh(4, 1);
    let time_config = TimeConfig {
        start_time: 0.0,
        end_time: 0.0,
        time_step: 1e-5,
        integration_scheme: TimeIntegrationScheme::NewmarkAverage,
    };
    let mut analysis = DynamicAnalysis::new(
        mesh,
        steel_no_poisson(),
        axial_bar_constraints(&grid),
        time_config,
        AnalysisConfig::default(),
    );
    let error = analysis
        .run()
        .expect_err("a zero-length time window must not yield a response");
    assert!(
        error.to_string().to_lowercase().contains("no steps"),
        "error should name the empty time window, got: {error}"
    );
}