//! Physical invariants of the element stiffness and mass matrices.
//!
//! These are the cheapest decisive checks that element matrix computation is
//! real. None of them needs a reference table: each is a property the exact
//! matrices must satisfy for any correct implementation, so a failure
//! localises to the quadrature, the Jacobian, the shape function derivatives
//! or the constitutive matrix rather than to an accuracy budget.
//!
//! They exist because `StandardFiniteElement::compute_element_matrices`
//! returned `DMatrix::zeros(..)` for every element in the mesh, which made
//! every global matrix `GlobalAssembler` produced zero, and every analysis
//! built on it a solve of a null system. A zero matrix passes a symmetry
//! check and a "does not crash" check, so the invariants below are chosen to
//! be ones a zero matrix fails.

use approx::assert_relative_eq;
use nalgebra::{DVector, Vector3};
use rtx_fea::elements::StandardFiniteElement;
use rtx_fea::materials::LinearElastic;
use rtx_fea::mesh::ElementType;

const E: f64 = 210e9;
const NU: f64 = 0.3;
const RHO: f64 = 7850.0;

fn steel() -> LinearElastic {
    LinearElastic::new(E, NU).with_density(RHO)
}

/// Unit square Quad4, counter-clockwise from the origin. Area = 1.
fn unit_square() -> StandardFiniteElement {
    StandardFiniteElement::new(
        ElementType::Quad4,
        vec![
            Vector3::new(0.0, 0.0, 0.0),
            Vector3::new(1.0, 0.0, 0.0),
            Vector3::new(1.0, 1.0, 0.0),
            Vector3::new(0.0, 1.0, 0.0),
        ],
    )
}

/// A 2 x 3 rectangle, to catch anything that only works when `det J == 1`.
fn rectangle_2x3() -> StandardFiniteElement {
    StandardFiniteElement::new(
        ElementType::Quad4,
        vec![
            Vector3::new(0.0, 0.0, 0.0),
            Vector3::new(2.0, 0.0, 0.0),
            Vector3::new(2.0, 3.0, 0.0),
            Vector3::new(0.0, 3.0, 0.0),
        ],
    )
}

/// The stub this suite exists to catch returned zeros everywhere.
#[test]
fn element_matrices_are_not_zero() {
    let m = unit_square()
        .compute_element_matrices(&steel(), 0.0)
        .unwrap();

    assert!(
        m.stiffness_matrix.amax() > 0.0,
        "stiffness matrix is entirely zero — element matrix computation is a stub"
    );
    let mass = m
        .mass_matrix
        .as_ref()
        .expect("a mass matrix is required for modal and dynamic analysis");
    assert!(
        mass.amax() > 0.0,
        "mass matrix is entirely zero — element matrix computation is a stub"
    );
}

/// `K` must be symmetric: it is `∫ Bᵀ D B dV` with `D` symmetric.
#[test]
fn stiffness_is_symmetric() {
    let m = rectangle_2x3()
        .compute_element_matrices(&steel(), 0.0)
        .unwrap();
    let k = &m.stiffness_matrix;

    for i in 0..k.nrows() {
        for j in 0..k.ncols() {
            assert_relative_eq!(k[(i, j)], k[(j, i)], epsilon = 1e-6 * k.amax());
        }
    }
}

/// A rigid-body translation stores no strain energy, so `K t = 0`.
///
/// This is the single most informative check on the strain-displacement
/// matrix `B`: it fails for a sign error, a mis-transformed derivative, or a
/// wrong DOF ordering, none of which a symmetry check detects. It is the
/// element-level half of the patch test.
#[test]
fn rigid_translation_produces_no_internal_force() {
    let element = rectangle_2x3();
    let m = element.compute_element_matrices(&steel(), 0.0).unwrap();
    let k = &m.stiffness_matrix;

    let num_nodes = 4;
    let dim = 2;
    assert_eq!(k.nrows(), num_nodes * dim);

    // One translation per spatial direction.
    for d in 0..dim {
        let mut t = DVector::zeros(num_nodes * dim);
        for node in 0..num_nodes {
            t[node * dim + d] = 1.0;
        }

        let f = k * &t;
        assert!(
            f.amax() < 1e-6 * k.amax(),
            "translating the element in direction {d} produced internal force {:.3e} \
             against a stiffness scale of {:.3e}",
            f.amax(),
            k.amax()
        );
    }
}

/// An unconstrained plane element has exactly three rigid-body modes — two
/// translations and one rotation — so `K` has a three-dimensional null space.
///
/// Fewer means the element is spuriously stiff; more means it is rank
/// deficient and admits a zero-energy deformation (hourglassing), which shows
/// up in a real analysis as a mode shape that is pure noise.
#[test]
fn stiffness_has_exactly_three_rigid_body_modes_in_2d() {
    let m = rectangle_2x3()
        .compute_element_matrices(&steel(), 0.0)
        .unwrap();
    let k = m.stiffness_matrix;

    let eigenvalues = k.clone().symmetric_eigenvalues();
    let scale = k.amax();
    let num_zero = eigenvalues
        .iter()
        .filter(|&&e| e.abs() < 1e-9 * scale)
        .count();

    assert_eq!(
        num_zero,
        3,
        "expected 3 rigid-body modes (2 translations + 1 rotation), found {num_zero}; \
         eigenvalues (scaled): {:?}",
        eigenvalues.iter().map(|e| e / scale).collect::<Vec<_>>()
    );
}

/// The consistent mass matrix integrates to the element's total mass.
///
/// `Σᵢⱼ Mᵢⱼ = ∫ ρ (Σᵢ Nᵢ)(Σⱼ Nⱼ) dV = ∫ ρ dV = ρV`, using the partition of
/// unity. This is decisive for the quadrature rule and the Jacobian
/// determinant together: get either wrong and the total mass is wrong by
/// exactly the factor of the error.
#[test]
fn consistent_mass_integrates_to_rho_times_volume() {
    for (element, area) in [(unit_square(), 1.0), (rectangle_2x3(), 6.0)] {
        let m = element.compute_element_matrices(&steel(), 0.0).unwrap();
        let mass = m.mass_matrix.expect("mass matrix required");

        // Vector-valued mass: each spatial direction carries the full mass,
        // so the whole matrix sums to `dim * rho * V`.
        let dim = 2;
        let total: f64 = mass.iter().sum();
        assert_relative_eq!(total, dim as f64 * RHO * area, max_relative = 1e-9);
    }
}

/// The mass matrix must be positive definite.
///
/// It is the metric in `K φ = λ M φ`; if it is not positive definite the
/// eigenproblem has no real spectrum and the Cholesky reduction in the
/// eigensolver fails outright.
#[test]
fn mass_matrix_is_positive_definite() {
    let m = rectangle_2x3()
        .compute_element_matrices(&steel(), 0.0)
        .unwrap();
    let mass = m.mass_matrix.expect("mass matrix required");

    let eigenvalues = mass.clone().symmetric_eigenvalues();
    let min = eigenvalues.iter().cloned().fold(f64::INFINITY, f64::min);
    assert!(
        min > 0.0,
        "mass matrix has a non-positive eigenvalue {min:.3e}; \
         eigenvalues: {eigenvalues:?}"
    );
}

/// Mass scales linearly with density, stiffness with elastic modulus, and
/// neither picks up the other's property.
///
/// A single hard-coded default leaking into the computation would break this,
/// and would otherwise be invisible in a suite that only ever uses one
/// material.
#[test]
fn matrices_scale_with_the_material_they_are_given() {
    let base = steel();
    let stiffer = LinearElastic::new(2.0 * E, NU).with_density(RHO);
    let denser = LinearElastic::new(E, NU).with_density(3.0 * RHO);

    let element = rectangle_2x3();
    let m0 = element.compute_element_matrices(&base, 0.0).unwrap();
    let m_stiff = element.compute_element_matrices(&stiffer, 0.0).unwrap();
    let m_dense = element.compute_element_matrices(&denser, 0.0).unwrap();

    assert_relative_eq!(
        m_stiff.stiffness_matrix.amax(),
        2.0 * m0.stiffness_matrix.amax(),
        max_relative = 1e-12
    );
    assert_relative_eq!(
        m_dense.mass_matrix.as_ref().unwrap().amax(),
        3.0 * m0.mass_matrix.as_ref().unwrap().amax(),
        max_relative = 1e-12
    );

    // Doubling E must not change the mass; tripling ρ must not change K.
    assert_relative_eq!(
        m_stiff.mass_matrix.as_ref().unwrap().amax(),
        m0.mass_matrix.as_ref().unwrap().amax(),
        max_relative = 1e-12
    );
    assert_relative_eq!(
        m_dense.stiffness_matrix.amax(),
        m0.stiffness_matrix.amax(),
        max_relative = 1e-12
    );
}

/// The same invariants must hold in 3D, where the null space is six
/// dimensional — three translations and three rotations.
#[test]
fn hexahedron_satisfies_the_same_invariants() {
    let element = StandardFiniteElement::new(
        ElementType::Hex8,
        vec![
            Vector3::new(0.0, 0.0, 0.0),
            Vector3::new(2.0, 0.0, 0.0),
            Vector3::new(2.0, 1.0, 0.0),
            Vector3::new(0.0, 1.0, 0.0),
            Vector3::new(0.0, 0.0, 3.0),
            Vector3::new(2.0, 0.0, 3.0),
            Vector3::new(2.0, 1.0, 3.0),
            Vector3::new(0.0, 1.0, 3.0),
        ],
    );
    let volume = 2.0 * 1.0 * 3.0;

    let m = element.compute_element_matrices(&steel(), 0.0).unwrap();
    let k = m.stiffness_matrix;
    let mass = m.mass_matrix.expect("mass matrix required");

    assert_eq!(k.nrows(), 8 * 3);

    // Rigid translations.
    for d in 0..3 {
        let mut t = DVector::zeros(24);
        for node in 0..8 {
            t[node * 3 + d] = 1.0;
        }
        let f = &k * &t;
        assert!(
            f.amax() < 1e-6 * k.amax(),
            "3-D rigid translation in direction {d} produced force {:.3e}",
            f.amax()
        );
    }

    let scale = k.amax();
    let num_zero = k
        .clone()
        .symmetric_eigenvalues()
        .iter()
        .filter(|&&e| e.abs() < 1e-9 * scale)
        .count();
    assert_eq!(
        num_zero, 6,
        "expected 6 rigid-body modes in 3-D, found {num_zero}"
    );

    let total: f64 = mass.iter().sum();
    assert_relative_eq!(total, 3.0 * RHO * volume, max_relative = 1e-9);
}