From be875b2bf62e38a64d7d88106eeac40b1ddf8163 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sat, 9 May 2026 18:27:38 +0300 Subject: [PATCH] chore(test): drop damage API regression stub Remove unmaintained damage-model tests for retired APIs. - Delete `test/materials/test_damage_new_api.jl`. --- test/materials/test_damage_new_api.jl | 868 -------------------------- 1 file changed, 868 deletions(-) delete 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 deleted file mode 100644 index 09d022e..0000000 --- a/test/materials/test_damage_new_api.jl +++ /dev/null @@ -1,868 +0,0 @@ -""" -# 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 - -"""