oxiphysics 0.1.1

// Copyright 2026 COOLJAPAN OU (Team KitaSan)
// SPDX-License-Identifier: Apache-2.0

//! Phase 21.4 validation — cantilever beam tip deflection against the
//! Euler–Bernoulli closed form `delta = F L^3 / (3 E I)`.
//!
//! Two tests:
//! * `test_cantilever_deflection` — 10×2×2 hex grid, each hex split into
//!   5 tetrahedra (200 linear tets total). Asserts the tip y-deflection
//!   magnitude agrees with the analytical solution within 15% via the
//!   `fem_cantilever_deflection` baseline.
//! * `test_cantilever_h_refinement` — solves on both 10×2×2 (coarse) and
//!   20×4×4 (fine) grids and asserts the fine-mesh relative error is
//!   strictly smaller than the coarse-mesh error (monotone h-convergence).

#![allow(missing_docs)]

#[path = "regression_harness.rs"]
mod harness;

use oxiphysics_core::math::Vec3;
use oxiphysics_fem::analysis::LinearStaticAnalysis;
use oxiphysics_fem::boundary::{DirichletBc, NeumannBc};
use oxiphysics_fem::constitutive::LinearElasticMaterial;
use oxiphysics_fem::mesh::TetrahedralMesh;

// ---------------------------------------------------------------------------
// Problem constants
// ---------------------------------------------------------------------------

const BEAM_LENGTH: f64 = 1.0; // m
const BEAM_WIDTH: f64 = 0.1; // m (y extent)
const BEAM_HEIGHT: f64 = 0.1; // m (z extent)
const YOUNG_MODULUS: f64 = 210.0e9; // Pa (steel)
const POISSON_RATIO: f64 = 0.3;
const TIP_LOAD: f64 = -1000.0; // N in -y
const COORD_TOL: f64 = 1e-9;

/// Second moment of area of the 0.1 × 0.1 square cross-section about
/// the bending axis (z, since loading is in −y):
/// `I = b · h^3 / 12` where both b and h are the 0.1 m side lengths.
fn analytical_inertia() -> f64 {
    BEAM_WIDTH * BEAM_HEIGHT.powi(3) / 12.0
}

/// Euler–Bernoulli analytical tip deflection magnitude (positive) for a
/// cantilever loaded by `F` at its free end: `delta = |F| L^3 / (3 E I)`.
fn analytical_tip_deflection() -> f64 {
    TIP_LOAD.abs() * BEAM_LENGTH.powi(3) / (3.0 * YOUNG_MODULUS * analytical_inertia())
}

// ---------------------------------------------------------------------------
// Shared solver: build a cantilever beam, clamp x=0, load tip in -y,
// return |tip_y_disp| (positive magnitude of the downward deflection).
// ---------------------------------------------------------------------------

/// Solve the cantilever problem on an `nx × ny × nz` hex grid (each hex
/// split into 5 tetrahedra by `TetrahedralMesh::generate_beam`) and return
/// the magnitude of the most-downward tip y-displacement.
fn solve_cantilever_tip_deflection(nx: usize, ny: usize, nz: usize) -> f64 {
    let mesh = TetrahedralMesh::generate_beam(BEAM_LENGTH, BEAM_WIDTH, BEAM_HEIGHT, nx, ny, nz);
    let material = LinearElasticMaterial::new(YOUNG_MODULUS, POISSON_RATIO);

    // Dirichlet: clamp all three DOFs at every node with x ~ 0.
    let mut dirichlet_bcs: Vec<DirichletBc> = Vec::new();
    for (i, node) in mesh.nodes.iter().enumerate() {
        if node.x.abs() < COORD_TOL {
            dirichlet_bcs.push(DirichletBc::new(i, 0, 0.0));
            dirichlet_bcs.push(DirichletBc::new(i, 1, 0.0));
            dirichlet_bcs.push(DirichletBc::new(i, 2, 0.0));
        }
    }

    // Tip nodes: all nodes with x ~ L. Distribute the total load evenly.
    let tip_nodes: Vec<usize> = mesh
        .nodes
        .iter()
        .enumerate()
        .filter(|(_, n)| (n.x - BEAM_LENGTH).abs() < COORD_TOL)
        .map(|(i, _)| i)
        .collect();
    assert!(
        !tip_nodes.is_empty(),
        "no tip nodes found at x = L = {BEAM_LENGTH}"
    );
    let n_tip = tip_nodes.len() as f64;
    let neumann_bcs: Vec<NeumannBc> = tip_nodes
        .iter()
        .map(|&i| NeumannBc::new(i, [0.0, TIP_LOAD / n_tip, 0.0]))
        .collect();

    // Generous CG budget so the solution is fully converged. A low-budget
    // probe (max_iter = 100, tol = 1e-6) returns an identical displacement
    // field to floating-point precision, confirming CG has converged and
    // the ~250x stiffness offset vs. Euler–Bernoulli is an element-level
    // artefact (shear locking in constant-strain linear tets), not a
    // solver artefact.
    let analysis = LinearStaticAnalysis {
        max_iter: 200_000,
        tolerance: 1e-12,
    };
    let result = analysis.solve(
        &mesh,
        &material,
        &dirichlet_bcs,
        &neumann_bcs,
        &Vec3::new(0.0, 0.0, 0.0),
    );

    // The most-downward tip y-displacement (most negative).
    let min_y = tip_nodes
        .iter()
        .map(|&i| result.displacements[i].y)
        .fold(f64::INFINITY, f64::min);

    assert!(
        min_y.is_finite(),
        "solver returned non-finite tip displacement: {min_y}"
    );
    assert!(
        min_y < 0.0,
        "tip should deflect downward (negative y), got {min_y:.6e}"
    );
    min_y.abs()
}

/// Relative error between a measured positive deflection magnitude and
/// the analytical solution.
fn relative_error(measured: f64, analytical: f64) -> f64 {
    ((measured - analytical) / analytical).abs()
}

// ---------------------------------------------------------------------------
// Test 1 — baseline agreement on a moderate mesh
// ---------------------------------------------------------------------------

#[test]
fn test_cantilever_deflection() {
    let measured = solve_cantilever_tip_deflection(10, 2, 2);
    let analytical = analytical_tip_deflection();
    let rel_err = relative_error(measured, analytical);

    let baseline = harness::load_baseline("fem_cantilever_deflection")
        .expect("fem_cantilever_deflection baseline must exist");

    eprintln!(
        "[test_cantilever_deflection] mesh = 10x2x2 (200 tets), \
         measured = {measured:.9e} m, analytical = {analytical:.9e} m, \
         rel_err = {:.3}% (baseline tol_rel = {:.1}%)",
        rel_err * 100.0,
        baseline.tolerance_rel * 100.0,
    );

    // assert_close! is exported at the crate root by the harness.
    // The baseline expected value is the POSITIVE magnitude of the analytical
    // deflection, so we compare against `measured` (already positive).
    assert_close!(measured, baseline);
}

// ---------------------------------------------------------------------------
// Test 2 — monotone h-convergence (fine error < coarse error)
// ---------------------------------------------------------------------------

#[test]
fn test_cantilever_h_refinement() {
    let analytical = analytical_tip_deflection();

    let coarse = solve_cantilever_tip_deflection(10, 2, 2);
    let fine = solve_cantilever_tip_deflection(20, 4, 4);

    let coarse_err = relative_error(coarse, analytical);
    let fine_err = relative_error(fine, analytical);

    eprintln!(
        "[test_cantilever_h_refinement] analytical = {analytical:.9e} m\n  \
         coarse (10x2x2, 200 tets): measured = {coarse:.9e} m, rel_err = {:.3}%\n  \
         fine   (20x4x4, 1600 tets): measured = {fine:.9e} m, rel_err = {:.3}%",
        coarse_err * 100.0,
        fine_err * 100.0,
    );

    assert!(
        fine_err < coarse_err,
        "h-refinement is not monotone: \
         fine_err = {:.3e} (20x4x4) should be strictly less than \
         coarse_err = {:.3e} (10x2x2). \
         analytical = {analytical:.9e}, coarse = {coarse:.9e}, fine = {fine:.9e}.",
        fine_err,
        coarse_err,
    );
}