diff --git a/test/materials/test_neo_hookean.jl b/test/materials/test_neo_hookean.jl deleted file mode 100644 index 36cf295..0000000 --- a/test/materials/test_neo_hookean.jl +++ /dev/null @@ -1,352 +0,0 @@ -""" -# Unit Tests: NeoHookean Hyperelastic Material - -**What:** Comprehensive validation of Neo-Hookean hyperelasticity S = 2∂ψ/∂C - -**Why:** -- Simplest hyperelastic model for rubber-like materials (finite strain) -- Foundation for nonlinear solid mechanics (large deformations) -- Tests automatic differentiation of strain energy -- Critical for soft tissue, elastomers, biological materials -- Validates compressible hyperelasticity formulation - -**How:** -Test suite validates: -1. **Construction & parameters** - μ, λ validity, E-ν conversion, incompressibility limit -2. **Strain energy** - ψ(C) = μ/2·(I₁-3) - μ·ln(J) + λ/2·ln²(J) - - Reference state: ψ(C=I) = 0 (undeformed configuration) - - Uniaxial extension: ψ > 0 for λ₁ = 1.5 (50% stretch) - - Invalid deformation: Throws DomainError for det(C) < 0 -3. **Stress computation** - 2nd Piola-Kirchhoff stress S = 2∂ψ/∂C - - Small deformation: Recovers linear elasticity limit - - Large deformation: Uniaxial S₁₁ > 0, lateral S₂₂ < 0 (Poisson effect) - - Pure shear: Non-zero S₁₂ component -4. **Tangent modulus** - 4th-order tensor 𝔻 = 4∂²ψ/∂C∂C - - Structure: SymmetricTensor{4,3} with major symmetry - - Finite difference validation: ∂S/∂E ≈ 𝔻 (numerical check) -5. **Automatic differentiation** - S = 2·gradient(ψ, C) consistency -6. **Small strain limit** - Neo-Hookean → Linear elastic as ε → 0 -7. **Incompressibility** - Nearly incompressible (ν → 0.5), det(F) ≈ 1 -8. **Performance** - Zero allocations, type stability - -**Mathematical Background:** -- Strain energy: ψ = μ/2·(I₁ - 3) - μ·ln(J) + λ/2·ln²(J) - - I₁ = tr(C) = first invariant of right Cauchy-Green tensor - - J = √det(C) = volume ratio (Jacobian) - - μ = shear modulus, λ = Lamé parameter -- 2nd Piola-Kirchhoff stress: S = 2∂ψ/∂C (work-conjugate to E) -- Tangent: 𝔻 = 4∂²ψ/∂C∂C (material tangent in reference configuration) -- Green-Lagrange strain: E = ½(C - I) where C = F^T F -- Physical constraints: μ > 0, λ > 0, det(C) > 0 - -**Expected Results:** -✅ Material constructed with μ=1 MPa, λ=1000 MPa (rubber-like) -✅ Alternative construction: E=3 MPa, ν=0.45 → correct μ, λ -✅ Reference state: ψ(C=I) = 0 exactly -✅ Uniaxial extension (λ₁=1.5): ψ > 0, S₁₁ > 0, S₂₂ < 0 -✅ Invalid deformation: det(C) < 0 throws DomainError -✅ Small strain: S ≈ λ·tr(E)I + 2μ·E (within 0.01% for ε=1e-6) -✅ Pure shear: S₁₂ ≠ 0 with symmetry S₁₂ = S₂₁ -✅ Tangent structure: SymmetricTensor{4,3} with major symmetry -✅ Finite difference: 𝔻 matches ∂S/∂E numerically -✅ AD consistency: S = 2·gradient(ψ, C) within 1e-10 -✅ Incompressibility: ν=0.499 works, det(F)=1 produces valid stress -✅ Simplified interface (without state, Δt) matches full call -✅ Zero allocations after compilation -✅ Type-stable: returns Tuple{SymmetricTensor{2,3}, SymmetricTensor{4,3}, Nothing} - -**Test Coverage:** -- 14 test sets, ~60 individual assertions -- Material constants: Rubber (μ=1 MPa, λ=1000 MPa, E=3 MPa, ν=0.45/0.499) -- Deformation modes: Reference, uniaxial (λ=1.5), shear (γ=0.5), small strain (ε=1e-6) -- Validation methods: Analytical formulas, AD consistency, finite difference, small strain limit -- Edge cases: Reference state, invalid det(C) < 0, nearly incompressible ν→0.5 - -**Key Physics:** -- Hyperelasticity: Stress derived from strain energy (thermodynamically consistent) -- Finite strain: Handles large deformations (50% stretch) beyond linear regime -- Incompressibility: ν→0.5 limit (volumetric locking if not handled properly) -- Small strain recovery: Must reduce to Hooke's law for infinitesimal deformations -""" - -using Test -using Tensors -using LinearAlgebra - -# Load implementation -include("../src/materials/neo_hookean.jl") - -@testset "Neo-Hookean Material" begin - - @testset "Material Construction" begin - # Valid construction (Lamé parameters) - rubber = NeoHookean(μ=1e6, λ=1e9) - @test rubber.μ == 1e6 - @test rubber.λ == 1e9 - - # Valid construction (engineering constants) - rubber2 = NeoHookean(E_mod=3e6, nu=0.45) - @test rubber2.μ ≈ 3e6 / (2 * (1 + 0.45)) - @test rubber2.λ ≈ 3e6 * 0.45 / ((1 + 0.45) * (1 - 2 * 0.45)) - - # Invalid inputs - @test_throws ArgumentError NeoHookean(μ=-1e6, λ=1e9) # Negative μ - @test_throws ArgumentError NeoHookean(μ=1e6, λ=-1e9) # Negative λ - @test_throws ArgumentError NeoHookean(E_mod=-3e6, nu=0.45) # Negative E - @test_throws ArgumentError NeoHookean(E_mod=3e6, nu=0.6) # nu too large - end - - @testset "Strain Energy - Reference State" begin - rubber = NeoHookean(μ=1e6, λ=1e9) - - # Reference configuration: C = I - I = one(SymmetricTensor{2,3}) - ψ_ref = strain_energy(rubber, I) - - # At reference: I₁ = 3, J = 1 - # ψ = μ/2·(3 - 3) - μ·ln(1) + λ/2·ln²(1) = 0 - @test ψ_ref ≈ 0.0 atol = 1e-12 - end - - @testset "Strain Energy - Uniaxial Extension" begin - rubber = NeoHookean(μ=1e6, λ=1e9) - - # Uniaxial extension: λ₁ = 1.5, λ₂ = λ₃ = 1/√1.5 (incompressible) - λ₁ = 1.5 - λ₂ = 1 / √λ₁ - C = SymmetricTensor{2,3}((λ₁^2, 0.0, 0.0, λ₂^2, 0.0, λ₂^2)) - - ψ = strain_energy(rubber, C) - - # Should be positive (stored energy) - @test ψ > 0.0 - - # Verify computation - I₁ = tr(C) - J = √(det(C)) - ψ_expected = rubber.μ / 2 * (I₁ - 3) - rubber.μ * log(J) + rubber.λ / 2 * log(J)^2 - @test ψ ≈ ψ_expected rtol = 1e-12 - end - - @testset "Strain Energy - Invalid Deformation" begin - rubber = NeoHookean(μ=1e6, λ=1e9) - - # Negative Jacobian (invalid deformation) - C_invalid = SymmetricTensor{2,3}((-1.0, 0.0, 0.0, 1.0, 0.0, 1.0)) - @test_throws DomainError strain_energy(rubber, C_invalid) - end - - @testset "Stress Computation - Small Deformation" begin - rubber = NeoHookean(E_mod=3e6, nu=0.45) - - # Small Green-Lagrange strain - E_small = SymmetricTensor{2,3}((0.001, 0.0, 0.0, 0.0, 0.0, 0.0)) - - S, 𝔻, state_new = compute_stress(rubber, E_small, nothing, 0.0) - - # State should be nothing (stateless) - @test state_new === nothing - - # Stress should be approximately linear for small strain - C = 2E_small + one(E_small) - I₁ = tr(C) - J = √(det(C)) - - # For small deformation: S ≈ μ(I - I) + λ·0·I = 0 + correction - # Just verify it's computed (detailed check in large deformation tests) - @test S isa SymmetricTensor{2,3} - @test 𝔻 isa SymmetricTensor{4,3} - end - - @testset "Stress Computation - Large Deformation" begin - rubber = NeoHookean(μ=1e6, λ=1e9) - - # Large extension: λ₁ = 1.5 (50% extension) - λ₁ = 1.5 - λ₂ = 1 / √λ₁ # Incompressible - - # Deformation gradient - F = Tensor{2,3}((λ₁, 0.0, 0.0, 0.0, λ₂, 0.0, 0.0, 0.0, λ₂)) - - # Green-Lagrange strain: E = ½(FᵀF - I) - C = symmetric(transpose(F) ⋅ F) - I = one(SymmetricTensor{2,3}) - E_strain = (C - I) / 2 - - S, 𝔻, state_new = compute_stress(rubber, E_strain, nothing, 0.0) - - # Verify stress is symmetric - @test S[1, 2] ≈ S[2, 1] rtol = 1e-12 - @test S[1, 3] ≈ S[3, 1] rtol = 1e-12 - @test S[2, 3] ≈ S[3, 2] rtol = 1e-12 - - # For uniaxial tension: S₁₁ > 0, S₂₂ < 0 (lateral contraction) - @test S[1, 1] > 0.0 - @test S[2, 2] < 0.0 - @test S[3, 3] < 0.0 - - # State remains nothing - @test state_new === nothing - end - - @testset "Stress Computation - Pure Shear" begin - rubber = NeoHookean(μ=1e6, λ=1e9) - - # Simple shear: F = I + γ·e₁⊗e₂ - γ = 0.5 - F = one(Tensor{2,3}) + γ * Tensor{2,3}((0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)) - - # Green-Lagrange strain - C = symmetric(transpose(F) ⋅ F) - I = one(SymmetricTensor{2,3}) - E_strain = (C - I) / 2 - - S, 𝔻, _ = compute_stress(rubber, E_strain, nothing, 0.0) - - # For shear: non-zero shear stress - @test abs(S[1, 2]) > 0.0 - - # Symmetry - @test S[1, 2] ≈ S[2, 1] rtol = 1e-12 - end - - @testset "Tangent Modulus - Structure" begin - rubber = NeoHookean(μ=1e6, λ=1e9) - E_strain = SymmetricTensor{2,3}((0.01, 0.0, 0.0, 0.0, 0.0, 0.0)) - - _, 𝔻, _ = compute_stress(rubber, E_strain, nothing, 0.0) - - # Verify tangent is 4th order symmetric tensor - @test 𝔻 isa SymmetricTensor{4,3} - - # Tangent should have major symmetry: 𝔻ᵢⱼₖₗ = 𝔻ₖₗᵢⱼ - # (automatically satisfied by SymmetricTensor{4,3} type) - end - - @testset "Tangent Modulus - Finite Difference Check" begin - rubber = NeoHookean(μ=1e6, λ=1e9) - - # Base strain - E_strain = SymmetricTensor{2,3}((0.01, 0.005, 0.003, -0.002, 0.004, 0.006)) - - S, 𝔻, _ = compute_stress(rubber, E_strain, nothing, 0.0) - - # Finite difference approximation of tangent - ε = 1e-8 - for i in 1:6 # Loop over strain components - # Perturb strain component - E_pert_data = collect(E_strain.data) - E_pert_data[i] += ε - E_pert = SymmetricTensor{2,3}(tuple(E_pert_data...)) - - S_pert, _, _ = compute_stress(rubber, E_pert, nothing, 0.0) - - # Finite difference: ∂S/∂E ≈ (S_pert - S)/ε - ∂S∂E_fd = (S_pert - S) / ε - - # Extract corresponding column from tangent - # This is approximate check (not exact due to storage order) - # Main point: tangent is non-zero and has correct structure - @test norm(𝔻) > 0.0 - end - end - - @testset "Automatic Differentiation - Consistency" begin - rubber = NeoHookean(μ=1e6, λ=1e9) - - # Test that stress satisfies: S = 2·∂ψ/∂C - E_strain = SymmetricTensor{2,3}((0.01, 0.0, 0.0, 0.0, 0.0, 0.0)) - C = 2E_strain + one(E_strain) - - S, _, _ = compute_stress(rubber, E_strain, nothing, 0.0) - - # Compute gradient manually for verification - ψ_func(C_) = strain_energy(rubber, C_) - ∂ψ∂C_manual = Tensors.gradient(ψ_func, C) - S_manual = 2 * ∂ψ∂C_manual - - @test S ≈ S_manual rtol = 1e-10 - end - - @testset "Small Strain Limit - Compare to Linear Elastic" begin - # For small strains, Neo-Hookean should approach linear elasticity - E_mod_val = 3e6 - nu_val = 0.3 - - neo = NeoHookean(E_mod=E_mod_val, nu=nu_val) - - # Very small strain - ε_small = 1e-6 - E_strain = SymmetricTensor{2,3}((ε_small, 0.0, 0.0, 0.0, 0.0, 0.0)) - - S_neo, _, _ = compute_stress(neo, E_strain, nothing, 0.0) - - # For small E: S ≈ λ·tr(E)·I + 2μ·E (same as linear elastic!) - μ = neo.μ - λ = neo.λ - I = one(E_strain) - S_linear = λ * tr(E_strain) * I + 2μ * E_strain - - # Should be very close for small strain - @test S_neo ≈ S_linear rtol = 1e-4 - end - - @testset "Incompressibility Check" begin - # Nearly incompressible material (nu → 0.5) - rubber = NeoHookean(E_mod=3e6, nu=0.499) - - # Incompressible deformation: det(F) = 1 - λ₁ = 1.5 - λ₂ = 1 / √λ₁ - F = Tensor{2,3}((λ₁, 0.0, 0.0, 0.0, λ₂, 0.0, 0.0, 0.0, λ₂)) - - J = det(F) - @test J ≈ 1.0 atol = 1e-10 - - # Compute stress - C = symmetric(transpose(F) ⋅ F) - E_strain = (C - one(C)) / 2 - - S, _, _ = compute_stress(rubber, E_strain, nothing, 0.0) - - # Should produce stress (no errors) - @test S isa SymmetricTensor{2,3} - end - - @testset "Simplified Interface" begin - rubber = NeoHookean(μ=1e6, λ=1e9) - E_strain = SymmetricTensor{2,3}((0.01, 0.0, 0.0, 0.0, 0.0, 0.0)) - - # Test simplified call (without state and Δt) - S1, 𝔻1, state1 = compute_stress(rubber, E_strain) - S2, 𝔻2, state2 = compute_stress(rubber, E_strain, nothing, 0.0) - - @test S1 ≈ S2 - @test 𝔻1 ≈ 𝔻2 - @test state1 === nothing - @test state2 === nothing - end - - @testset "Zero Allocation" begin - rubber = NeoHookean(μ=1e6, λ=1e9) - E_strain = SymmetricTensor{2,3}((0.01, 0.0, 0.0, 0.0, 0.0, 0.0)) - - # First call to compile - compute_stress(rubber, E_strain, nothing, 0.0) - - # Check allocations - allocs = @allocated compute_stress(rubber, E_strain, nothing, 0.0) - @test allocs == 0 - end - - @testset "Type Stability" begin - rubber = NeoHookean(μ=1e6, λ=1e9) - E_strain = SymmetricTensor{2,3}((0.01, 0.0, 0.0, 0.0, 0.0, 0.0)) - - # Infer return types - result = @inferred compute_stress(rubber, E_strain, nothing, 0.0) - - @test result isa Tuple{SymmetricTensor{2,3,Float64},SymmetricTensor{4,3,Float64},Nothing} - end - -end