diff --git a/test/materials/test_hyperelasticity_new_api.jl b/test/materials/test_hyperelasticity_new_api.jl deleted file mode 100644 index 46141f6..0000000 --- a/test/materials/test_hyperelasticity_new_api.jl +++ /dev/null @@ -1,693 +0,0 @@ -""" -# Hyperelasticity - NEW API (Test-Driven Development) - -**What:** Shows how hyperelastic materials SHOULD work with the NEW API - -**Why:** -- **Finite strain** - Large deformations (rubber, soft tissue, biomechanics) -- **Energy-based** - Strain energy function Ψ(F) -- **Frame-invariant** - Material objectivity (rotation independence) -- **Multiple models** - Neo-Hookean, Mooney-Rivlin, Ogden, etc. -- **Incompressibility** - Nearly incompressible (ν ≈ 0.5) - -**NEW API Concepts:** -1. **Hyperelastic material types** - Neo-Hookean, Mooney-Rivlin, Ogden -2. **Strain energy function** - Ψ(F) and derivatives -3. **Push-forward stress** - σ = (1/J) P F^T (Cauchy from 1st Piola-Kirchhoff) -4. **Tangent moduli** - C_ijkl = ∂²Ψ/∂F_ij∂F_kl -5. **Incompressibility constraint** - det(F) ≈ 1 - -**Test Problems:** - -## Test 1: Neo-Hookean (Simplest Hyperelastic) -- Ψ = μ/2 (I₁ - 3) - μ ln(J) + λ/2 (ln J)² -- Validates stress computation from energy -- Tests incompressibility limit - -## Test 2: Mooney-Rivlin (Two-Parameter) -- Ψ = C₁(I₁ - 3) + C₂(I₂ - 3) + κ/2 (J - 1)² -- Better for rubber than Neo-Hookean -- Validates second invariant I₂ - -## Test 3: Ogden Model (Multi-Term) -- Ψ = Σᵢ μᵢ/αᵢ (λ₁^αᵢ + λ₂^αᵢ + λ₃^αᵢ - 3) -- Uses principal stretches λᵢ -- Most accurate for rubber - -## Test 4: Uniaxial Tension Test -- Compare to experimental data -- Validates material parameters -- Tests large strain (λ > 2) - -**Expected Behavior (when implemented):** -✅ Stress computed from ∂Ψ/∂F correctly -✅ Tangent moduli symmetric and positive-definite -✅ Incompressibility enforced (J ≈ 1) -✅ Frame-invariant (rotations don't change Ψ) -✅ Matches experimental stress-strain curves -✅ Works with Newton-Krylov solver - -**Status:** 🚧 VISIONARY TEST - Implementation in progress -""" - -using Test -using JuliaFEM -using Tensors -using LinearAlgebra -using Statistics - -@testset "Hyperelasticity - NEW API (TDD)" begin - - # ============================================================================= - # NEO-HOOKEAN MODEL (SIMPLEST) - # ============================================================================= - - @testset "Neo-Hookean Material (Visionary)" begin - @test_skip begin # Skip until implemented - - # Material parameters - μ = 1000.0 # Shear modulus - λ = 2000.0 # Lame parameter (nearly incompressible) - - # NEW: Hyperelastic material type - material = NeoHookean( - μ=μ, - λ=λ, - formulation=:compressible # or :incompressible - ) - - # Test deformation gradient - F = Tensor{2,3}(( - 1.2, 0.1, 0.0, - 0.0, 0.9, 0.0, - 0.0, 0.0, 1.0 - )) - - # Compute strain energy - Ψ = strain_energy(material, F) - - # Neo-Hookean energy: - # Ψ = μ/2 (I₁ - 3) - μ ln(J) + λ/2 (ln J)² - C = tdot(F) # Right Cauchy-Green: C = F^T F - I₁ = tr(C) # First invariant - J = det(F) # Volume ratio - - Ψ_analytical = μ / 2 * (I₁ - 3) - μ * log(J) + λ / 2 * log(J)^2 - - @test isapprox(Ψ, Ψ_analytical, rtol=1e-10) - - # First Piola-Kirchhoff stress: P = ∂Ψ/∂F - P = first_piola_kirchhoff_stress(material, F) - - # Analytical P for Neo-Hookean - F_inv = inv(F) - P_analytical = μ * (F - tdot(F_inv)) + λ * log(J) * tdot(F_inv) - - @test P ≈ P_analytical rtol = 1e-10 - - # Cauchy stress: σ = (1/J) P F^T - σ = cauchy_stress(material, F) - σ_from_P = (1 / J) * P ⊡ transpose(F) - - @test σ ≈ σ_from_P rtol = 1e-10 - - end - end - - # ============================================================================= - # INCOMPRESSIBILITY CONSTRAINT - # ============================================================================= - - @testset "Incompressibility (Nearly) (Visionary)" begin - @test_skip begin - - # Nearly incompressible (ν → 0.5) - E = 1000.0 - ν = 0.499 # Nearly incompressible - - μ = E / (2 * (1 + ν)) - λ = E * ν / ((1 + ν) * (1 - 2ν)) # Very large! - - material = NeoHookean( - μ=μ, - λ=λ, - formulation=:compressible - ) - - # Deformation (should preserve volume) - F = Tensor{2,3}(( - 1.5, 0.0, 0.0, - 0.0, 1 / sqrt(1.5), 0.0, - 0.0, 0.0, 1 / sqrt(1.5) - )) - - J = det(F) - - # For incompressible: J = 1 (volume preserving) - @test isapprox(J, 1.0, atol=1e-3) - - # Hydrostatic pressure enforces incompressibility - σ = cauchy_stress(material, F) - p = -tr(σ) / 3 # Hydrostatic pressure - - # For nearly incompressible, pressure should be large - @test abs(p) > 1000.0 # Significant pressure - - end - end - - # ============================================================================= - # MOONEY-RIVLIN MODEL - # ============================================================================= - - @testset "Mooney-Rivlin Material (Visionary)" begin - @test_skip begin - - # Material parameters (typical for rubber) - C₁ = 0.5 # MPa - C₂ = 0.1 # MPa - κ = 100.0 # Bulk modulus (incompressibility) - - # NEW: Mooney-Rivlin material - material = MooneyRivlin( - C1=C₁, - C2=C₂, - bulk_modulus=κ - ) - - # Test deformation - F = Tensor{2,3}(( - 1.5, 0.2, 0.0, - 0.1, 0.8, 0.0, - 0.0, 0.0, 1.0 - )) - - # Strain energy: Ψ = C₁(I₁ - 3) + C₂(I₂ - 3) + κ/2 (J - 1)² - C = tdot(F) - I₁ = tr(C) - I₂ = 0.5 * (tr(C)^2 - tr(C ⊡ C)) # Second invariant - J = det(F) - - Ψ = strain_energy(material, F) - Ψ_analytical = C₁ * (I₁ - 3) + C₂ * (I₂ - 3) + κ / 2 * (J - 1)^2 - - @test isapprox(Ψ, Ψ_analytical, rtol=1e-10) - - # Stress computation - P = first_piola_kirchhoff_stress(material, F) - σ = cauchy_stress(material, F) - - # Validate symmetry of Cauchy stress - @test isapprox(σ, transpose(σ), atol=1e-10) - - end - end - - # ============================================================================= - # OGDEN MODEL (PRINCIPAL STRETCHES) - # ============================================================================= - - @testset "Ogden Material (Visionary)" begin - @test_skip begin - - # Ogden parameters (multi-term) - μ_terms = [1.0, 0.5, 0.2] # Shear moduli - α_terms = [2.0, 3.0, -2.0] # Exponents - κ = 100.0 - - # NEW: Ogden material - material = Ogden( - mu=μ_terms, - alpha=α_terms, - bulk_modulus=κ - ) - - # Deformation - F = Tensor{2,3}(( - 2.0, 0.0, 0.0, - 0.0, 0.6, 0.0, - 0.0, 0.0, 0.8 - )) - - # Compute principal stretches - C = tdot(F) - eigenvalues_C = eigvals(C) - λ = sqrt.(eigenvalues_C) # Principal stretches - - # Strain energy: Ψ = Σᵢ μᵢ/αᵢ (λ₁^αᵢ + λ₂^αᵢ + λ₃^αᵢ - 3) - Ψ = strain_energy(material, F) - - J = det(F) - Ψ_analytical = sum( - μ_terms[i] / α_terms[i] * (sum(λ .^ α_terms[i]) - 3) - for i in 1:3 - ) + κ / 2 * (J - 1)^2 - - @test isapprox(Ψ, Ψ_analytical, rtol=1e-10) - - end - end - - # ============================================================================= - # UNIAXIAL TENSION TEST - # ============================================================================= - - @testset "Uniaxial Tension (Large Strain) (Visionary)" begin - @test_skip begin - - # Geometry: Unit cube under tension - mesh = generate_mesh( - geometry=UnitCube(), - element_type=Hex8, - n_elements=(4, 4, 4) - ) - - # Material (Neo-Hookean rubber) - material = NeoHookean( - μ=1.0, # MPa - λ=10.0, # Nearly incompressible - density=1000.0 - ) - - # Physics - elastic_physics = ContinuumPhysics{Displacement}( - material=material, - formulation=FullThreeD(), - finite_strain=true # CRITICAL! - ) - - domain = Domain( - name="RUBBER", - elements=mesh.elements, - physics=elastic_physics - ) - - # Boundary conditions - bc_fixed = DirichletBC( - nodes=mesh.node_sets["LEFT"], - dof=:displacement, - values=[0.0, 0.0, 0.0] - ) - - # Applied stretch (λ = 2.0, 100% strain!) - u_applied = 1.0 # Stretch from 1.0 to 2.0 - - bc_stretch = DirichletBC( - nodes=mesh.node_sets["RIGHT"], - dof=:displacement, - component=:x, - value=u_applied - ) - - # Nonlinear problem (finite strain) - problem = NonlinearProblem( - domains=[domain], - boundary_conditions=[bc_fixed, bc_stretch] - ) - - # Newton-Krylov solver - solver = NewtonKrylovSolver( - max_iterations=20, - convergence_tol=1e-6, - krylov_solver=GMRES(restart=30), - line_search=BacktrackingLineSearch() - ) - - solution = solve!(problem, solver) - - # Extract stress-stretch curve - λ = 1.0 + u_applied # Stretch ratio - - # Compute engineering stress - F_avg = compute_average_deformation_gradient(solution) - P_avg = first_piola_kirchhoff_stress(material, F_avg) - - # Engineering stress: σ_eng = P_11 (1st Piola-Kirchhoff in x) - σ_eng = P_avg[1, 1] - - # For Neo-Hookean uniaxial: - # σ_eng = μ(λ - 1/λ²) - σ_analytical = material.μ * (λ - 1 / λ^2) - - @test isapprox(σ_eng, σ_analytical, rtol=0.05) # 5% tolerance - - end - end - - # ============================================================================= - # FRAME INVARIANCE (OBJECTIVITY) - # ============================================================================= - - @testset "Frame Invariance (Visionary)" begin - @test_skip begin - - material = NeoHookean(μ=1000.0, λ=2000.0) - - # Deformation gradient - F = Tensor{2,3}(( - 1.2, 0.1, 0.0, - 0.0, 0.9, 0.0, - 0.0, 0.0, 1.0 - )) - - # Rotation tensor (90° about z-axis) - θ = π / 2 - Q = Tensor{2,3}(( - cos(θ), -sin(θ), 0.0, - sin(θ), cos(θ), 0.0, - 0.0, 0.0, 1.0 - )) - - # Rotated deformation: F' = Q F - F_rotated = Q ⊡ F - - # Strain energy should be INVARIANT - Ψ = strain_energy(material, F) - Ψ_rotated = strain_energy(material, F_rotated) - - @test isapprox(Ψ, Ψ_rotated, rtol=1e-10) - - # Cauchy stress should transform: σ' = Q σ Q^T - σ = cauchy_stress(material, F) - σ_rotated = cauchy_stress(material, F_rotated) - - σ_transformed = Q ⊡ σ ⊡ transpose(Q) - - @test σ_rotated ≈ σ_transformed rtol = 1e-10 - - end - end - - # ============================================================================= - # TANGENT MODULI (FOR NEWTON) - # ============================================================================= - - @testset "Tangent Moduli (Visionary)" begin - @test_skip begin - - material = NeoHookean(μ=1000.0, λ=2000.0) - - F = Tensor{2,3}(( - 1.2, 0.1, 0.0, - 0.0, 0.9, 0.0, - 0.0, 0.0, 1.0 - )) - - # Material tangent: C_ijkl = ∂²Ψ/∂F_ij∂F_kl - C = material_tangent(material, F) - - # Validate major symmetry: C_ijkl = C_klij - # (Minor symmetries don't hold for finite strain) - for i in 1:3, j in 1:3, k in 1:3, l in 1:3 - @test isapprox(C[i, j, k, l], C[k, l, i, j], atol=1e-10) - end - - # Validate positive-definiteness - # For small perturbation δF, δ²Ψ = C_ijkl δF_ij δF_kl > 0 - δF = 0.01 * rand(Tensor{2,3}) - δ²Ψ = dcontract(dcontract(C, δF), δF) - - @test δ²Ψ > 0 # Positive-definite - - end - end - - # ============================================================================= - # PSEUDO-CODE: HYPERELASTIC ASSEMBLY - # ============================================================================= - - @testset "Hyperelastic Assembly Pattern (Visionary)" begin - # Pseudo-code showing finite strain assembly - - println("\n" * "="^70) - println("HYPERELASTIC ASSEMBLY (FINITE STRAIN)") - println("="^70) - - assembly_pseudo = """ - # For hyperelastic materials, assembly uses current configuration - - function compute_residual_hyperelastic(u, material, elements) - r = zeros(length(u)) - - for elem in elements - # Current deformation gradient: F = I + ∇u - X = reference_coordinates(elem) # Undeformed - x = X + u[elem.nodes] # Deformed - - for ip in integration_points(elem) - # Jacobian in reference config - J₀, dN_dX = jacobian(elem, ip, X) - - # Deformation gradient: F = ∂x/∂X - F = compute_deformation_gradient(x, dN_dX) - - # 1st Piola-Kirchhoff stress: P = ∂Ψ/∂F - P = first_piola_kirchhoff_stress(material, F) - - # Residual: r = ∫_Ω₀ P : ∇_X(δu) dΩ₀ - # = ∫_Ω₀ P_iJ (dN_I/dX_J) dΩ₀ - for I in 1:n_nodes - for i in 1:3 - for J in 1:3 - r[3*(I-1)+i] += P[i,J] * dN_dX[I,J] * J₀ * ip.weight - end - end - end - end - end - - return r - end - - # Tangent stiffness (for Newton): - # K_t = ∫_Ω₀ (dN_I/dX_K) C_iJkL (dN_J/dX_L) dΩ₀ - # where C_iJkL = ∂²Ψ/∂F_iJ∂F_kL - """ - - println(assembly_pseudo) - println("="^70) - println("✓ Uses reference configuration Ω₀ (not current!)") - println("✓ Deformation gradient F = ∂x/∂X") - println("✓ 1st Piola-Kirchhoff stress P = ∂Ψ/∂F") - println("✓ Tangent moduli C = ∂²Ψ/∂F²") - println("✓ Works with nodal assembly (same pattern!)") - println("="^70) - end - - # ============================================================================= - # KEY ARCHITECTURAL INSIGHTS - # ============================================================================= - - println("\n" * "="^70) - println("HYPERELASTICITY ARCHITECTURE INSIGHTS (NEW API)") - println("="^70) - println("✓ Hyperelastic materials: NeoHookean, MooneyRivlin, Ogden") - println("✓ Strain energy function Ψ(F) is fundamental") - println("✓ Stress from energy: P = ∂Ψ/∂F, σ = (1/J) P F^T") - println("✓ Tangent from energy: C = ∂²Ψ/∂F²") - println("✓ Frame-invariant: Rotations don't change Ψ") - println("✓ Incompressibility: det(F) ≈ 1 for rubber") - println("✓ Works with Newton-Krylov (unsymmetric OK)") - println("✓ Assembly in reference config (not current!)") - println("✓ Large strains: λ > 2 (100%+ strain)") - println("="^70) - -end - -""" -# IMPLEMENTATION NOTES - -## Hyperelastic Material Models - -### Neo-Hookean (Simplest) - -**Strain energy:** -Ψ = μ/2 (I₁ - 3) - μ ln(J) + λ/2 (ln J)² - -where: -- I₁ = tr(C) = tr(F^T F) (first invariant) -- J = det(F) (volume ratio) -- μ, λ = Lame parameters - -**1st Piola-Kirchhoff stress:** -P = ∂Ψ/∂F = μ(F - F^{-T}) + λ ln(J) F^{-T} - -where F^{-T} = (F^{-1})^T. - -**Cauchy stress:** -σ = (1/J) P F^T = μ/J (B - I) + λ/J ln(J) I - -where B = F F^T (left Cauchy-Green). - -**Use case:** Rubber at moderate strains (< 100%). - -### Mooney-Rivlin (Two-Parameter) - -**Strain energy:** -Ψ = C₁(I₁ - 3) + C₂(I₂ - 3) + κ/2 (J - 1)² - -where: -- I₁ = tr(C) -- I₂ = 1/2 [(tr C)² - tr(C²)] -- J = det(F) -- C₁, C₂ = material parameters -- κ = bulk modulus - -**Better fit for rubber** than Neo-Hookean. - -**Relation to Neo-Hookean:** C₂ = 0 → Neo-Hookean. - -### Ogden Model (Multi-Term) - -**Strain energy:** -Ψ = Σᵢ μᵢ/αᵢ (λ₁^{αᵢ} + λ₂^{αᵢ} + λ₃^{αᵢ} - 3) + κ/2 (J - 1)² - -where: -- λ₁, λ₂, λ₃ = principal stretches (eigenvalues of F) -- μᵢ, αᵢ = material parameters (typically 3-6 terms) - -**Most accurate for rubber** (fits experimental data well). - -**Implementation:** Requires spectral decomposition of C. - -## Stress Measures - -### First Piola-Kirchhoff (P) - -**Definition:** P = ∂Ψ/∂F - -**Properties:** -- Non-symmetric -- Force per reference area -- Work-conjugate to F - -**Use:** Weak form in reference config. - -### Cauchy Stress (σ) - -**Definition:** σ = (1/J) P F^T - -**Properties:** -- Symmetric -- True stress (force per current area) -- What we measure - -**Use:** Post-processing, failure criteria. - -### Second Piola-Kirchhoff (S) - -**Definition:** S = F^{-1} P = J F^{-1} σ F^{-T} - -**Properties:** -- Symmetric -- Work-conjugate to E (Green-Lagrange strain) -- Energy-conjugate - -**Use:** Theoretical derivations. - -## Incompressibility - -**Constraint:** det(F) = J = 1 (volume preserving) - -**Nearly incompressible:** ν → 0.5, λ → ∞ - -**Enforcement:** - -1. **Penalty:** Add κ/2 (J - 1)² to Ψ (large κ) -2. **Lagrange multiplier:** Introduce pressure p -3. **Mixed formulation:** (u, p) unknowns - -**JuliaFEM approach:** Penalty for compressible materials, mixed for truly incompressible. - -## Frame Invariance (Objectivity) - -**Definition:** Strain energy invariant under rigid rotations. - -**Mathematical:** Ψ(Q F) = Ψ(F) for all rotations Q. - -**Why:** Material doesn't "know" about global rotations. - -**Implementation:** Use invariants (I₁, I₂, I₃) or principal stretches (λᵢ). - -**Validation:** -```julia -Q = rotation_matrix(θ) -F_rotated = Q * F -@test strain_energy(F_rotated) ≈ strain_energy(F) -``` - -## Tangent Moduli - -**Material tangent:** -C_{iJkL} = ∂²Ψ/∂F_{iJ}∂F_{kL} - -**Spatial tangent:** -c_{ijkl} = (1/J) F_{iI} F_{jJ} F_{kK} F_{lL} C_{IJKL} - -**Symmetries:** -- Major: C_{iJkL} = C_{kLiJ} (from Ψ) -- Minor: Generally NOT symmetric in hyperelasticity - -**Use in Newton:** -δP = C : δF - -## Assembly (Finite Strain) - -**Residual (weak form):** -r = ∫_{Ω₀} P : ∇_X(δu) dΩ₀ - f_ext - -**In components:** -r_I^i = ∫_{Ω₀} P_{iJ} ∂N_I/∂X_J dΩ₀ - f_I^i - -**Tangent stiffness:** -K_{IJ}^{ik} = ∫_{Ω₀} ∂N_I/∂X_K C_{iJkL} ∂N_J/∂X_L dΩ₀ - -**Key differences from small strain:** -1. Integrate over Ω₀ (reference), not Ω (current) -2. Use ∂N/∂X (reference gradients), not ∂N/∂x -3. P (1st Piola-Kirchhoff), not σ (Cauchy) - -## Nodal Assembly (Hyperelastic) - -```julia -function tangent_matvec_hyperelastic!(w, v, u_current, material, elements, node_to_elements) - Threads.@threads for node_i in 1:n_nodes - w_local = zero(Vec{3}) - - for elem in node_to_elements[node_i] - # Current deformation - F = deformation_gradient(elem, u_current) - - # Material tangent - C = material_tangent(material, F) - - for node_j in elem.nodes - # Tangent block: K_t_ij - K_t_ij = compute_hyperelastic_tangent_block(elem, node_i, node_j, C, F) - - v_j = Vec{3}(v[3*(node_j-1)+1:3*node_j]) - w_local += K_t_ij ⊡ v_j - end - end - - w[3*(node_i-1)+1:3*node_i] = w_local - end -end -``` - -**Same pattern as linear elasticity!** Only the tangent computation changes. - -## Next Steps - -1. Implement `NeoHookean` material type -2. Implement `strain_energy` function -3. Implement `first_piola_kirchhoff_stress` -4. Implement `cauchy_stress` (push-forward) -5. Implement `material_tangent` (C_ijkl) -6. Implement `MooneyRivlin` material -7. Implement `Ogden` material (spectral decomposition) -8. Validate against analytical solutions -9. Validate against experimental data -10. Performance benchmarks - -"""