rustytorch/crates/specialized/rtx-fea/examples/cantilever_beam.rs

// Copyright (c) 2024 RustyTorch++ Team
// Licensed under the Apache License, Version 2.0

//! Comprehensive cantilever beam analysis example.
//!
//! This example demonstrates a complete finite element analysis workflow
//! using the RTX-FEA crate, including:
//! - Mesh generation
//! - Material assignment
//! - Boundary condition application
//! - Static linear analysis
//! - Results visualization

use rtx_fea::prelude::*;
use rtx_fea::{
    analysis::{AnalysisConfig, StaticLinearAnalysis},
    assembly::DofComponent,
    boundary::{BoundaryConditionSet, DirichletPatterns, NeumannPatterns},
    materials::{LinearElastic, MaterialDatabase},
    mesh::MaterialId,
    utils::{BenchmarkUtils, StringUtils, VtkWriter},
};

fn main() -> FeaResult<()> {
    println!("RTX-FEA Cantilever Beam Analysis Example");
    println!("=========================================");

    // Step 1: Create the beam geometry and mesh
    println!("\n1. Creating beam mesh...");
    let beam_length = 10.0; // meters
    let beam_height = 1.0; // meters
    let nx_elements = 21; // Number of nodes in x-direction (20 elements)
    let ny_elements = 5; // Number of nodes in y-direction (4 elements)

    // Use Mesh::generate_rectangle instead of Rectangle::generate_quad_mesh
    let mesh = Mesh::generate_rectangle(beam_length, beam_height, nx_elements, ny_elements)?;

    println!(
        "  - Mesh created with {} nodes and {} elements",
        mesh.num_nodes(),
        mesh.num_elements()
    );

    // Step 2: Set up materials
    println!("\n2. Setting up materials...");
    let mut materials = MaterialDatabase::new();

    // Steel properties
    let elastic_modulus = 200e9; // Pa (200 GPa)
    let poisson_ratio = 0.3;
    let _density = 7850.0; // kg/m³ (not used in linear elastic analysis)

    let steel = LinearElastic::new(elastic_modulus, poisson_ratio);
    let steel_id = MaterialId(0);
    materials.add_material(steel_id, steel, Some("Steel".to_string()));

    println!(
        "  - Added steel material: E = {:.0} GPa, ν = {:.2}",
        elastic_modulus / 1e9,
        poisson_ratio
    );

    // Step 3: Set up boundary conditions
    println!("\n3. Setting up boundary conditions...");
    let mut boundary_conditions = BoundaryConditionSet::new();

    // Find nodes at the left end (x = 0) for fixed support
    let mut fixed_nodes = Vec::new();
    for (&node_id, node) in &mesh.nodes {
        if node.position().x < 1e-6 {
            // At left end
            fixed_nodes.push(node_id);
        }
    }

    // Find nodes at the right end (x = beam_length) for applied load
    let mut loaded_nodes = Vec::new();
    for (&node_id, node) in &mesh.nodes {
        if (node.position().x - beam_length).abs() < 1e-6 {
            // At right end
            loaded_nodes.push(node_id);
        }
    }

    // Apply fixed support at left end
    let fixed_support = DirichletPatterns::cantilever_beam(fixed_nodes.clone());
    for bc in fixed_support {
        boundary_conditions.add_condition(rtx_fea::boundary::BoundaryCondition::Dirichlet(bc));
    }

    // Apply downward force at right end
    let applied_force = -1000.0; // N (downward)
    let tip_load = NeumannPatterns::distributed_line_load(
        loaded_nodes.clone(),
        DofComponent::DisplacementY,
        applied_force,
    );
    boundary_conditions.add_condition(rtx_fea::boundary::BoundaryCondition::Neumann(tip_load));

    println!(
        "  - Fixed support applied to {} nodes at left end",
        fixed_nodes.len()
    );
    println!(
        "  - Distributed load of {:.0} N applied to {} nodes at right end",
        applied_force,
        loaded_nodes.len()
    );

    // Step 4: Configure analysis
    println!("\n4. Configuring analysis...");
    let mut config = AnalysisConfig::default();
    config.name = "Cantilever Beam Analysis".to_string();
    config.description = "Static analysis of a cantilever beam under tip loading".to_string();

    // Solver options
    config.solver_options.tolerance = 1e-8;
    config.solver_options.use_gpu = true;
    config.solver_options.max_iterations = 1000;

    // Assembly options
    config.assembly_options.use_gpu = true;
    config.assembly_options.tolerance = 1e-12;

    println!("  - Analysis configured with GPU acceleration");
    println!(
        "  - Solver tolerance: {:.0e}",
        config.solver_options.tolerance
    );

    // Step 5: Run the analysis
    println!("\n5. Running static linear analysis...");

    let (analysis_result, analysis_time) = BenchmarkUtils::time_function(|| {
        let mut analysis =
            StaticLinearAnalysis::new(mesh.clone(), materials.clone(), boundary_conditions, config);

        analysis.run()
    });

    let results = analysis_result?;

    println!(
        "  - Analysis completed in {}",
        StringUtils::format_duration(analysis_time)
    );
    println!(
        "  - Convergence: {}",
        if results.convergence.converged {
            "SUCCESS"
        } else {
            "FAILED"
        }
    );
    println!("  - Solver iterations: {}", results.convergence.iterations);

    // Step 6: Post-process results
    println!("\n6. Post-processing results...");

    let max_displacement = results.max_displacement();
    println!("  - Maximum displacement: {:.6} m", max_displacement);

    // Calculate theoretical tip deflection for comparison
    let moment_of_inertia = beam_height.powi(3) / 12.0; // For unit width
    let theoretical_deflection =
        (applied_force.abs() * beam_length.powi(3)) / (3.0 * elastic_modulus * moment_of_inertia);

    println!(
        "  - Theoretical tip deflection: {:.6} m",
        theoretical_deflection
    );

    let error_percentage =
        ((max_displacement - theoretical_deflection) / theoretical_deflection * 100.0).abs();
    println!("  - Error compared to theory: {:.2}%", error_percentage);

    // Find maximum stress if available
    if let Some(ref _stresses) = results.stresses {
        let max_stress = results.max_stress().unwrap_or(0.0);
        println!("  - Maximum stress: {:.2} MPa", max_stress / 1e6);
    }

    // Step 7: Output results
    println!("\n7. Writing output files...");

    // Prepare data for VTK output
    let nodes: Vec<(f64, f64, f64)> = mesh
        .nodes
        .values()
        .map(|node| {
            let pos = node.position();
            (pos.x, pos.y, pos.z)
        })
        .collect();

    let elements: Vec<Vec<usize>> = mesh
        .elements
        .values()
        .map(|element| element.nodes.iter().map(|id| id.0).collect())
        .collect();

    // Write VTK file for visualization
    VtkWriter::write_results(
        "cantilever_beam_results.vtk",
        &nodes,
        &elements,
        Some(&results.displacements),
        results.stresses.as_deref(),
    )?;

    println!("  - VTK file written: cantilever_beam_results.vtk");
    println!("  - Open with ParaView for visualization");

    // Step 8: Display timing information
    println!("\n8. Performance Summary:");
    println!("{}", results.timing);

    // Step 9: Validation
    println!("\n9. Validation:");
    if error_percentage < 5.0 {
        println!("  ✓ Results within 5% of theoretical solution");
    } else {
        println!("  ⚠ Results differ significantly from theoretical solution");
    }

    if results.convergence.converged {
        println!("  ✓ Analysis converged successfully");
    } else {
        println!("  ✗ Analysis failed to converge");
    }

    println!("\nAnalysis complete! 🎉");

    Ok(())
}

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

    #[test]
    fn test_cantilever_analysis() {
        // Simplified test version
        let mesh = Mesh::generate_rectangle(2.0, 0.2, 5, 3).unwrap();

        let mut materials = MaterialDatabase::new();
        let steel = LinearElastic::new(200e9, 0.3);
        materials.add_material(MaterialId(0), steel, Some("Steel".to_string()));

        let boundary_conditions = BoundaryConditionSet::new();
        let config = AnalysisConfig::default();

        let mut analysis = StaticLinearAnalysis::new(mesh, materials, boundary_conditions, config);

        // The analysis should be creatable without errors
        assert_eq!(analysis.analysis_type(), "Static Linear");
    }

    #[test]
    fn test_theoretical_deflection_calculation() {
        let force = 1000.0;
        let length = 10.0;
        let height = 1.0;
        let elastic_modulus = 200e9;

        let moment_of_inertia = height.powi(3) / 12.0;
        let deflection = (force * length.powi(3)) / (3.0 * elastic_modulus * moment_of_inertia);

        assert!(deflection > 0.0);
        assert!(deflection < 1.0); // Reasonable deflection
    }
}