From c514c1bdcc133bc923a17fa993dd52892723cdd2 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Mon, 15 Dec 2025 08:41:24 +0200 Subject: [PATCH] test(materials): add damage mechanics new API visionary test New 868-line test file for damage mechanics NEW API (test-driven development): - Visionary tests for isotropic damage model - Tests for ductile damage (Lemaitre model) - Tests for crack band regularization (mesh independence) - Tests for damage unloading and irreversibility - Extensive documentation of intended API design - Tests currently skipped (@test_skip) until implementation complete - Documents damage variable, effective stress, damage evolution - Documents regularization and coupled damage-plasticity This test file serves as both test suite and API design documentation for the new damage mechanics material model interface. --- test/materials/test_damage_new_api.jl | 868 ++++++++++++++++++++++++++ 1 file changed, 868 insertions(+) create mode 100644 test/materials/test_damage_new_api.jl diff --git a/test/materials/test_damage_new_api.jl b/test/materials/test_damage_new_api.jl new file mode 100644 index 0000000..09d022e --- /dev/null +++ b/test/materials/test_damage_new_api.jl @@ -0,0 +1,868 @@ +""" +# Damage Mechanics - NEW API (Test-Driven Development) + +**What:** Shows how damage models SHOULD work with the NEW API + +**Why:** +- **Stiffness degradation** - Material weakens (cracks, voids) +- **Irreversible** - Damage cannot heal (unlike plasticity unloading) +- **Mesh independence** - Regularization required (crack band, nonlocal) +- **Failure prediction** - Crack initiation, propagation + +**NEW API Concepts:** +1. **Damage variable** - d ∈ [0,1] where 0=intact, 1=failed +2. **Effective stress** - σ̄ = σ/(1-d) (undamaged configuration) +3. **Damage evolution** - ḋ = f(ε, ε_max, damage parameters) +4. **Regularization** - Length scale to avoid mesh sensitivity +5. **Coupled damage-plasticity** - Combined degradation + permanent deformation + +**Test Problems:** + +## Test 1: Isotropic Damage +- Stiffness reduction: E_eff = (1-d) E +- Damage driven by strain energy +- Validates crack initiation + +## Test 2: Ductile Damage (Lemaitre) +- Coupled damage-plasticity +- Damage from plastic dissipation +- Validates void growth + +## Test 3: Crack Band Regularization +- Mesh-independent energy dissipation +- Length scale h = element size +- Validates objectivity + +## Test 4: Damage Unloading +- Permanent stiffness loss +- No damage healing +- Validates irreversibility + +**Expected Behavior (when implemented):** +✅ Damage variable d ∈ [0,1] +✅ Stiffness degrades smoothly +✅ Mesh-independent fracture energy +✅ No healing on unload +✅ Crack localization captured +✅ Failure criterion satisfied + +**Status:** 🚧 VISIONARY TEST - Implementation in progress +""" + +using Test +using JuliaFEM +using Tensors +using LinearAlgebra +using Statistics + +@testset "Damage Mechanics - NEW API (TDD)" begin + + # ============================================================================= + # ISOTROPIC DAMAGE + # ============================================================================= + + @testset "Isotropic Damage - Strain-Based (Visionary)" begin + @test_skip begin # Skip until implemented + + # Material parameters + E = 200e3 # Young's modulus (MPa) + ν = 0.3 # Poisson's ratio + ε_d0 = 0.001 # Damage threshold strain + ε_f = 0.01 # Failure strain + + # NEW: Isotropic damage material + material = IsotropicDamage( + E=E, + ν=ν, + damage_threshold=ε_d0, + failure_strain=ε_f, + evolution_law=:exponential # or :linear, :power + ) + + # Strain history (uniaxial tension) + ε_max = 0.015 # Beyond failure + n_steps = 200 + ε_history = range(0, ε_max, length=n_steps) + + σ_history = [] + d_history = [] # Damage variable + + # Internal state + state = DamageState( + d=0.0, # Damage variable (0=intact, 1=failed) + ε_eq_max=0.0 # Maximum equivalent strain (history) + ) + + for ε in ε_history + # Strain tensor (uniaxial tension) + ε_tensor = SymmetricTensor{2,3}(( + ε, 0.0, 0.0, + 0.0, -ν * ε, 0.0, + 0.0, 0.0, -ν * ε + )) + + # Compute stress (with damage) + σ, state_new = compute_stress_damage(material, ε_tensor, state) + + push!(σ_history, σ[1, 1]) + push!(d_history, state_new.d) + + state = state_new + end + + # Validate elastic region (no damage) + elastic_indices = findall(ε_history .<= ε_d0) + for i in elastic_indices + @test isapprox(σ_history[i], E * ε_history[i], rtol=0.01) + @test d_history[i] == 0.0 + end + + # Validate damage growth + damage_indices = findall((ε_history .> ε_d0) .& (ε_history .< ε_f)) + for i in damage_indices + @test 0.0 < d_history[i] < 1.0 + + # Effective stiffness reduces + E_eff = E * (1 - d_history[i]) + @test E_eff < E + end + + # Validate failure + failure_indices = findall(ε_history .>= ε_f) + for i in failure_indices + @test d_history[i] >= 0.99 # Nearly complete damage + @test σ_history[i] < 0.1 * maximum(σ_history) # Stress vanishes + end + + end + end + + # ============================================================================= + # DAMAGE EVOLUTION LAWS + # ============================================================================= + + @testset "Damage Evolution Laws (Visionary)" begin + @test_skip begin + + E = 200e3 + ν = 0.3 + ε_d0 = 0.001 + ε_f = 0.01 + + # Test different evolution laws + laws = [:linear, :exponential, :power] + + for law in laws + material = IsotropicDamage( + E=E, + ν=ν, + damage_threshold=ε_d0, + failure_strain=ε_f, + evolution_law=law + ) + + state = DamageState(d=0.0, ε_eq_max=0.0) + + # Strain at 50% between threshold and failure + ε_mid = (ε_d0 + ε_f) / 2 + ε_tensor = SymmetricTensor{2,3}((ε_mid, 0.0, 0.0, 0.0, 0.0, 0.0)) + + σ, state_new = compute_stress_damage(material, ε_tensor, state) + + # Damage should be growing + @test 0.0 < state_new.d < 1.0 + + # Different laws give different d values + println("Law: $law, d = $(state_new.d)") + end + + end + end + + # ============================================================================= + # DUCTILE DAMAGE (LEMAITRE MODEL) + # ============================================================================= + + @testset "Ductile Damage - Coupled Plasticity (Visionary)" begin + @test_skip begin + + # Coupled damage-plasticity + material = DuctileDamage( + # Elastic properties + E=200e3, + ν=0.3, + + # Plasticity + yield_stress=250.0, + hardening=IsotropicHardening(H=2000.0), + + # Damage (Lemaitre) + damage_threshold=0.001, # Plastic strain threshold + S_crit=1.0, # Critical damage value + s_damage=1.5, # Triaxiality sensitivity + damage_exponent=2.0 + ) + + # Strain history (tension to failure) + ε_max = 0.05 + n_steps = 200 + ε_history = range(0, ε_max, length=n_steps) + + σ_history = [] + d_history = [] + ε_p_history = [] + + state = DuctileDamageState( + ε_p=zero(SymmetricTensor{2,3}), + ε_p_eq=0.0, + d=0.0, + α=zero(SymmetricTensor{2,3}) + ) + + for ε in ε_history + ε_tensor = SymmetricTensor{2,3}((ε, 0.0, 0.0, 0.0, 0.0, 0.0)) + σ, state = compute_stress_ductile_damage(material, ε_tensor, state) + + push!(σ_history, σ[1, 1]) + push!(d_history, state.d) + push!(ε_p_history, state.ε_p_eq) + end + + # Validate coupling: Damage grows with plastic strain + @test all(diff(d_history[ε_p_history.>material.damage_threshold]) .>= 0) + + # Validate softening: Peak stress followed by descent + σ_max_idx = argmax(σ_history) + @test σ_max_idx < length(σ_history) # Not at end + + # After peak, stress decreases (softening) + @test σ_history[end] < σ_history[σ_max_idx] + + # Damage increases monotonically + @test all(diff(d_history) .>= 0) + + end + end + + # ============================================================================= + # CRACK BAND REGULARIZATION + # ============================================================================= + + @testset "Crack Band Regularization (Visionary)" begin + @test_skip begin + + # Fracture energy per unit area (N/mm) + G_f = 0.1 # Fracture energy + + # Two different mesh sizes + h_coarse = 10.0 # mm + h_fine = 2.0 # mm + + # Crack band materials (adjust ε_f based on mesh) + # Energy = G_f = ∫ σ dε * h + # For linear softening: G_f ≈ (1/2) σ_max ε_f * h + + σ_max = 10.0 # Tensile strength (MPa) + + # Coarse mesh: larger ε_f + ε_f_coarse = 2 * G_f / (σ_max * h_coarse) + material_coarse = IsotropicDamage( + E=200e3, + ν=0.3, + damage_threshold=σ_max / 200e3, + failure_strain=ε_f_coarse, + evolution_law=:linear, + crack_band_width=h_coarse + ) + + # Fine mesh: smaller ε_f + ε_f_fine = 2 * G_f / (σ_max * h_fine) + material_fine = IsotropicDamage( + E=200e3, + ν=0.3, + damage_threshold=σ_max / 200e3, + failure_strain=ε_f_fine, + evolution_law=:linear, + crack_band_width=h_fine + ) + + # Compute energy dissipation for both + function compute_dissipation(material, ε_max, n_steps) + ε_history = range(0, ε_max, length=n_steps) + σ_history = [] + + state = DamageState(d=0.0, ε_eq_max=0.0) + + for ε in ε_history + ε_tensor = SymmetricTensor{2,3}((ε, 0.0, 0.0, 0.0, 0.0, 0.0)) + σ, state = compute_stress_damage(material, ε_tensor, state) + push!(σ_history, σ[1, 1]) + end + + # Integrate σ dε (trapezoid rule) + W = sum((σ_history[i] + σ_history[i+1]) / 2 * (ε_history[i+1] - ε_history[i]) + for i in 1:length(ε_history)-1) + + return W + end + + W_coarse = compute_dissipation(material_coarse, 1.2 * ε_f_coarse, 500) + W_fine = compute_dissipation(material_fine, 1.2 * ε_f_fine, 500) + + # Fracture energy per volume + G_v_coarse = W_coarse * h_coarse + G_v_fine = W_fine * h_fine + + # Should be mesh-independent! + @test isapprox(G_v_coarse, G_v_fine, rtol=0.1) + @test isapprox(G_v_coarse, G_f, rtol=0.1) + + end + end + + # ============================================================================= + # DAMAGE UNLOADING (IRREVERSIBILITY) + # ============================================================================= + + @testset "Damage Unloading - Irreversible (Visionary)" begin + @test_skip begin + + material = IsotropicDamage( + E=200e3, + ν=0.3, + damage_threshold=0.001, + failure_strain=0.01, + evolution_law=:linear + ) + + # Load-unload cycle + ε_max = 0.005 # Partial damage + + # Loading + ε_loading = range(0, ε_max, length=100) + σ_loading = [] + d_loading = [] + + state = DamageState(d=0.0, ε_eq_max=0.0) + + for ε in ε_loading + ε_tensor = SymmetricTensor{2,3}((ε, 0.0, 0.0, 0.0, 0.0, 0.0)) + σ, state = compute_stress_damage(material, ε_tensor, state) + + push!(σ_loading, σ[1, 1]) + push!(d_loading, state.d) + end + + d_max = state.d # Damage at peak load + + # Unloading + ε_unloading = range(ε_max, 0, length=100) + σ_unloading = [] + d_unloading = [] + + for ε in ε_unloading + ε_tensor = SymmetricTensor{2,3}((ε, 0.0, 0.0, 0.0, 0.0, 0.0)) + σ, state = compute_stress_damage(material, ε_tensor, state) + + push!(σ_unloading, σ[1, 1]) + push!(d_unloading, state.d) + end + + # Validate irreversibility + @test all(d_unloading .≈ d_max) # Damage does not heal! + + # Validate reduced stiffness + E_damaged = (1 - d_max) * material.E + + # Unloading slope should match damaged stiffness + # (linear regression on unloading curve) + ε_unload_vals = collect(ε_unloading) + slope = (σ_unloading[1] - σ_unloading[end]) / (ε_unload_vals[1] - ε_unload_vals[end]) + + @test isapprox(slope, E_damaged, rtol=0.1) + + end + end + + # ============================================================================= + # CONSISTENT TANGENT (DAMAGE) + # ============================================================================= + + @testset "Consistent Tangent - Damage (Visionary)" begin + @test_skip begin + + material = IsotropicDamage( + E=200e3, + ν=0.3, + damage_threshold=0.001, + failure_strain=0.01, + evolution_law=:exponential + ) + + # Strain state (damaged) + ε = SymmetricTensor{2,3}((0.003, 0.001, 0.0, 0.001, 0.002, 0.0)) + state = DamageState(d=0.3, ε_eq_max=0.003) + + # Compute stress and tangent + σ, state_new, C_damage = compute_stress_tangent_damage(material, ε, state) + + # Validate tangent via finite difference + δε = 1e-8 + + for i in 1:6 # Voigt notation + ε_pert = ε + δε * basis_symmetric_tensor(i) + σ_pert, _ = compute_stress_damage(material, ε_pert, state) + + dσ_numerical = (σ_pert - σ) / δε + dσ_tangent = C_damage ⊡ basis_symmetric_tensor(i) + + @test isapprox(dσ_numerical, dσ_tangent, rtol=0.01) + end + + # Validate symmetry + for i in 1:6, j in 1:6 + @test isapprox(C_damage[i, j], C_damage[j, i], atol=1e-10) + end + + # Validate degradation + C_elastic = compute_elastic_stiffness(material) + + # Damaged stiffness should be less + @test norm(C_damage) < norm(C_elastic) + + end + end + + # ============================================================================= + # PSEUDO-CODE: DAMAGE INTEGRATION + # ============================================================================= + + @testset "Damage Integration Pattern (Visionary)" begin + # Pseudo-code showing damage evolution + + println("\n" * "="^70) + println("DAMAGE INTEGRATION") + println("="^70) + + integration_pseudo = """ + # Damage evolution (strain-based isotropic) + + function compute_stress_damage(material, ε, state_old) + # 1. Compute equivalent strain + ε_eq = compute_equivalent_strain(ε) # e.g., sqrt(ε:ε) + + # 2. Update history (loading surface) + ε_eq_max = max(state_old.ε_eq_max, ε_eq) + + # 3. Check damage threshold + if ε_eq_max <= material.ε_d0 + # No damage + d = 0.0 + else + # Damage evolution + d = compute_damage(material, ε_eq_max) + end + + # 4. Effective stress + # Strain energy equivalence: W = W̄ + # σ : ε = σ̄ : ε in undamaged configuration + + # Elastic stress (undamaged) + C_elastic = compute_elastic_stiffness(material) + σ_undamaged = C_elastic ⊡ ε + + # Apply damage + σ = (1 - d) * σ_undamaged + + # 5. Update state + state_new = DamageState(d, ε_eq_max) + + return σ, state_new + end + + # Damage evolution laws + function compute_damage(material, ε_eq_max) + ε_d0 = material.damage_threshold + ε_f = material.failure_strain + + if material.evolution_law == :linear + # Linear: d = (ε - ε_d0) / (ε_f - ε_d0) + d = (ε_eq_max - ε_d0) / (ε_f - ε_d0) + + elseif material.evolution_law == :exponential + # Exponential: d = 1 - exp(-α(ε - ε_d0)) + α = -log(0.01) / (ε_f - ε_d0) # d(ε_f) ≈ 0.99 + d = 1 - exp(-α * (ε_eq_max - ε_d0)) + + elseif material.evolution_law == :power + # Power law: d = ((ε - ε_d0)/(ε_f - ε_d0))^n + n = 2.0 + d = ((ε_eq_max - ε_d0) / (ε_f - ε_d0))^n + end + + return clamp(d, 0.0, 0.99) # Numerical: never fully failed + end + """ + + println(integration_pseudo) + println("="^70) + println("✓ Equivalent strain: History variable") + println("✓ Loading surface: ε_eq_max = max(ε_eq_max_old, ε_eq)") + println("✓ Damage evolution: d = f(ε_eq_max)") + println("✓ Effective stress: σ = (1-d) σ_undamaged") + println("✓ Irreversible: d never decreases") + println("="^70) + end + + # ============================================================================= + # KEY ARCHITECTURAL INSIGHTS + # ============================================================================= + + println("\n" * "="^70) + println("DAMAGE MECHANICS ARCHITECTURE INSIGHTS (NEW API)") + println("="^70) + println("✓ Isotropic damage: Stiffness degradation (1-d)E") + println("✓ Damage variable: d ∈ [0,1] (0=intact, 1=failed)") + println("✓ Irreversible: d never decreases (no healing)") + println("✓ History: ε_eq_max (maximum strain ever reached)") + println("✓ Evolution laws: Linear, exponential, power") + println("✓ Crack band: Mesh-independent G_f via length scale h") + println("✓ Ductile damage: Coupled with plasticity (Lemaitre)") + println("✓ Consistent tangent: C_damage = ∂σ/∂ε (with damage)") + println("✓ Regularization: REQUIRED for mesh objectivity") + println("✓ Works with Newton-Krylov (tangent from damage law)") + println("="^70) + +end + +""" +# IMPLEMENTATION NOTES + +## Isotropic Damage + +### Damage Variable + +**Definition:** d ∈ [0,1] +- d = 0: Intact material +- d = 1: Completely damaged (failed) + +**Effective stress concept:** +σ̄ = σ / (1 - d) + +where σ̄ = stress in undamaged (effective) configuration. + +**Strain energy equivalence:** +W(σ, ε) = W̄(σ̄, ε) + +Implies: +σ = (1 - d) σ̄ + +where σ̄ = C_elastic : ε. + +### Equivalent Strain + +**For isotropic damage:** Need scalar measure of strain state. + +**Tension-driven:** +ε_eq = √(<ε_1>² + <ε_2>² + <ε_3>²) + +where <·> = positive part, ε_i = principal strains. + +**Reason:** Damage in tension (cracks open), not compression. + +**Alternative (modified von Mises):** +ε_eq = κ I_1 / (1-2ν) + √(3J_2) / (1+ν) + +where κ weighs volumetric vs deviatoric. + +### Damage Evolution Laws + +**Linear:** +d = (ε_eq - ε_d0) / (ε_f - ε_d0) for ε_eq ∈ [ε_d0, ε_f] + +**Exponential (smoother):** +d = 1 - exp(-α(ε_eq - ε_d0)) + +where α chosen such that d(ε_f) ≈ 0.99. + +**Power law:** +d = ((ε_eq - ε_d0) / (ε_f - ε_d0))^n + +where n controls softening rate. + +### History Variable + +**Loading surface:** ε_eq_max = max(ε_eq_history) + +**Damage depends on history:** +d = f(ε_eq_max) NOT f(ε_eq) + +**Irreversibility:** ε_eq_max only increases. + +**Update:** +```julia +ε_eq_max_new = max(ε_eq_max_old, ε_eq_current) +``` + +## Crack Band Regularization + +### Mesh Sensitivity Problem + +**Without regularization:** Fracture energy depends on mesh size! + +G_num = ∫ σ dε * h + +where h = element size. + +**Finer mesh → less energy dissipation → spurious brittleness.** + +### Crack Band Model + +**Idea:** Fracture happens over a band of width h. + +**Energy balance:** +G_f = ∫_0^{ε_f} σ dε * h + +where G_f = fracture energy per unit area (material property). + +**Adjust failure strain:** +ε_f = G_f / (∫_0^{ε_f} σ dε * h) + +**For linear softening:** +ε_f = 2 G_f / (σ_max h) + +where σ_max = tensile strength. + +**Result:** Mesh-independent fracture energy! + +### Implementation + +```julia +struct IsotropicDamage + E::Float64 + ν::Float64 + damage_threshold::Float64 + failure_strain::Float64 # Computed from G_f and h! + crack_band_width::Float64 # h (element size) +end + +function IsotropicDamage(; E, ν, G_f, σ_max, h) + ε_d0 = σ_max / E + ε_f = ε_d0 + 2*G_f / (σ_max * h) # Linear softening + + return IsotropicDamage(E, ν, ε_d0, ε_f, h) +end +``` + +## Ductile Damage (Lemaitre Model) + +### Coupling: Damage + Plasticity + +**Damage drives plasticity:** +σ_y_eff = σ_y / (1 - d) + +**Plasticity drives damage:** +ḋ = f(plastic dissipation) + +### Lemaitre Damage Evolution + +**Damage rate:** +ḋ = (Y / S)^s ε̇_p_eq + +where: +- Y = damage energy release rate = (σ_eq²) / (2E(1-d)²) +- S = material damage strength +- s = damage exponent +- ε̇_p_eq = equivalent plastic strain rate + +**Damage threshold:** +d = 0 until ε_p_eq > ε_p_threshold + +**Triaxiality influence:** +Y = Y(σ_eq, σ_m / σ_eq) + +where σ_m = mean stress (pressure). + +**High triaxiality → void growth → more damage.** + +### Integration + +```julia +function compute_ductile_damage(material, σ, ε_p_eq, state) + if ε_p_eq < material.ε_p_threshold + return 0.0 + end + + # Damage energy release rate + σ_eq = von_mises_stress(σ) + Y = σ_eq^2 / (2 * material.E * (1 - state.d)^2) + + # Triaxiality (optional) + σ_m = trace(σ) / 3 + η = σ_m / σ_eq + + # Damage increment + Δε_p = ε_p_eq - state.ε_p_eq_old + Δd = (Y / material.S)^material.s * Δε_p + + d_new = state.d + Δd + + return clamp(d_new, 0.0, 0.99) +end +``` + +## Consistent Tangent (Damage) + +**For Newton:** Need C_damage = dσ/dε. + +**Elastic damage:** +σ = (1 - d) C_elastic : ε + +**Tangent:** +C_damage = (1 - d) C_elastic + ∂d/∂ε ⊗ σ_elastic + +where ⊗ = outer product. + +**Derivative of damage:** +∂d/∂ε = (∂d/∂ε_eq) (∂ε_eq/∂ε) + +**Chain rule through damage evolution law.** + +**For exponential:** +∂d/∂ε_eq = α exp(-α(ε_eq - ε_d0)) + +**For linear:** +∂d/∂ε_eq = 1 / (ε_f - ε_d0) + +### Symmetry + +**Major symmetry:** C_damage may NOT be symmetric if ∂d/∂ε ⊗ σ not symmetric. + +**Options:** +1. Symmetrize: C_sym = (C + C^T) / 2 +2. Use unsymmetric solver (GMRES handles it!) + +## Internal State Storage + +**Per integration point:** +```julia +struct DamageState + d::Float64 # Damage variable + ε_eq_max::Float64 # Maximum equivalent strain (history) +end + +# Coupled damage-plasticity +struct DuctileDamageState{dim} + ε_p::SymmetricTensor{2,dim} + ε_p_eq::Float64 + d::Float64 + α::SymmetricTensor{2,dim} +end +``` + +**Element-level:** +```julia +struct DamageElement + topology::AbstractTopology + basis::AbstractBasis + nodes::NTuple{N,Int} + state::Vector{DamageState} # Per integration point +end +``` + +## Nodal Assembly (Damage) + +```julia +function tangent_matvec_damage!(w, v, u_current, material, elements, states) + Threads.@threads for node_i in 1:n_nodes + w_local = zero(Vec{3}) + + for elem in node_to_elements[node_i] + for (ip_idx, ip) in enumerate(integration_points(elem)) + # Current state + state = states[elem][ip_idx] + + # Strain + ε = compute_strain(elem, ip, u_current) + + # Consistent tangent (with damage) + σ, state_new, C_damage = compute_stress_tangent_damage(material, ε, state) + + for node_j in elem.nodes + # Tangent block (damaged stiffness) + K_t_ij = compute_damage_tangent_block(elem, node_i, node_j, C_damage, ip) + + v_j = Vec{3}(v[3*(node_j-1)+1:3*node_j]) + w_local += K_t_ij ⊡ v_j + end + + # Update state + states[elem][ip_idx] = state_new + end + end + + w[3*(node_i-1)+1:3*node_i] = w_local + end +end +``` + +**Key:** Damage state updated each iteration! + +## Regularization Techniques + +### 1. Crack Band (Local) + +**Pros:** Simple, fast +**Cons:** Still some mesh sensitivity + +### 2. Nonlocal Damage + +**Averaged equivalent strain:** +ε̄_eq(x) = (1/V_R) ∫_{B_R(x)} α(||y-x||) ε_eq(y) dy + +where: +- B_R(x) = ball of radius R around x +- α = weight function (Gaussian) + +**Damage driven by ε̄_eq instead of ε_eq.** + +**Pros:** Mesh-independent +**Cons:** Expensive (nonlocal averaging) + +### 3. Gradient Damage + +**Higher-order PDE:** +ε̄_eq - c ∇²ε̄_eq = ε_eq + +where c = internal length scale. + +**Requires additional DOF or coupled system.** + +**Pros:** Mesh-independent, smooth localization +**Cons:** Complex implementation + +### 4. Phase Field (Future) + +**Crack as diffuse interface:** +φ(x) ∈ [0,1] where φ=1 is crack. + +**Coupled:** +- Elasticity with φ-dependent stiffness +- Allen-Cahn or Ginzburg-Landau equation for φ + +**Pros:** Arbitrary crack topology, no remeshing +**Cons:** Very expensive + +## Next Steps + +1. Implement `IsotropicDamage` material type +2. Implement `DamageState` struct +3. Implement damage evolution laws +4. Implement crack band regularization +5. Implement `DuctileDamage` (coupled) +6. Implement consistent tangent +7. Validate mesh independence +8. Validate against experiments +9. Performance benchmarks + +"""