//! Code verification by the Method of Manufactured Solutions.
//!
//! Every other test in this crate asserts that some quantity is *close enough*
//! to a value someone believed. MMS asserts something stronger and more
//! objective: that the discretisation converges to the exact solution at the
//! rate the theory predicts.
//!
//! The method needs no benchmark data. Choose any sufficiently smooth field
//! `u`, substitute it into the governing equations, and whatever they fail to
//! balance is the body force `f = -∇·σ(u)` that makes `u` the exact solution of
//! the problem with that load. Solve with `f` applied and `u` imposed on the
//! boundary, measure the error, refine, and read off the observed order
//! `log(e_h / e_{h/2}) / log(2)`.
//!
//! This matters here specifically. Every defect found in this crate returned a
//! plausible wrong answer and passed the tests written for it — element
//! matrices of zeros, a quadrature rule with no points, an assembly transposed
//! per node. None of them survives an order-of-accuracy measurement: a wrong
//! discretisation either fails to converge or converges at the wrong rate, and
//! neither can be tuned away.
//!
//! # What is measured
//!
//! Four element types, each against **its own** theoretical L2 displacement
//! rate — `h^(p+1)` for a space complete to polynomial degree `p`:
//!
//! | element | space                       | expected L2 order |
//! |---------|-----------------------------|-------------------|
//! | `Tri3`  | linear triangle, 2-D        | 2                 |
//! | `Quad4` | bilinear quadrilateral, 2-D | 2                 |
//! | `Quad8` | serendipity quad, 2-D       | 3                 |
//! | `Hex8`  | trilinear hexahedron, 3-D   | 2                 |
//!
//! A blanket "order ≥ 1.8" across all four would pass for `Quad8` while it was
//! silently no better than `Quad4`, so each element asserts its own rate.
//!
//! # The manufactured solutions
//!
//! Two dimensions (plane stress), shared by `Tri3`, `Quad4` and `Quad8`:
//!
//! ```text
//!     u(x, y) = sin(pi x) sin(pi y)
//!     v(x, y) = x^2 (1 - x) y (1 - y)
//! ```
//!
//! Three dimensions (`Hex8`):
//!
//! ```text
//!     u(x, y, z) = sin(pi x) sin(pi y) sin(pi z)
//!     v(x, y, z) = cos(pi x) sin(2 pi y) z^2
//!     w(x, y, z) = x^2 y sin(pi z)
//! ```
//!
//! Both are chosen so the field is smooth but **not** in any of the element
//! spaces, so the error measured is genuinely discretisation error rather than
//! round-off. The components are deliberately different in form and every
//! mixed derivative `du_i/dx_j` is non-zero, so all three shear rows and the
//! full off-diagonal `lambda` block of the 3-D constitutive matrix are
//! exercised rather than silently skipped.
//!
//! Neither field vanishes on the boundary. That is deliberate: the exact
//! solution is imposed as a **non-zero** Dirichlet condition, which is the
//! case a solver that drops the `K_fc u_c` coupling term gets wrong.

use nalgebra::{DVector, Vector3};
use rtx_fea::assembly::{
    AdvancedDofNumbering, AssemblyOptions, DofComponent, DofMappingStrategy, GlobalAssembler,
};
use rtx_fea::elements::{ElementMatrixComputer, FiniteElement, StandardFiniteElement};
use rtx_fea::materials::{LinearElastic, MaterialDatabase};
use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
use rtx_fea::solvers::{SolverFactory, SolverOptions};
use std::f64::consts::PI;

const E: f64 = 1.0;
const NU: f64 = 0.3;

/// Lamé's first parameter, as `elasticity_matrix` builds it for `spatial_dim = 3`.
const LAMBDA: f64 = E * NU / ((1.0 + NU) * (1.0 - 2.0 * NU));
/// Shear modulus. Also the plane-stress shear entry `c (1 - nu) / 2`.
const MU: f64 = E / (2.0 * (1.0 + NU));

// ---------------------------------------------------------------------------
// Two-dimensional manufactured solution (plane stress)
// ---------------------------------------------------------------------------

/// Exact displacement field.
fn exact(x: f64, y: f64) -> (f64, f64) {
    (
        (PI * x).sin() * (PI * y).sin(),
        x * x * (1.0 - x) * y * (1.0 - y),
    )
}

/// Body force `f = -div(sigma(u))` for plane stress.
///
/// Derived by hand from the manufactured field above with
/// `sigma_xx = c(e_xx + nu e_yy)`, `sigma_yy = c(nu e_xx + e_yy)`,
/// `sigma_xy = mu gamma_xy`, where `c = E / (1 - nu^2)` and
/// `mu = E / (2(1 + nu))` — the same plane-stress constitutive matrix
/// `ElementMatrixComputer::elasticity_matrix` builds.
///
/// Strains:
/// ```text
///     e_xx = pi cos(pi x) sin(pi y)
///     e_yy = x^2 (1 - x) (1 - 2y)
///     gamma_xy = pi sin(pi x) cos(pi y) + (2x - 3x^2) y (1 - y)
/// ```
fn body_force(x: f64, y: f64) -> (f64, f64) {
    let c = E / (1.0 - NU * NU);
    let mu = MU;

    let (sx, cx) = ((PI * x).sin(), (PI * x).cos());
    let (sy, cy) = ((PI * y).sin(), (PI * y).cos());

    // d(e_xx)/dx, d(e_yy)/dx, d(e_xx)/dy, d(e_yy)/dy
    let dexx_dx = -PI * PI * sx * sy;
    let deyy_dx = (2.0 * x - 3.0 * x * x) * (1.0 - 2.0 * y);
    let dexx_dy = PI * PI * cx * cy;
    let deyy_dy = -2.0 * x * x * (1.0 - x);

    // d(gamma_xy)/dy and d(gamma_xy)/dx
    let dgxy_dy = -PI * PI * sx * sy + (2.0 * x - 3.0 * x * x) * (1.0 - 2.0 * y);
    let dgxy_dx = PI * PI * cx * cy + (2.0 - 6.0 * x) * y * (1.0 - y);

    let dsxx_dx = c * (dexx_dx + NU * deyy_dx);
    let dsyy_dy = c * (NU * dexx_dy + deyy_dy);
    let dsxy_dy = mu * dgxy_dy;
    let dsxy_dx = mu * dgxy_dx;

    (-(dsxx_dx + dsxy_dy), -(dsxy_dx + dsyy_dy))
}

/// The 2-D exact field as a spatial function.
fn exact_2d(p: Vector3<f64>) -> Vector3<f64> {
    let (u, v) = exact(p.x, p.y);
    Vector3::new(u, v, 0.0)
}

/// The 2-D body force as a spatial function.
fn force_2d(p: Vector3<f64>) -> Vector3<f64> {
    let (fx, fy) = body_force(p.x, p.y);
    Vector3::new(fx, fy, 0.0)
}

// ---------------------------------------------------------------------------
// Three-dimensional manufactured solution (isotropic, lambda/mu form)
// ---------------------------------------------------------------------------

/// Exact 3-D displacement field.
///
/// ```text
///     u1 = sin(pi x) sin(pi y) sin(pi z)
///     u2 = cos(pi x) sin(2 pi y) z^2
///     u3 = x^2 y sin(pi z)
/// ```
fn exact_3d(p: Vector3<f64>) -> Vector3<f64> {
    let (x, y, z) = (p.x, p.y, p.z);
    Vector3::new(
        (PI * x).sin() * (PI * y).sin() * (PI * z).sin(),
        (PI * x).cos() * (2.0 * PI * y).sin() * z * z,
        x * x * y * (PI * z).sin(),
    )
}

/// Body force `f = -div(sigma(u))` for 3-D isotropic elasticity.
///
/// For `sigma = lambda tr(eps) I + 2 mu eps` — which is exactly the 6x6 matrix
/// `elasticity_matrix` builds for `spatial_dim = 3`, with `lambda + 2 mu` on
/// the first three diagonal entries, `lambda` off-diagonal and `mu` on the
/// three shear entries — the divergence collapses to the Navier form
///
/// ```text
///     div sigma = (lambda + mu) grad(div u) + mu laplacian(u)
/// ```
///
/// so `f = -[(lambda + mu) grad(div u) + mu laplacian(u)]`. Both terms are
/// written out below directly from the field above; nothing here is
/// differentiated numerically, and `body_force_3d_matches_stress_divergence`
/// checks the whole expression against an independent numerical divergence of
/// `sigma`.
fn body_force_3d(p: Vector3<f64>) -> Vector3<f64> {
    let (x, y, z) = (p.x, p.y, p.z);

    let (sx, cx) = ((PI * x).sin(), (PI * x).cos());
    let (sy, cy) = ((PI * y).sin(), (PI * y).cos());
    let (sz, cz) = ((PI * z).sin(), (PI * z).cos());
    let (s2y, c2y) = ((2.0 * PI * y).sin(), (2.0 * PI * y).cos());

    // Laplacians.
    //   u1 is a product of three sines of argument pi * (coordinate), so each
    //   second derivative contributes -pi^2 u1.
    let lap_u1 = -3.0 * PI * PI * sx * sy * sz;
    //   u2 = cos(pi x) sin(2 pi y) z^2:  -pi^2 u2 - 4 pi^2 u2 + 2 cos sin
    let lap_u2 = cx * s2y * (2.0 - 5.0 * PI * PI * z * z);
    //   u3 = x^2 y sin(pi z):  2 y sin(pi z) + 0 - pi^2 x^2 y sin(pi z)
    let lap_u3 = y * sz * (2.0 - PI * PI * x * x);

    // grad(div u), where
    //   div u = pi cos(pi x) sin(pi y) sin(pi z)
    //         + 2 pi cos(pi x) cos(2 pi y) z^2
    //         + pi x^2 y cos(pi z)
    let ddiv_dx =
        -PI * PI * sx * sy * sz - 2.0 * PI * PI * sx * c2y * z * z + 2.0 * PI * x * y * cz;
    let ddiv_dy = PI * PI * cx * cy * sz - 4.0 * PI * PI * cx * s2y * z * z + PI * x * x * cz;
    let ddiv_dz = PI * PI * cx * sy * cz + 4.0 * PI * cx * c2y * z - PI * PI * x * x * y * sz;

    Vector3::new(
        -((LAMBDA + MU) * ddiv_dx + MU * lap_u1),
        -((LAMBDA + MU) * ddiv_dy + MU * lap_u2),
        -((LAMBDA + MU) * ddiv_dz + MU * lap_u3),
    )
}

// ---------------------------------------------------------------------------
// Meshes
// ---------------------------------------------------------------------------

/// Unit square meshed with `n` by `n` `Quad4` elements.
fn quad4_mesh(n: usize) -> Mesh {
    let mut mesh = Mesh::new(2).unwrap();
    let mut grid = vec![vec![NodeId(0); n + 1]; n + 1];
    for (i, column) in grid.iter_mut().enumerate() {
        for (j, slot) in column.iter_mut().enumerate() {
            let x = i as f64 / n as f64;
            let y = j as f64 / n as f64;
            *slot = mesh.add_node(Node::new_2d(x, y));
        }
    }
    for i in 0..n {
        for j in 0..n {
            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
}

/// Unit square meshed with `2 n^2` `Tri3` elements.
///
/// Each cell of the same `n` by `n` grid is split along its diagonal, so the
/// node set — and therefore the DOF count — is identical to the `Quad4` mesh
/// at the same `n`. Both triangles are listed counter-clockwise, which keeps
/// the Jacobian determinant positive; a clockwise listing gives a negative
/// determinant that `det.abs()` would hide while the strain-displacement
/// matrix silently changed sign.
fn tri3_mesh(n: usize) -> Mesh {
    let mut mesh = Mesh::new(2).unwrap();
    let mut grid = vec![vec![NodeId(0); n + 1]; n + 1];
    for (i, column) in grid.iter_mut().enumerate() {
        for (j, slot) in column.iter_mut().enumerate() {
            let x = i as f64 / n as f64;
            let y = j as f64 / n as f64;
            *slot = mesh.add_node(Node::new_2d(x, y));
        }
    }
    for i in 0..n {
        for j in 0..n {
            for nodes in [
                vec![grid[i][j], grid[i + 1][j], grid[i + 1][j + 1]],
                vec![grid[i][j], grid[i + 1][j + 1], grid[i][j + 1]],
            ] {
                mesh.add_element(Element::new(ElementType::Tri3, nodes, MaterialId(0)).unwrap())
                    .unwrap();
            }
        }
    }
    mesh
}

/// Unit square meshed with `n` by `n` `Quad8` serendipity elements.
///
/// Nodes sit on a `(2n+1)` by `(2n+1)` lattice with the cell centres left out,
/// which is what makes the element serendipity rather than `Quad9`. Ordering
/// per element is the standard one the shape functions assume: four corners
/// counter-clockwise from `(-1, -1)`, then the four mid-edge nodes starting
/// with the edge between corners 0 and 1.
fn quad8_mesh(n: usize) -> Mesh {
    let mut mesh = Mesh::new(2).unwrap();
    let lattice = 2 * n + 1;
    let mut grid = vec![vec![None; lattice]; lattice];
    for (i, column) in grid.iter_mut().enumerate() {
        for (j, slot) in column.iter_mut().enumerate() {
            if i % 2 == 1 && j % 2 == 1 {
                continue; // cell centre: no node in a serendipity element
            }
            let x = i as f64 / (2 * n) as f64;
            let y = j as f64 / (2 * n) as f64;
            *slot = Some(mesh.add_node(Node::new_2d(x, y)));
        }
    }
    for i in 0..n {
        for j in 0..n {
            let (a, b) = (2 * i, 2 * j);
            let nodes = vec![
                grid[a][b].unwrap(),
                grid[a + 2][b].unwrap(),
                grid[a + 2][b + 2].unwrap(),
                grid[a][b + 2].unwrap(),
                grid[a + 1][b].unwrap(),
                grid[a + 2][b + 1].unwrap(),
                grid[a + 1][b + 2].unwrap(),
                grid[a][b + 1].unwrap(),
            ];
            mesh.add_element(Element::new(ElementType::Quad8, nodes, MaterialId(0)).unwrap())
                .unwrap();
        }
    }
    mesh
}

/// Unit cube meshed with `n` by `n` by `n` `Hex8` elements.
///
/// Node ordering per element is the standard trilinear one: the four nodes of
/// the `zeta = -1` face counter-clockwise from `(-1, -1, -1)`, then the four
/// of the `zeta = +1` face in the same order.
fn hex8_mesh(n: usize) -> Mesh {
    let mut mesh = Mesh::new(3).unwrap();
    let mut grid = vec![vec![vec![NodeId(0); n + 1]; n + 1]; n + 1];
    for (i, plane) in grid.iter_mut().enumerate() {
        for (j, column) in plane.iter_mut().enumerate() {
            for (k, slot) in column.iter_mut().enumerate() {
                let x = i as f64 / n as f64;
                let y = j as f64 / n as f64;
                let z = k as f64 / n as f64;
                *slot = mesh.add_node(Node::new_3d(x, y, z));
            }
        }
    }
    for i in 0..n {
        for j in 0..n {
            for k in 0..n {
                let nodes = vec![
                    grid[i][j][k],
                    grid[i + 1][j][k],
                    grid[i + 1][j + 1][k],
                    grid[i][j + 1][k],
                    grid[i][j][k + 1],
                    grid[i + 1][j][k + 1],
                    grid[i + 1][j + 1][k + 1],
                    grid[i][j + 1][k + 1],
                ];
                mesh.add_element(Element::new(ElementType::Hex8, nodes, MaterialId(0)).unwrap())
                    .unwrap();
            }
        }
    }
    mesh
}

// ---------------------------------------------------------------------------
// Solve and measure
// ---------------------------------------------------------------------------

/// Whether a node of the unit square/cube lies on its boundary.
fn on_boundary(p: Vector3<f64>, dim: usize) -> bool {
    (0..dim).any(|d| p[d].abs() < 1e-12 || (p[d] - 1.0).abs() < 1e-12)
}

/// Solve the manufactured problem on `mesh` and return the L2 displacement
/// error, integrated with the elements' own quadrature.
///
/// `load_quadrature` and `error_quadrature` are passed straight through to
/// `quadrature_rule`. They are explicit because the right rule is a property
/// of the element: the consistent load vector must be integrated at least as
/// accurately as the stiffness or the load error, not the element, sets the
/// observed order.
fn solve_mms(
    mesh: &Mesh,
    exact_field: &dyn Fn(Vector3<f64>) -> Vector3<f64>,
    force_field: &dyn Fn(Vector3<f64>) -> Vector3<f64>,
    load_quadrature: Option<usize>,
    error_quadrature: Option<usize>,
) -> f64 {
    let dim = mesh.spatial_dimension;

    let mut materials = MaterialDatabase::new();
    materials.add_material(MaterialId(0), LinearElastic::new(E, NU), None);

    let mut dof_numbering =
        AdvancedDofNumbering::displacement_only(mesh, DofMappingStrategy::Sequential).unwrap();

    let components = [
        DofComponent::DisplacementX,
        DofComponent::DisplacementY,
        DofComponent::DisplacementZ,
    ];

    // Impose the exact solution on the whole boundary.
    let mut prescribed = Vec::new();
    for (&node_id, node) in &mesh.nodes {
        let position = node.position();
        if !on_boundary(position, dim) {
            continue;
        }
        let value = exact_field(position);
        for (component_index, &component) in components[..dim].iter().enumerate() {
            let dof = dof_numbering.get_dof(node_id, component).unwrap();
            dof_numbering.constrain_dof(dof).unwrap();
            prescribed.push((dof, value[component_index]));
        }
    }

    let simple = dof_numbering.to_dof_numbering();
    let assembler = GlobalAssembler::new(AssemblyOptions::default(), materials);
    let mut system = assembler.assemble_system(mesh, &simple, 0.0).unwrap();

    // Consistent body force: integral of N_i f over each element.
    for element in mesh.elements.values() {
        let node_coords: Vec<Vector3<f64>> = element
            .nodes
            .iter()
            .map(|id| mesh.get_node(*id).unwrap().position())
            .collect();
        let fe = StandardFiniteElement::new(element.element_type, node_coords.clone());
        let local = ElementMatrixComputer::compute_body_force_vector(
            &fe,
            &node_coords,
            force_field,
            load_quadrature,
        )
        .unwrap();

        let mut global_dofs = Vec::new();
        for node in &element.nodes {
            global_dofs.extend_from_slice(simple.get_node_dofs(*node).unwrap());
        }
        for (local_index, &global_dof) in global_dofs.iter().enumerate() {
            system.force_vector[global_dof] += local[local_index];
        }
    }

    for (dof, value) in prescribed {
        system.set_prescribed_value(dof, value).unwrap();
    }

    let (free_stiffness, free_force) = system.extract_free_system().unwrap();
    let properties = SolverFactory::analyze_matrix(&free_stiffness);
    let options = SolverOptions::default();
    let mut solver = SolverFactory::create_linear_solver(&properties, &options);
    let (free_solution, _) = solver
        .solve(&free_stiffness, &free_force, &options)
        .unwrap();

    // Expand to the full DOF vector, with the prescribed values in place.
    let mut solution = DVector::zeros(system.dof_numbering.total_dofs);
    for (i, &dof) in system.dof_numbering.free_dofs.iter().enumerate() {
        solution[dof] = free_solution[i];
    }
    for (&dof, &value) in system.prescribed_values() {
        solution[dof] = value;
    }

    // L2 error by the elements' own quadrature.
    let mut squared = 0.0;
    for element in mesh.elements.values() {
        let node_coords: Vec<Vector3<f64>> = element
            .nodes
            .iter()
            .map(|id| mesh.get_node(*id).unwrap().position())
            .collect();
        let fe = StandardFiniteElement::new(element.element_type, node_coords.clone());
        let rule = fe.quadrature_rule(error_quadrature).unwrap();

        let dofs: Vec<usize> = element
            .nodes
            .iter()
            .flat_map(|node| simple.get_node_dofs(*node).unwrap().to_vec())
            .collect();

        for point in &rule.points {
            let shape = fe.shape_functions(&point.coords).unwrap();
            let jacobian = fe.jacobian(&point.coords, &node_coords).unwrap();
            let physical = fe.map_to_physical(&point.coords, &node_coords).unwrap();

            let mut uh: Vector3<f64> = Vector3::zeros();
            for node_index in 0..element.nodes.len() {
                let n = shape.value(node_index).unwrap();
                for d in 0..dim {
                    uh[d] += n * solution[dofs[node_index * dim + d]];
                }
            }

            let exact_value = exact_field(physical.coords);
            let weight = point.weight * jacobian.determinant().abs();
            let mut pointwise = 0.0;
            for d in 0..dim {
                pointwise += (uh[d] - exact_value[d]).powi(2);
            }
            squared += pointwise * weight;
        }
    }

    squared.sqrt()
}

/// L2 error of the **nodal interpolant** of the exact field on `mesh`.
///
/// This is the best the element space can do at this resolution, up to a
/// constant: no choice of nodal values reproduces the exact field more closely
/// than a few times this. It touches no DOF numbering, no assembly and no
/// solver, so comparing the computed error against it isolates "how well can
/// this space represent the answer" from "did the machinery find it".
fn interpolation_error(
    mesh: &Mesh,
    exact_field: &dyn Fn(Vector3<f64>) -> Vector3<f64>,
    error_quadrature: Option<usize>,
) -> f64 {
    let dim = mesh.spatial_dimension;
    let mut squared = 0.0;

    for element in mesh.elements.values() {
        let node_coords: Vec<Vector3<f64>> = element
            .nodes
            .iter()
            .map(|id| mesh.get_node(*id).unwrap().position())
            .collect();
        let fe = StandardFiniteElement::new(element.element_type, node_coords.clone());
        let rule = fe.quadrature_rule(error_quadrature).unwrap();

        for point in &rule.points {
            let shape = fe.shape_functions(&point.coords).unwrap();
            let jacobian = fe.jacobian(&point.coords, &node_coords).unwrap();
            let physical = fe.map_to_physical(&point.coords, &node_coords).unwrap();

            let mut interpolated: Vector3<f64> = Vector3::zeros();
            for (node_index, coord) in node_coords.iter().enumerate() {
                let n = shape.value(node_index).unwrap();
                let nodal = exact_field(*coord);
                for d in 0..dim {
                    interpolated[d] += n * nodal[d];
                }
            }

            let exact_value = exact_field(physical.coords);
            let weight = point.weight * jacobian.determinant().abs();
            let mut pointwise = 0.0;
            for d in 0..dim {
                pointwise += (interpolated[d] - exact_value[d]).powi(2);
            }
            squared += pointwise * weight;
        }
    }

    squared.sqrt()
}

/// Observed order of accuracy between consecutive refinement levels.
///
/// Computed from the mesh-size ratio rather than assuming successive halving,
/// so a non-doubling ladder still reads correctly.
fn observed_orders(resolutions: &[usize], errors: &[f64]) -> Vec<f64> {
    errors
        .windows(2)
        .zip(resolutions.windows(2))
        .map(|(e, n)| (e[0] / e[1]).ln() / (n[1] as f64 / n[0] as f64).ln())
        .collect()
}

/// Print the whole ladder so the measurement is auditable rather than merely
/// asserted.
fn report(label: &str, resolutions: &[usize], errors: &[f64], rates: &[f64]) {
    println!("  {label}");
    for (i, &n) in resolutions.iter().enumerate() {
        let rate = if i == 0 {
            String::from("    -")
        } else {
            format!("{:5.2}", rates[i - 1])
        };
        println!(
            "    n = {n:3}   L2 error = {:.6e}   observed order = {rate}",
            errors[i]
        );
    }
}

/// Run a refinement ladder and assert every observed order lies in `window`.
fn assert_order(
    label: &str,
    resolutions: &[usize],
    errors: &[f64],
    expected: f64,
    window: (f64, f64),
) {
    let rates = observed_orders(resolutions, errors);
    report(label, resolutions, errors, &rates);

    for (i, &rate) in rates.iter().enumerate() {
        assert!(
            rate > window.0 && rate < window.1,
            "{label}: refinement {} -> {}: observed order {rate:.3}, expected {expected} \
             (accepted {:.2}..{:.2}). Errors: {errors:?}",
            resolutions[i],
            resolutions[i + 1],
            window.0,
            window.1,
        );
    }
}

// ---------------------------------------------------------------------------
// Fixture verification: the body force really is -div(sigma(u))
// ---------------------------------------------------------------------------
//
// A wrong body force is the failure mode that looks most like a wrong solver:
// the solution converges, just not to the field being compared against, and
// the observed order drops. These two tests pin the fixture independently of
// every line of library code, by differentiating the *stress* numerically and
// comparing with the analytic expressions above. Fourth-order central
// differences at h = 1e-2 leave a truncation error near 1e-5 on quantities of
// order 10-100 and a round-off error near 1e-12, so a genuine mismatch — a
// dropped term, a sign, a wrong Lamé constant — is orders of magnitude larger
// than the noise floor.

/// Fourth-order central difference of `f` along axis `dir`.
fn d_dx(f: &dyn Fn(Vector3<f64>) -> f64, p: Vector3<f64>, dir: usize, h: f64) -> f64 {
    let at = |steps: f64| {
        let mut q = p;
        q[dir] += steps * h;
        f(q)
    };
    (at(-2.0) - 8.0 * at(-1.0) + 8.0 * at(1.0) - at(2.0)) / (12.0 * h)
}

/// `-div(sigma)` computed by numerical differentiation of a stress function.
fn numerical_force(
    stress: &dyn Fn(Vector3<f64>, usize, usize) -> f64,
    p: Vector3<f64>,
    dim: usize,
    h: f64,
) -> Vector3<f64> {
    let mut f = Vector3::zeros();
    for i in 0..dim {
        let mut divergence = 0.0;
        for j in 0..dim {
            divergence += d_dx(&|q| stress(q, i, j), p, j, h);
        }
        f[i] = -divergence;
    }
    f
}

#[test]
fn body_force_2d_matches_stress_divergence() {
    let c = E / (1.0 - NU * NU);
    let strain = |q: Vector3<f64>, i: usize, j: usize| {
        let ui = |r: Vector3<f64>| exact_2d(r)[i];
        let uj = |r: Vector3<f64>| exact_2d(r)[j];
        0.5 * (d_dx(&ui, q, j, 1e-2) + d_dx(&uj, q, i, 1e-2))
    };
    // Plane stress, exactly as `elasticity_matrix` builds it for dim 2.
    let stress = |q: Vector3<f64>, i: usize, j: usize| {
        let (exx, eyy) = (strain(q, 0, 0), strain(q, 1, 1));
        match (i, j) {
            (0, 0) => c * (exx + NU * eyy),
            (1, 1) => c * (NU * exx + eyy),
            _ => 2.0 * MU * strain(q, i, j),
        }
    };

    let mut worst: f64 = 0.0;
    for &x in &[0.13, 0.37, 0.5, 0.71, 0.94] {
        for &y in &[0.09, 0.28, 0.5, 0.66, 0.87] {
            let p = Vector3::new(x, y, 0.0);
            let numeric = numerical_force(&stress, p, 2, 1e-2);
            let analytic = force_2d(p);
            worst = worst.max((numeric - analytic).norm());
        }
    }

    assert!(
        worst < 1e-4,
        "analytic plane-stress body force disagrees with the numerical \
         divergence of sigma by {worst:.3e}; the manufactured fixture is wrong"
    );
}

#[test]
fn body_force_3d_matches_stress_divergence() {
    let strain = |q: Vector3<f64>, i: usize, j: usize| {
        let ui = |r: Vector3<f64>| exact_3d(r)[i];
        let uj = |r: Vector3<f64>| exact_3d(r)[j];
        0.5 * (d_dx(&ui, q, j, 1e-2) + d_dx(&uj, q, i, 1e-2))
    };
    // sigma = lambda tr(eps) I + 2 mu eps, the lambda/mu form of the 6x6
    // matrix `elasticity_matrix` builds for spatial_dim = 3.
    let stress = |q: Vector3<f64>, i: usize, j: usize| {
        let trace: f64 = (0..3).map(|k| strain(q, k, k)).sum();
        let volumetric = if i == j { LAMBDA * trace } else { 0.0 };
        volumetric + 2.0 * MU * strain(q, i, j)
    };

    let mut worst: f64 = 0.0;
    for &x in &[0.17, 0.5, 0.83] {
        for &y in &[0.11, 0.5, 0.79] {
            for &z in &[0.23, 0.5, 0.91] {
                let p = Vector3::new(x, y, z);
                let numeric = numerical_force(&stress, p, 3, 1e-2);
                let analytic = body_force_3d(p);
                worst = worst.max((numeric - analytic).norm());
            }
        }
    }

    assert!(
        worst < 1e-3,
        "analytic 3-D body force disagrees with the numerical divergence of \
         sigma by {worst:.3e}; the manufactured fixture is wrong"
    );
}

/// The triangle rule actually used below integrates the reference triangle.
///
/// `QuadratureRule::triangle(2)` is the rule the `Tri3` ladder integrates its
/// error with, and its weights must sum to the reference area 1/2. This is
/// checked because the *order-3* triangle rule in the same function does not:
/// its weights sum to 1/4, half of what they should be. That defect is
/// reported rather than fixed here — this file owns no library code — and the
/// `Tri3` ladder avoids order 3 for that reason.
#[test]
fn triangle_error_quadrature_integrates_the_reference_area() {
    let element = StandardFiniteElement::new(
        ElementType::Tri3,
        vec![
            Vector3::new(0.0, 0.0, 0.0),
            Vector3::new(1.0, 0.0, 0.0),
            Vector3::new(0.0, 1.0, 0.0),
        ],
    );
    let rule = element.quadrature_rule(Some(2)).unwrap();
    let area: f64 = rule.points.iter().map(|p| p.weight).sum();
    assert!(
        (area - 0.5).abs() < 1e-12,
        "the 3-point triangle rule integrates the reference triangle to {area}, \
         not 1/2, so every integral taken with it is scaled wrongly"
    );
}

// ---------------------------------------------------------------------------
// Order-of-accuracy ladders
// ---------------------------------------------------------------------------

/// `Quad4`: bilinear, theoretical L2 displacement rate 2.
///
/// A rate near 1 means something in the chain is only first-order accurate; a
/// rate near 0 means the solution is not converging to the exact one at all,
/// which is what every stub found in this crate would produce.
#[test]
fn quad4_displacement_converges_at_second_order() {
    // 8 -> 16 -> 32 gives two independent rate estimates, which is enough to
    // establish the order. A 64 x 64 rung adds a third at roughly ten times
    // the cost, and it reads 2.00 as well.
    let resolutions = [8usize, 16, 32];
    let errors: Vec<f64> = resolutions
        .iter()
        .map(|&n| solve_mms(&quad4_mesh(n), &exact_2d, &force_2d, None, Some(3)))
        .collect();

    assert_order(
        "Quad4 (bilinear quad, 2-D)",
        &resolutions,
        &errors,
        2.0,
        (1.8, 2.2),
    );
}

/// `Tri3`: linear, theoretical L2 displacement rate 2.
///
/// The same node set as the `Quad4` ladder, so the two are directly
/// comparable: a constant-strain triangle should reach the same rate with a
/// larger constant.
///
/// The load vector is integrated with the 3-point (degree-2) triangle rule
/// rather than the 1-point default. `Tri3` stiffness is exact at one point
/// because `B` is constant, but the consistent load `∫ N_i f` is not: the
/// centroid rule is exact only through degree 1, and a load integrated less
/// accurately than the element interpolates is a variational crime that caps
/// the observed order below the element's own.
#[test]
fn tri3_displacement_converges_at_second_order() {
    let resolutions = [8usize, 16, 32];
    let errors: Vec<f64> = resolutions
        .iter()
        .map(|&n| solve_mms(&tri3_mesh(n), &exact_2d, &force_2d, Some(2), Some(2)))
        .collect();

    assert_order(
        "Tri3 (linear triangle, 2-D)",
        &resolutions,
        &errors,
        2.0,
        (1.8, 2.2),
    );
}

/// `Quad8`: quadratic serendipity, theoretical L2 displacement rate 3.
///
/// This is the rung that a blanket "second order" assertion would not catch.
/// The serendipity space is complete to degree 2, so the L2 displacement error
/// falls as `h^3`; an eight-node element that had, say, mis-ordered its
/// mid-edge nodes would still converge — at rate 2, or worse — and every
/// symmetry and patch check would still pass.
#[test]
fn quad8_displacement_converges_at_third_order() {
    // Third order burns through the error fast: by n = 16 the discretisation
    // error is near 1e-6, so a fourth rung would start competing with the
    // solver's own round-off rather than measuring the element.
    let resolutions = [4usize, 8, 16];
    let errors: Vec<f64> = resolutions
        .iter()
        .map(|&n| solve_mms(&quad8_mesh(n), &exact_2d, &force_2d, Some(4), Some(4)))
        .collect();

    assert_order(
        "Quad8 (serendipity quad, 2-D)",
        &resolutions,
        &errors,
        3.0,
        (2.7, 3.3),
    );
}

/// `Hex8`: trilinear, theoretical L2 displacement rate 2, in three dimensions.
///
/// This exercises the parts of the chain the 2-D ladders cannot reach: the
/// 6-component strain-displacement matrix, the 6x6 `lambda`/`mu` constitutive
/// matrix with its full off-diagonal block, the 3x3 Jacobian and its inverse,
/// and three displacement components per node through DOF numbering and
/// assembly.
#[test]
fn hex8_displacement_converges_at_second_order() {
    // Four rungs, because the cost is dominated by the last one and the cheap
    // rungs are nearly free. `SolverFactory` picks a **dense** LDL^T
    // factorisation of the free-DOF block, which in three dimensions grows as
    // 3 (n-1)^3: 81 unknowns at n = 4, 1029 at n = 8, 2187 at n = 10 — and
    // 3993 at n = 12, which alone takes 73 s against 6 s for this entire
    // ladder. That cubic wall, not the measurement, is what caps the mesh; the
    // observed order is already within 2% of the theoretical value by n = 10.
    let resolutions = [4usize, 6, 8, 10];
    let errors: Vec<f64> = resolutions
        .iter()
        .map(|&n| solve_mms(&hex8_mesh(n), &exact_3d, &body_force_3d, Some(3), Some(3)))
        .collect();

    assert_order(
        "Hex8 (trilinear hexahedron, 3-D)",
        &resolutions,
        &errors,
        2.0,
        (1.8, 2.2),
    );
}

// ---------------------------------------------------------------------------
// Magnitude checks
// ---------------------------------------------------------------------------

/// The error must actually be small, not merely converging at the right rate.
///
/// A scheme converging at second order toward the wrong answer would satisfy
/// the rate test alone.
#[test]
fn quad4_error_is_small_on_a_fine_mesh() {
    let error = solve_mms(&quad4_mesh(32), &exact_2d, &force_2d, None, Some(3));
    assert!(
        error < 1e-3,
        "L2 error {error:.3e} on a 32x32 mesh is too large for a second-order \
         method on a smooth solution"
    );
}

/// The quadratic element must be materially more accurate than the linear one
/// on the same number of elements.
///
/// Rate alone does not establish this: an element could converge at order 3
/// from a starting error so large that it is still worse than `Quad4` on every
/// mesh anyone would use. On a 8 x 8 grid the serendipity element should be at
/// least an order of magnitude closer.
#[test]
fn quad8_beats_quad4_on_the_same_grid() {
    let serendipity = solve_mms(&quad8_mesh(8), &exact_2d, &force_2d, Some(4), Some(4));
    let bilinear = solve_mms(&quad4_mesh(8), &exact_2d, &force_2d, None, Some(3));
    println!(
        "  Quad8 = {serendipity:.6e}   Quad4 = {bilinear:.6e}   ratio = {:.1}x",
        bilinear / serendipity
    );
    assert!(
        serendipity < bilinear / 10.0,
        "Quad8 L2 error {serendipity:.3e} is not an order of magnitude below \
         Quad4's {bilinear:.3e} on the same 8x8 grid"
    );
}

/// The computed 3-D solution must be near the best its own space allows.
///
/// A rate test alone is satisfied by a scheme converging at second order
/// toward the wrong field. The 2-D ladder pins that down with an absolute
/// threshold, but an absolute threshold in three dimensions is just a number
/// someone believed — this manufactured field is larger and rougher than the
/// 2-D one, and 12 x 12 x 12 is a much coarser mesh than 32 x 32.
///
/// The theoretical statement is stronger and needs no such constant. Céa's
/// lemma plus the Aubin–Nitsche duality argument bound the computed error by a
/// **mesh-independent** multiple of the best approximation in the space, of
/// which the nodal interpolant is a representative:
///
/// ```text
///     ||u - u_h||_L2  <=  C ||u - I_h u||_L2,     C independent of h
/// ```
///
/// So this checks the ratio itself, at two resolutions, and requires that it
/// stay small and stop growing. A solver converging to the wrong field, or one
/// that only limps toward the right one, shows up as a ratio that climbs with
/// refinement; a zeroed element matrix or a dropped Dirichlet coupling shows
/// up as a ratio far above 1 at both.
#[test]
fn hex8_error_stays_near_the_best_its_space_allows() {
    let mut ratios = Vec::new();
    for n in [4usize, 8] {
        let mesh = hex8_mesh(n);
        let computed = solve_mms(&mesh, &exact_3d, &body_force_3d, Some(3), Some(3));
        let best = interpolation_error(&mesh, &exact_3d, Some(3));
        println!(
            "  n = {n:3}   computed = {computed:.6e}   interpolant = {best:.6e}   \
             ratio = {:.3}",
            computed / best
        );
        ratios.push(computed / best);
    }

    for (n, &ratio) in [4usize, 8].iter().zip(&ratios) {
        assert!(
            ratio < 3.0,
            "Hex8 at n = {n}: computed L2 error is {ratio:.2}x the nodal \
             interpolation error. Quasi-optimality bounds this by a constant of \
             order one; a value this large means the solve, not the element \
             space, is the limit"
        );
    }
    assert!(
        ratios[1] < ratios[0] * 1.2,
        "Hex8: the computed-to-interpolant error ratio grew from {:.3} to \
         {:.3} under refinement. That constant is h-independent in theory, so a \
         growing one means the discretisation is losing ground as the mesh is \
         refined",
        ratios[0],
        ratios[1]
    );
}