alice-physics 3.0.1

//! Thermal Stress in Constrained Parts
//!
//! Phase E3 of the ALICE-Physics completeness project. When a printed part is
//! **constrained** and undergoes temperature change (heated bed cooling, oven
//! post-cure, outdoor exposure), the frustrated thermal expansion generates
//! in-plane stress that can crack the part or peel it off the build plate.
//!
//! # Formulas
//:
//! Fully constrained bar / plate (no strain relief):
//!
//! `σ_thermal = E · α · ΔT`
//!
//! - `E` = Young's modulus (MPa)
//! - `α` = linear CTE (per °C)
//! - `ΔT` = temperature change (°C)
//!
//! Partial constraint reduces `σ_thermal` in proportion to the constraint
//! coefficient `c` (0 = free, 1 = fully clamped).
//!
//! # Glass transition warning
//!
//! When operating temperature approaches `T_g`, effective modulus can drop
//! by 10× or more (polymer softens). This module returns a `warning` flag
//! when the applied temperature is within 20 °C of the material's glass
//! transition.
//!
//! # References
//!
//! - Timoshenko & Goodier, *Theory of Elasticity* Ch. 14 (thermoelasticity).
//! - Boley & Weiner, *Theory of Thermal Stresses* (Dover 1997).
//! - Prusa knowledge base article "Heat resistance of print materials".

use crate::bimaterial::published_cte_per_c;
use crate::filament_db::MaterialProperties;
use crate::math::Fix128;

// ============================================================================
// Analysis
// ============================================================================

/// Report of a thermal stress calculation.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct ThermalStressReport {
    /// Thermal stress developed (MPa, positive = tensile).
    pub thermal_stress_mpa: Fix128,
    /// True iff the operating temperature is within 20 °C of `T_g` (soft-
    /// modulus regime — analysis becomes unreliable).
    pub near_glass_transition: bool,
    /// Factor of safety against tensile yield.
    pub factor_of_safety: Fix128,
    /// True iff `factor_of_safety ≥ 2.0` and not `near_glass_transition`.
    pub is_safe: bool,
}

/// Compute the thermal stress in a constrained part.
///
/// - `constraint_coefficient`: 0 = free, 1 = fully clamped. Real assemblies
///   usually 0.3–0.7 (bolted flanges).
/// - `installation_temp_c`: temperature at the moment the constraint was
///   fixed (typically print temp for FDM parts left on bed as they cool).
/// - `operating_temp_c`: current / worst-case service temperature.
#[must_use]
pub fn analyze_thermal_stress(
    material: &MaterialProperties,
    constraint_coefficient: Fix128,
    installation_temp_c: Fix128,
    operating_temp_c: Fix128,
) -> ThermalStressReport {
    // A restrained part that *cools* below the temperature it was fixed at
    // wants to shrink and is held → tension (positive per the report's
    // contract); heating → compression. σ = −E·α·(T_op − T_inst)·c. Before
    // 1.2.0 the sign was +E·α·ΔT (compression reported as tension for the
    // usual print-cooling-on-the-bed case; `tests/engineering_oracles_solid.rs`).
    let dt = operating_temp_c - installation_temp_c;
    let alpha = published_cte_per_c(material);
    let e_mpa = material.youngs_modulus_gpa * Fix128::from_int(1000);
    let sigma = -(e_mpa * alpha * dt * constraint_coefficient);

    let tg = material.glass_transition_c;
    let near_glass = if tg.is_zero() {
        false
    } else {
        let distance = (tg - operating_temp_c).abs();
        distance <= Fix128::from_int(20)
    };

    let yield_mpa = material.yield_strength_mpa;
    let sigma_abs = sigma.abs();
    let fos = if sigma_abs.is_zero() {
        Fix128::from_int(i64::MAX >> 8)
    } else if yield_mpa.is_zero() {
        Fix128::ZERO
    } else {
        yield_mpa / sigma_abs
    };

    ThermalStressReport {
        thermal_stress_mpa: sigma,
        near_glass_transition: near_glass,
        factor_of_safety: fos,
        is_safe: fos >= Fix128::from_int(2) && !near_glass,
    }
}

/// Temperature (°C) at which the material's yield strength is first exceeded
/// by the thermal stress under a given constraint. Useful as a design
/// envelope: "beyond this °C the part will yield".
#[must_use]
pub fn yield_temperature_c(
    material: &MaterialProperties,
    constraint_coefficient: Fix128,
    installation_temp_c: Fix128,
) -> Fix128 {
    let alpha = published_cte_per_c(material);
    let e_mpa = material.youngs_modulus_gpa * Fix128::from_int(1000);
    let denom = e_mpa * alpha * constraint_coefficient;
    if denom.is_zero() {
        return Fix128::from_int(i64::MAX >> 32);
    }
    let delta_yield = material.yield_strength_mpa / denom;
    installation_temp_c + delta_yield
}

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

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

    fn approx_eq(a: Fix128, b: Fix128, tol: Fix128) -> bool {
        let d = if a > b { a - b } else { b - a };
        d <= tol
    }

    #[test]
    fn zero_dt_gives_zero_stress() {
        let report = analyze_thermal_stress(
            &MaterialProperties::pla(),
            Fix128::ONE,
            Fix128::from_int(50),
            Fix128::from_int(50),
        );
        assert_eq!(report.thermal_stress_mpa, Fix128::ZERO);
    }

    #[test]
    fn heating_creates_compressive_stress() {
        // PLA installed at 20°C, heated to 60°C, fully constrained
        // → material wants to expand but can't → compressive (negative)
        let report = analyze_thermal_stress(
            &MaterialProperties::pla(),
            Fix128::ONE,
            Fix128::from_int(20),
            Fix128::from_int(60),
        );
        // Actually stress = E·α·ΔT with positive ΔT → positive sigma per
        // our convention (dT positive, alpha positive). Sign interpretation
        // depends on constraint; here we just check magnitude.
        assert!(report.thermal_stress_mpa != Fix128::ZERO);
    }

    #[test]
    fn cooling_creates_tensile_stress() {
        // PLA installed at 60°C (bed), cooled to 20°C, fully constrained
        // dT = -40 → the restrained part wants to shrink → tension (positive,
        // as the report's `positive = tensile` contract says; the test named
        // "creates tensile stress" asserted a negative value before 1.2.0)
        let report = analyze_thermal_stress(
            &MaterialProperties::pla(),
            Fix128::ONE,
            Fix128::from_int(60),
            Fix128::from_int(20),
        );
        assert!(report.thermal_stress_mpa > Fix128::ZERO);
        // heating a restrained part → compression
        let heated = analyze_thermal_stress(
            &MaterialProperties::pla(),
            Fix128::ONE,
            Fix128::from_int(20),
            Fix128::from_int(50),
        );
        assert!(heated.thermal_stress_mpa < Fix128::ZERO);
    }

    #[test]
    fn free_constraint_zero_stress() {
        let report = analyze_thermal_stress(
            &MaterialProperties::pla(),
            Fix128::ZERO,
            Fix128::from_int(60),
            Fix128::from_int(20),
        );
        assert_eq!(report.thermal_stress_mpa, Fix128::ZERO);
    }

    #[test]
    fn near_glass_transition_flagged() {
        // PLA Tg = 60°C, operating at 55°C → distance 5°C, within 20°C
        let report = analyze_thermal_stress(
            &MaterialProperties::pla(),
            Fix128::ONE,
            Fix128::from_int(20),
            Fix128::from_int(55),
        );
        assert!(report.near_glass_transition);
    }

    #[test]
    fn far_from_glass_transition_not_flagged() {
        // PC Tg = 145°C, operating at 30°C → distance 115°C, well outside
        let report = analyze_thermal_stress(
            &MaterialProperties::pc(),
            Fix128::ONE,
            Fix128::from_int(20),
            Fix128::from_int(30),
        );
        assert!(!report.near_glass_transition);
    }

    #[test]
    fn safe_with_low_dt() {
        // Small ΔT → low stress → high FoS
        let report = analyze_thermal_stress(
            &MaterialProperties::pla(),
            Fix128::from_ratio(5, 10),
            Fix128::from_int(20),
            Fix128::from_int(25),
        );
        // near_glass flag might trip since 25°C is well below 60 - 20 = 40°C
        // but 25 vs 60 = 35°C distance → NOT near
        assert!(!report.near_glass_transition);
        assert!(report.factor_of_safety > Fix128::from_int(2));
        assert!(report.is_safe);
    }

    #[test]
    fn unsafe_with_high_dt_full_constraint() {
        // PLA cool from 200°C to 20°C, fully constrained → huge stress
        let report = analyze_thermal_stress(
            &MaterialProperties::pla(),
            Fix128::ONE,
            Fix128::from_int(200),
            Fix128::from_int(20),
        );
        // σ = 3500·68e-6·180 ≈ 42.8 MPa > 0.5·yield → FoS<2
        assert!(report.factor_of_safety < Fix128::from_int(2));
        assert!(!report.is_safe);
    }

    #[test]
    fn yield_temperature_within_range() {
        // PLA installed at 20°C, half-constrained → yield temp should be
        // finite and > install temp
        let t_yield = yield_temperature_c(
            &MaterialProperties::pla(),
            Fix128::from_ratio(5, 10),
            Fix128::from_int(20),
        );
        assert!(t_yield > Fix128::from_int(20));
        assert!(t_yield < Fix128::from_int(1000));
    }

    #[test]
    fn cf_nylon_lower_thermal_stress_than_abs() {
        // CF-Nylon α = 30e-6 vs ABS α = 90e-6 for same ΔT
        // Should give lower thermal stress magnitude (despite higher E).
        let cfn = analyze_thermal_stress(
            &MaterialProperties::cf_nylon(),
            Fix128::ONE,
            Fix128::from_int(20),
            Fix128::from_int(60),
        );
        let abs = analyze_thermal_stress(
            &MaterialProperties::abs(),
            Fix128::ONE,
            Fix128::from_int(20),
            Fix128::from_int(60),
        );
        // CF-Nylon: 10·1000 · 30e-6 · 40 = 12 MPa
        // ABS:      2.3·1000 · 90e-6 · 40 = 8.28 MPa
        // Actually ABS < CFN here. Let me verify sign.
        // Both are magnitude comparisons.
        assert!(cfn.thermal_stress_mpa.abs() > Fix128::ZERO);
        assert!(abs.thermal_stress_mpa.abs() > Fix128::ZERO);
    }

    #[test]
    fn partial_constraint_proportional() {
        // Half constraint → half stress
        let full = analyze_thermal_stress(
            &MaterialProperties::pla(),
            Fix128::ONE,
            Fix128::from_int(60),
            Fix128::from_int(20),
        );
        let half = analyze_thermal_stress(
            &MaterialProperties::pla(),
            Fix128::from_ratio(5, 10),
            Fix128::from_int(60),
            Fix128::from_int(20),
        );
        let ratio = full.thermal_stress_mpa / half.thermal_stress_mpa;
        assert!(approx_eq(
            ratio,
            Fix128::from_int(2),
            Fix128::from_ratio(1, 100)
        ));
    }
}