//! Sample Data for StructuralPINN Demo
//!
//! This module provides sample geometries, materials, boundary conditions,
//! and loads for testing and demonstration purposes.

use structural_shared::{
    sample_plate_with_hole as create_plate_with_hole, sample_rectangular_mesh, BoundaryCondition,
    Element2D, Geometry2D, Load, Material, Point2D,
};

// ============================================================================
// Sample Materials
// ============================================================================

/// Structural steel material (AISI 1020).
pub fn steel_material() -> Material {
    Material::steel()
}

/// Aluminum alloy material (6061-T6).
pub fn aluminum_material() -> Material {
    Material::aluminum()
}

/// Titanium alloy material (Ti-6Al-4V).
pub fn titanium_material() -> Material {
    Material::titanium()
}

/// Concrete material.
pub fn concrete_material() -> Material {
    Material::concrete()
}

/// Custom material with specified properties.
pub fn custom_material(
    youngs_modulus: f64,
    poisson_ratio: f64,
    density: f64,
    yield_stress: f64,
) -> Material {
    Material::new(youngs_modulus, poisson_ratio, density, yield_stress)
}

// ============================================================================
// Sample Geometries
// ============================================================================

/// Create a cantilever beam geometry.
/// Beam is 1.0m long and 0.1m tall.
pub fn cantilever_beam() -> Geometry2D {
    sample_rectangular_mesh(1.0, 0.1, 20, 4)
}

/// Create a cantilever beam with finer mesh.
pub fn cantilever_beam_fine() -> Geometry2D {
    sample_rectangular_mesh(1.0, 0.1, 40, 8)
}

/// Create a simply supported beam geometry.
/// Beam is 2.0m long and 0.1m tall.
pub fn simply_supported_beam() -> Geometry2D {
    sample_rectangular_mesh(2.0, 0.1, 40, 4)
}

/// Create a plate with a circular hole.
/// Plate is 1.0m x 1.0m with a 0.2m radius hole.
pub fn plate_with_hole() -> Geometry2D {
    create_plate_with_hole(1.0, 1.0, 0.2, 20, 20)
}

/// Create a rectangular plate for uniform loading.
pub fn rectangular_plate(width: f64, height: f64) -> Geometry2D {
    sample_rectangular_mesh(width, height, 20, 20)
}

/// Create an L-shaped bracket geometry.
pub fn l_bracket() -> Geometry2D {
    let mut vertices = Vec::new();
    let mut elements = Vec::new();

    // L-bracket dimensions: 1.0 x 1.0 overall, 0.3m thickness
    let thickness = 0.3;
    let length = 1.0;
    let height = 1.0;
    let n = 10;

    // Bottom horizontal part
    for j in 0..=n / 3 {
        for i in 0..=n {
            let x = i as f64 / n as f64 * length;
            let y = j as f64 / (n / 3) as f64 * thickness;
            vertices.push(Point2D::new(x, y));
        }
    }

    // Vertical part
    let offset = vertices.len();
    for j in 1..=n {
        for i in 0..=n / 3 {
            let x = i as f64 / (n / 3) as f64 * thickness;
            let y = thickness + j as f64 / n as f64 * (height - thickness);
            vertices.push(Point2D::new(x, y));
        }
    }

    // Create elements for bottom part
    let cols1 = n + 1;
    let rows1 = n / 3;
    for j in 0..rows1 {
        for i in 0..n {
            let n0 = j * cols1 + i;
            let n1 = n0 + 1;
            let n2 = n0 + cols1;
            let n3 = n2 + 1;
            elements.push(Element2D::Triangle([n0, n1, n2]));
            elements.push(Element2D::Triangle([n1, n3, n2]));
        }
    }

    // Create elements for vertical part (simplified)
    let cols2 = n / 3 + 1;
    for j in 0..n - 1 {
        for i in 0..n / 3 {
            let n0 = offset + j * cols2 + i;
            let n1 = n0 + 1;
            let n2 = n0 + cols2;
            let n3 = n2 + 1;
            if n3 < vertices.len() {
                elements.push(Element2D::Triangle([n0, n1, n2]));
                elements.push(Element2D::Triangle([n1, n3, n2]));
            }
        }
    }

    Geometry2D::new(vertices, elements)
}

/// Create a thin-walled cylinder approximation (2D cross-section).
pub fn pressure_vessel_2d() -> Geometry2D {
    // Annular cross-section
    let inner_radius = 0.4;
    let outer_radius = 0.5;
    let n_radial = 4;
    let n_circumferential = 32;

    let mut vertices = Vec::new();
    let mut elements = Vec::new();

    // Generate vertices
    for i in 0..=n_radial {
        let r = inner_radius + (i as f64 / n_radial as f64) * (outer_radius - inner_radius);
        for j in 0..n_circumferential {
            let theta = 2.0 * std::f64::consts::PI * j as f64 / n_circumferential as f64;
            vertices.push(Point2D::new(r * theta.cos(), r * theta.sin()));
        }
    }

    // Generate elements
    for i in 0..n_radial {
        for j in 0..n_circumferential {
            let n0 = i * n_circumferential + j;
            let n1 = i * n_circumferential + (j + 1) % n_circumferential;
            let n2 = (i + 1) * n_circumferential + j;
            let n3 = (i + 1) * n_circumferential + (j + 1) % n_circumferential;
            elements.push(Element2D::Triangle([n0, n1, n2]));
            elements.push(Element2D::Triangle([n1, n3, n2]));
        }
    }

    Geometry2D::new(vertices, elements)
}

// ============================================================================
// Boundary Conditions and Loads
// ============================================================================

/// Cantilever beam boundary conditions and loads.
/// Fixed at left end, point load at right end.
pub fn cantilever_beam_bcs_loads() -> (Vec<BoundaryCondition>, Vec<Load>) {
    let geometry = cantilever_beam();
    let n_vertices = geometry.num_vertices();
    let nx = 21; // Number of vertices in x-direction

    // Fixed nodes at x = 0 (left edge)
    let fixed_nodes: Vec<usize> = (0..n_vertices).filter(|&i| i % nx == 0).collect();

    // Loaded nodes at x = 1.0 (right edge)
    let loaded_nodes: Vec<usize> = (0..n_vertices).filter(|&i| i % nx == nx - 1).collect();

    let bcs = vec![BoundaryCondition::fixed(fixed_nodes, 2)];

    let loads = vec![Load::Point {
        nodes: loaded_nodes,
        force: vec![0.0, -10000.0], // 10 kN downward
    }];

    (bcs, loads)
}

/// Simply supported beam boundary conditions and loads.
/// Pin at left, roller at right, distributed load on top.
pub fn simply_supported_beam_bcs_loads() -> (Vec<BoundaryCondition>, Vec<Load>) {
    let geometry = simply_supported_beam();
    let n_vertices = geometry.num_vertices();
    let nx = 41; // Number of vertices in x-direction
    let ny = 5; // Number of vertices in y-direction

    // Pin at left (x = 0, y = 0): fixed in x and y
    let left_pin = vec![0];

    // Roller at right (x = 2.0, y = 0): fixed in y only
    let right_roller = vec![nx - 1];

    let bcs = vec![
        BoundaryCondition::fixed(left_pin, 2),
        BoundaryCondition::Dirichlet {
            nodes: right_roller,
            displacement: vec![0.0, 0.0],
            constrained: vec![false, true], // Only y constrained
        },
    ];

    // Distributed load on top surface
    let top_nodes: Vec<usize> = (0..n_vertices).filter(|&i| i / nx == ny - 1).collect();

    let loads = vec![Load::Distributed {
        elements: top_nodes,
        intensity: vec![0.0, -5000.0], // 5 kN/m downward
        direction: vec![0.0, -1.0],
    }];

    (bcs, loads)
}

/// Plate with hole boundary conditions and loads.
/// Fixed at left edge, tension at right edge.
pub fn plate_with_hole_bcs_loads() -> (Vec<BoundaryCondition>, Vec<Load>) {
    let geometry = plate_with_hole();
    let n_vertices = geometry.num_vertices();

    // Fixed nodes near x = 0
    let fixed_nodes: Vec<usize> = (0..n_vertices)
        .filter(|&i| geometry.vertices[i].x < 0.1)
        .collect();

    // Loaded nodes near x = 1.0
    let loaded_nodes: Vec<usize> = (0..n_vertices)
        .filter(|&i| geometry.vertices[i].x > 0.9)
        .collect();

    let bcs = vec![BoundaryCondition::fixed(fixed_nodes, 2)];

    let loads = vec![Load::Point {
        nodes: loaded_nodes,
        force: vec![50000.0, 0.0], // 50 kN tension
    }];

    (bcs, loads)
}

/// L-bracket boundary conditions and loads.
/// Fixed at bottom, load at top of vertical section.
pub fn l_bracket_bcs_loads() -> (Vec<BoundaryCondition>, Vec<Load>) {
    let geometry = l_bracket();
    let n_vertices = geometry.num_vertices();

    // Fixed nodes at bottom (y = 0)
    let fixed_nodes: Vec<usize> = (0..n_vertices)
        .filter(|&i| geometry.vertices[i].y < 0.05)
        .collect();

    // Loaded nodes at top of vertical section
    let loaded_nodes: Vec<usize> = (0..n_vertices)
        .filter(|&i| geometry.vertices[i].y > 0.9)
        .collect();

    let bcs = vec![BoundaryCondition::fixed(fixed_nodes, 2)];

    let loads = vec![Load::Point {
        nodes: loaded_nodes,
        force: vec![20000.0, 0.0], // 20 kN horizontal
    }];

    (bcs, loads)
}

/// Pressure vessel boundary conditions and loads.
/// Internal pressure with symmetric boundary conditions.
pub fn pressure_vessel_bcs_loads() -> (Vec<BoundaryCondition>, Vec<Load>) {
    let geometry = pressure_vessel_2d();
    let n_vertices = geometry.num_vertices();
    let n_circumferential = 32;

    // Inner surface nodes for pressure
    let inner_nodes: Vec<usize> = (0..n_circumferential).collect();

    // Symmetry at x = 0 plane (nodes near theta = pi/2 or 3pi/2)
    let sym_nodes: Vec<usize> = (0..n_vertices)
        .filter(|&i| geometry.vertices[i].x.abs() < 0.05)
        .collect();

    let bcs = vec![BoundaryCondition::Symmetry {
        nodes: sym_nodes,
        normal: vec![1.0, 0.0],
    }];

    let loads = vec![Load::Pressure {
        faces: inner_nodes,
        magnitude: 1.0e6, // 1 MPa internal pressure
    }];

    (bcs, loads)
}

/// Gravity load on a structure.
pub fn gravity_load() -> Load {
    Load::gravity()
}

/// Thermal load with uniform temperature change.
pub fn thermal_load(num_nodes: usize, delta_t: f64) -> Load {
    Load::Thermal {
        temperature: vec![delta_t; num_nodes],
        reference_temp: 293.15, // 20 C
    }
}

// ============================================================================
// Complete Problem Configurations
// ============================================================================

/// Complete cantilever beam problem.
pub fn cantilever_problem() -> (Geometry2D, Material, Vec<BoundaryCondition>, Vec<Load>) {
    let geometry = cantilever_beam();
    let material = steel_material();
    let (bcs, loads) = cantilever_beam_bcs_loads();
    (geometry, material, bcs, loads)
}

/// Complete simply supported beam problem.
pub fn simply_supported_problem() -> (Geometry2D, Material, Vec<BoundaryCondition>, Vec<Load>) {
    let geometry = simply_supported_beam();
    let material = steel_material();
    let (bcs, loads) = simply_supported_beam_bcs_loads();
    (geometry, material, bcs, loads)
}

/// Complete plate with hole problem.
pub fn plate_with_hole_problem() -> (Geometry2D, Material, Vec<BoundaryCondition>, Vec<Load>) {
    let geometry = plate_with_hole();
    let material = aluminum_material();
    let (bcs, loads) = plate_with_hole_bcs_loads();
    (geometry, material, bcs, loads)
}

/// Complete pressure vessel problem.
pub fn pressure_vessel_problem() -> (Geometry2D, Material, Vec<BoundaryCondition>, Vec<Load>) {
    let geometry = pressure_vessel_2d();
    let material = steel_material();
    let (bcs, loads) = pressure_vessel_bcs_loads();
    (geometry, material, bcs, loads)
}

// ============================================================================
// Analytical Solutions for Verification
// ============================================================================

/// Analytical maximum deflection for cantilever beam with end load.
/// delta_max = P * L^3 / (3 * E * I)
pub fn cantilever_analytical_deflection(
    load: f64,   // Point load (N)
    length: f64, // Beam length (m)
    e: f64,      // Young's modulus (Pa)
    i: f64,      // Second moment of area (m^4)
) -> f64 {
    load * length.powi(3) / (3.0 * e * i)
}

/// Analytical maximum stress for cantilever beam with end load.
/// sigma_max = M * c / I = P * L * c / I
pub fn cantilever_analytical_stress(
    load: f64,   // Point load (N)
    length: f64, // Beam length (m)
    height: f64, // Beam height (m)
    i: f64,      // Second moment of area (m^4)
) -> f64 {
    let moment = load * length;
    let c = height / 2.0;
    moment * c / i
}

/// Second moment of area for rectangular cross-section.
/// I = b * h^3 / 12
pub fn rectangular_moment_of_inertia(width: f64, height: f64) -> f64 {
    width * height.powi(3) / 12.0
}

/// Analytical hoop stress in thin-walled pressure vessel.
/// sigma_hoop = p * r / t
pub fn pressure_vessel_hoop_stress(pressure: f64, radius: f64, thickness: f64) -> f64 {
    pressure * radius / thickness
}

/// Stress concentration factor for plate with circular hole.
/// K_t approximately 3.0 for small hole in wide plate.
pub fn plate_with_hole_scf(hole_radius: f64, plate_width: f64) -> f64 {
    // Peterson's approximation
    let r_w = hole_radius / plate_width;
    if r_w < 0.5 {
        3.0 - 3.13 * r_w + 3.66 * r_w.powi(2) - 1.53 * r_w.powi(3)
    } else {
        3.0 // Approximate
    }
}

// ============================================================================
// Tests
// ============================================================================

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_steel_material() {
        let steel = steel_material();
        assert!((steel.youngs_modulus - 200.0e9).abs() < 1.0);
    }

    #[test]
    fn test_aluminum_material() {
        let al = aluminum_material();
        assert!(al.youngs_modulus < steel_material().youngs_modulus);
    }

    #[test]
    fn test_cantilever_beam() {
        let geom = cantilever_beam();
        assert!(geom.num_vertices() > 0);
        assert!(geom.num_elements() > 0);
        // Beam should be longer than tall
        assert!(geom.bounds.width() > geom.bounds.height());
    }

    #[test]
    fn test_simply_supported_beam() {
        let geom = simply_supported_beam();
        assert!(geom.num_vertices() > 0);
        assert!((geom.bounds.width() - 2.0).abs() < 0.01);
    }

    #[test]
    fn test_plate_with_hole() {
        let geom = plate_with_hole();
        assert!(geom.num_vertices() > 0);
    }

    #[test]
    fn test_l_bracket() {
        let geom = l_bracket();
        assert!(geom.num_vertices() > 0);
    }

    #[test]
    fn test_pressure_vessel() {
        let geom = pressure_vessel_2d();
        assert!(geom.num_vertices() > 0);
        // Check circular shape - centroid should be near origin
        let cx: f64 = geom.vertices.iter().map(|v| v.x).sum::<f64>() / geom.num_vertices() as f64;
        let cy: f64 = geom.vertices.iter().map(|v| v.y).sum::<f64>() / geom.num_vertices() as f64;
        assert!(cx.abs() < 0.1);
        assert!(cy.abs() < 0.1);
    }

    #[test]
    fn test_cantilever_bcs() {
        let (bcs, loads) = cantilever_beam_bcs_loads();
        assert!(!bcs.is_empty());
        assert!(!loads.is_empty());
    }

    #[test]
    fn test_simply_supported_bcs() {
        let (bcs, loads) = simply_supported_beam_bcs_loads();
        assert!(bcs.len() >= 2); // Pin and roller
        assert!(!loads.is_empty());
    }

    #[test]
    fn test_cantilever_analytical() {
        let p = 10000.0; // 10 kN
        let l = 1.0; // 1 m
        let e = 200.0e9; // Steel
        let h = 0.1; // 10 cm
        let b = 0.01; // 1 cm (assuming unit depth)
        let i = rectangular_moment_of_inertia(b, h);

        let deflection = cantilever_analytical_deflection(p, l, e, i);
        assert!(deflection > 0.0);
        assert!(deflection < 1.0); // Reasonable deflection

        let stress = cantilever_analytical_stress(p, l, h, i);
        assert!(stress > 0.0);
    }

    #[test]
    fn test_moment_of_inertia() {
        let b = 0.1; // 10 cm
        let h = 0.2; // 20 cm
        let i = rectangular_moment_of_inertia(b, h);
        // I = 0.1 * 0.2^3 / 12 = 6.67e-5 m^4
        assert!((i - 6.67e-5).abs() < 1e-6);
    }

    #[test]
    fn test_hoop_stress() {
        let p = 1.0e6; // 1 MPa
        let r = 0.45; // Mean radius
        let t = 0.1; // Wall thickness
        let stress = pressure_vessel_hoop_stress(p, r, t);
        // sigma = 1e6 * 0.45 / 0.1 = 4.5 MPa
        assert!((stress - 4.5e6).abs() < 1e4);
    }

    #[test]
    fn test_scf_small_hole() {
        let r = 0.1;
        let w = 1.0;
        let scf = plate_with_hole_scf(r, w);
        // For small hole (r/w = 0.1), SCF from Peterson's formula is ~2.7
        assert!(scf > 2.5 && scf < 3.1);
    }

    #[test]
    fn test_complete_problems() {
        let (geom, mat, bcs, loads) = cantilever_problem();
        assert!(geom.num_vertices() > 0);
        assert!(mat.youngs_modulus > 0.0);
        assert!(!bcs.is_empty());
        assert!(!loads.is_empty());
    }

    #[test]
    fn test_gravity_load() {
        let load = gravity_load();
        match load {
            Load::Body { acceleration } => {
                assert!((acceleration[1] + 9.81).abs() < 0.01);
            }
            _ => panic!("Expected Body load"),
        }
    }

    #[test]
    fn test_thermal_load() {
        let load = thermal_load(100, 50.0);
        match load {
            Load::Thermal {
                temperature,
                reference_temp,
            } => {
                assert_eq!(temperature.len(), 100);
                assert!((temperature[0] - 50.0).abs() < 0.01);
                assert!((reference_temp - 293.15).abs() < 0.01);
            }
            _ => panic!("Expected Thermal load"),
        }
    }
}