diff --git a/test/dofs/test_thmecd_hexa_physics.jl b/test/dofs/test_thmecd_hexa_physics.jl
deleted file mode 100644
index cbcaf4e..0000000
--- a/test/dofs/test_thmecd_hexa_physics.jl
+++ /dev/null
@@ -1,746 +0,0 @@
-"""
-🌟 THE ULTIMATE FINALE: Thermo-Hydro-Mechanical-Electric-Chemical-Damage (THM-ECD)
-
-This implements **SIX-FIELD** coupled physics - literally EVERYTHING!
-
-Field Variables:
-- T: Temperature (Float64) at VERTICES - continuous H¹ field
-- u: Displacement (Vec{3}) at VERTICES - continuous H¹ vector field
-- p: Pore pressure (Float64) at CELLS - discontinuous L² field
-- φ: Electric potential (Float64) at EDGES - H(curl) field
-- c: Chemical concentration (Float64) at VERTICES - continuous H¹ field
-- d: Damage variable (Float64) at VERTICES - continuous H¹ field **[NEW!]**
-
-═══════════════════════════════════════════════════════════════════════
-COMPLETE PHYSICS FORMULATION - FULLY COUPLED THM-ECD SYSTEM
-═══════════════════════════════════════════════════════════════════════
-
-1️⃣  THERMAL (Heat with Damage-Enhanced Transport):
-   ρcₚ ∂T/∂t - ∇·(κ(d)∇T) = Q + α_T·T₀·E/(1-2ν) ∇·∂u/∂t + β_T·∂p/∂t + S·∇·J + Q_chem
-   
-   κ(d) = κ₀·(1 + β_κ·d) - Cracks INCREASE thermal conductivity
-   
-   Why? Cracks create preferential heat paths (convection in voids)
-
-2️⃣  MECHANICAL (Elasticity with Damage Degradation):
-   ρ ∂²u/∂t² - ∇·σ = f
-   
-   σ = (1-d)·C : ε(u) - α_T·(T-T₀)·I - α_p·p·I - e^T·E - α_c·c·I
-       \_____/
-        Damage reduces stiffness!
-   
-   Classic Kachanov damage: E_damaged = (1-d)·E₀
-   d = 0: intact material
-   d = 1: complete failure
-
-3️⃣  HYDRAULIC (Flow with Damage-Enhanced Permeability):
-   S_s ∂p/∂t + α_p ∂(∇·u)/∂t + β_T ∂T/∂t - ∇·(k(d)/μ_f ∇p) = q - ζ·∇·J + q_chem
-   
-   k(d) = k₀·exp(β_k·d) - Cracks EXPONENTIALLY increase permeability!
-   
-   Why exponential? Cubic law: k ~ w³ where w = crack width ~ d
-   
-   Examples:
-   - Hydraulic fracturing: d increases → k increases 1000×
-   - Rock damage: k(intact) = 10⁻²⁰ m², k(damaged) = 10⁻¹⁵ m²
-
-4️⃣  ELECTRIC (Charge Transport with Damage):
-   ∇·D = ρ_e
-   J = σ_e(d)·E + S·(-κ∇T) + z·F·D_m·∇c
-   
-   σ_e(d) = σ_e0·(1 - β_σ·d) - Cracks REDUCE electrical conductivity
-   
-   Why? Cracks are insulators (unless filled with electrolyte!)
-
-5️⃣  CHEMICAL (Species Transport with Enhanced Diffusion):
-   ∂c/∂t + ∇·J_c = R(c,T) + S_crack·d·∂d/∂t
-   
-   J_c = -D_eff(p,T,d)·∇c + u̇·c + μ_m·c·E - D_T·c·∇T
-   
-   D_eff(d) = D₀·(1 + β_D·d) - Cracks increase diffusivity
-   S_crack·d·∂d/∂t - Fresh crack surfaces provide reactive sites!
-   
-   Applications:
-   - Stress corrosion cracking: cracks expose fresh metal
-   - Concrete spalling: cracks accelerate chloride ingress
-   - Shale gas: fractures enable gas diffusion
-
-6️⃣  DAMAGE EVOLUTION (NEW FIELD!):
-   ∂d/∂t = f_damage(Y, d, T, c, p, φ)
-   
-   Y = ½ε:(C:ε) - Elastic energy density (damage driving force)
-   
-   **Damage evolution law (unified):**
-   
-   ∂d/∂t = <Y - Y₀> / η · g(d) · h_T(T) · h_c(c) · h_p(p) · h_φ(φ)
-            \______/   \___/   \____/   \____/   \____/   \____/
-              Rate     Growth   Thermal  Chemical Pressure Electric
-   
-   Where:
-   - Y₀: Damage threshold [J/m³]
-   - η: Viscosity parameter [J·s/m³]
-   - g(d) = (1-d)^m: Damage evolution function
-   - h_T(T) = exp(β_T^d · T): Thermal activation (creep)
-   - h_c(c) = (1 + γ_c · c): Chemically-assisted damage (SCC)
-   - h_p(p) = (1 + γ_p · p): Pressure-assisted damage
-   - h_φ(φ) = (1 + γ_φ · |∇φ|): Electric field damage
-   
-   **Physical mechanisms:**
-   
-   a) **Mechanical damage**: Y > Y₀ → cracks grow
-   
-   b) **Thermal damage**: High T → accelerated creep → damage
-      - Concrete: Thermal spalling at T > 400°C
-      - Metals: Creep damage at T > 0.4·T_melt
-   
-   c) **Stress corrosion cracking (SCC)**: c + stress → accelerated damage
-      - Chloride SCC in stainless steel
-      - Hydrogen embrittlement
-      - Environmentally assisted cracking
-   
-   d) **Pressure damage**: High p → pore pressure fracture
-      - Hydraulic fracturing
-      - Overpressured reservoirs
-   
-   e) **Electrochemical damage**: φ gradients → corrosion → damage
-      - Galvanic corrosion
-      - Cathodic disbondment
-      - Electromigration in conductors
-
-═══════════════════════════════════════════════════════════════════════
-COUPLING MATRIX (30 OFF-DIAGONAL BLOCKS! 15 BIDIRECTIONAL PAIRS!)
-═══════════════════════════════════════════════════════════════════════
-
-        │  T          u          p          φ          c          d
-   ─────┼─────────────────────────────────────────────────────────────────
-     T  │ K_TT      K_Tu       K_Tp       K_Tφ       K_Tc       K_Td
-        │ κ(d)∇∇    (α_T)      (β_T)      (S)        (H_rxn)    (β_κ·∇T)
-   ─────┼─────────────────────────────────────────────────────────────────
-     u  │ K_uT      K_uu       K_up       K_uφ       K_uc       K_ud
-        │ (α_T)    (1-d)C:ε:ε  (α_p)      (e_kij)    (α_c)      (-C:ε:ε)
-   ─────┼─────────────────────────────────────────────────────────────────
-     p  │ K_pT      K_pu       K_pp       K_pφ       K_pc       K_pd
-        │ (β_T)     (α_p)     k(d)∇∇      (ζ)        (ν_f)      (β_k·k·∇p)
-   ─────┼─────────────────────────────────────────────────────────────────
-     φ  │ K_φT      K_φu       K_φp       K_φφ       K_φc       K_φd
-        │ (S)       (e_kij)    (ζ)       σ(d)∇∇   (z·F·D_m)   (-β_σ·σ·∇φ)
-   ─────┼─────────────────────────────────────────────────────────────────
-     c  │ K_cT      K_cu       K_cp       K_cφ       K_cc       K_cd
-        │ (D_T)     (adv)      (D_eff)    (μ_m)     D(d)∇∇      (β_D·D·∇c)
-   ─────┼─────────────────────────────────────────────────────────────────
-     d  │ K_dT      K_du       K_dp       K_dφ       K_dc       K_dd
-        │ (β_T^d)   (∂Y/∂ε)    (γ_p)      (γ_φ)      (γ_c)      (viscous)
-
-   **NEW DAMAGE COUPLINGS (12 blocks!):**
-   
-   K_Td: Thermal conductivity change with damage
-   K_dT: Thermal activation of damage (creep)
-   
-   K_ud: Stiffness degradation (main damage effect!)
-   K_du: Elastic energy drives damage
-   
-   K_pd: Permeability change with damage (HUGE effect!)
-   K_dp: Pressure-assisted damage
-   
-   K_φd: Conductivity change with damage
-   K_dφ: Electric field damage
-   
-   K_cd: Diffusivity change with damage
-   K_dc: Chemical damage (SCC, hydrogen embrittlement)
-   
-   K_dd: Rate-dependent damage evolution (viscoplasticity)
-
-═══════════════════════════════════════════════════════════════════════
-MATERIAL PARAMETERS (THM-ECD):
-═══════════════════════════════════════════════════════════════════════
-
-**Original THM-EC parameters:**
-(Same as before - see test_thmec_penta_physics.jl)
-
-**NEW Damage parameters:**
-
-Damage evolution:
-- Y₀ = 1e6 [J/m³]: Damage threshold
-- η = 1e12 [J·s/m³]: Viscosity parameter
-- m = 2.0: Damage evolution exponent
-
-Property degradation coefficients:
-- β_κ = 2.0: Thermal conductivity increase (cracks → convection)
-- β_k = 10.0: Permeability increase (exponential!)
-- β_σ = 0.8: Electrical conductivity decrease
-- β_D = 3.0: Diffusivity increase
-
-Coupled damage coefficients:
-- β_T^d = 0.001 [1/K]: Thermal damage activation
-- γ_c = 1e-3 [m³/mol]: Chemical damage enhancement (SCC)
-- γ_p = 1e-9 [1/Pa]: Pressure damage enhancement
-- γ_φ = 1e-8 [1/(V/m)]: Electric field damage
-
-═══════════════════════════════════════════════════════════════════════
-REAL-WORLD APPLICATIONS (Where ALL 6 Fields Matter):
-═══════════════════════════════════════════════════════════════════════
-
-1. **Geothermal Reservoir Stimulation**
-   - Inject cold water (T↓) → thermal stress → damage (d↑)
-   - Damage → permeability increase (k↑) → better flow (p)
-   - Mineral dissolution (c) at crack surfaces
-   - Electrokinetic effects (φ) from fluid flow
-   - Result: Enhanced geothermal system (EGS)
-
-2. **Nuclear Waste Canister Corrosion**
-   - Heat from decay (T) → thermal expansion (u)
-   - Groundwater pressure (p) → stress
-   - Corrosion reactions (c) → volume expansion → stress
-   - Galvanic currents (φ) → accelerated corrosion
-   - Stress + corrosion → damage (d) → canister failure
-   - Damage → permeability → radionuclide release
-
-3. **Hydraulic Fracturing (Fracking)**
-   - High pressure injection (p) → crack opening (d↑)
-   - Damage → permeability increase (k = k₀·exp(10·d)) → gas flow
-   - Thermal effects from deep formations (T)
-   - Chemical reactions with formation water (c)
-   - Electrokinetic effects from shale (φ)
-   - Result: Economic gas production
-
-4. **Reinforced Concrete Corrosion**
-   - Chloride ingress (c) through cracks (d)
-   - Rebar corrosion (c + φ → Fe²⁺) → expansion
-   - Expansion → cracking (d↑) → more chloride (c↑)
-   - Thermal cycles (T) → additional cracking
-   - Saturated pores (p) → freeze-thaw damage
-   - Feedback loop: Progressive deterioration
-
-5. **Stress Corrosion Cracking (SCC) in Pipelines**
-   - External pressure (p) + internal stress (u) → strain energy (Y)
-   - Corrosive environment (c) reduces threshold: Y₀(c) = Y₀·(1-γ_c·c)
-   - Damage evolution: ∂d/∂t ~ Y/η · (1+γ_c·c)
-   - Temperature fluctuations (T) accelerate creep
-   - Stray currents (φ) enhance corrosion
-   - Result: Sudden pipeline failure
-
-6. **Battery Degradation (Capacity Fade)**
-   - Li-ion diffusion (c) → concentration gradients
-   - Volume changes from intercalation (u) → particle stress
-   - Mechanical stress → particle cracking (d↑)
-   - Cracks → impedance increase + side reactions (c)
-   - Heat generation (T) from cycling
-   - Electric field damage (φ) at high charge rates
-   - Result: Capacity fade, thermal runaway risk
-
-7. **Rock Salt Cavern Storage (H₂, CO₂, Natural Gas)**
-   - Creep damage from storage pressure (p → d)
-   - Salt dissolution at interfaces (c) if brine present
-   - Thermal stress from gas temperature (T ≠ T_rock)
-   - Damage → permeability → leakage risk
-   - Electrochemical effects from brines (φ)
-
-═══════════════════════════════════════════════════════════════════════
-THE KEY INSIGHT: DAMAGE IS THE ULTIMATE COUPLING FIELD
-═══════════════════════════════════════════════════════════════════════
-
-Damage doesn't just RESPOND to other fields - it FUNDAMENTALLY CHANGES
-the material properties that govern all other fields!
-
-Traditional approach: Properties are constants
-   κ = 2.0 W/(m·K) ❌ WRONG for damaged material!
-   k = 1e-15 m² ❌ Can change by 5 orders of magnitude!
-   E = 30 GPa ❌ Drops to zero at failure!
-
-JuliaFEM multi-field approach: Properties evolve with damage
-   κ(d) = κ₀·(1 + β_κ·d) ✓ Damage tracked explicitly
-   k(d) = k₀·exp(β_k·d) ✓ Exponential permeability increase
-   C(d) = (1-d)·C₀ ✓ Classic damage mechanics
-
-This is IMPOSSIBLE in traditional single-field or iteratively-coupled FEM!
-
-═══════════════════════════════════════════════════════════════════════
-"""
-
-using Test
-using JuliaFEM
-using LinearAlgebra
-using SparseArrays
-using StaticArrays
-using Tensors
-
-@testset "🌟 THM-ECD: HEXA-PHYSICS (6 Fields!) on All Entity Types" begin
-
-# ═══════════════════════════════════════════════════════════════════════
-# STEP 1: Create 6-field element specification
-# ═══════════════════════════════════════════════════════════════════════
-
-S = @DOFSet{
-    T::DOF{Temperature, Vertex},           # Temperature at vertices
-    u::DOF{Displacement{3}, Vertex},       # Displacement at vertices
-    p::DOF{Pressure, Cell},                # Pressure at cells
-    φ::DOF{ElectricPotential, Edge},       # Electric potential at edges
-    c::DOF{ChemicalConcentration, Vertex},  # Chemical concentration at vertices
-    d::DOF{Damage, Vertex}                 # Damage variable at vertices **[NEW!]**
-}
-
-println("\n" * "="^75)
-println("🌟 THM-ECD: The Ultimate 6-Field Multi-Physics System!")
-println("="^75)
-println("Field specification S includes:")
-println("  1️⃣  T (Temperature) - Float64 at Vertices")
-println("  2️⃣  u (Displacement) - Vec{3,Float64} at Vertices")
-println("  3️⃣  p (Pore Pressure) - Float64 at Cells")
-println("  4️⃣  φ (Electric Potential) - Float64 at Edges")
-println("  5️⃣  c (Chemical Concentration) - Float64 at Vertices")
-println("  6️⃣  d (Damage) - Float64 at Vertices **[NEW!]**")
-println("="^75)
-
-# ═══════════════════════════════════════════════════════════════════════
-# STEP 2: Create element and verify DOF structure
-# ═══════════════════════════════════════════════════════════════════════
-
-# Create a tetrahedron element
-mesh = create_simple_tet_mesh()
-dof_mgr = DOFManager(mesh)
-register_fields!(dof_mgr, S)
-
-elements = create_elements!(dof_mgr, Element{Tetrahedron{4}, Lagrange{1}, S})
-@test length(elements) == 2
-elem = elements[1]
-
-println("\n📦 Element 1: Total local DOFs: ", ndofs(elem))
-
-# Verify field ranges
-T_local = field_dof_range(elem, :T)
-u_local = field_dof_range(elem, :u)
-p_local = field_dof_range(elem, :p)
-φ_local = field_dof_range(elem, :φ)
-c_local = field_dof_range(elem, :c)
-d_local = field_dof_range(elem, :d)  # NEW!
-
-println("  T local range: $T_local ($(length(T_local)) DOFs)")
-println("  u local range: $u_local ($(length(u_local)) DOFs)")
-println("  p local range: $p_local ($(length(p_local)) DOFs)")
-println("  φ local range: $φ_local ($(length(φ_local)) DOFs)")
-println("  c local range: $c_local ($(length(c_local)) DOFs)")
-println("  d local range: $d_local ($(length(d_local)) DOFs) **[NEW!]**")
-
-@test length(T_local) == 4  # 4 vertices
-@test length(u_local) == 12  # 4 vertices × 3 components
-@test length(p_local) == 1  # 1 cell
-@test length(φ_local) == 6  # 6 edges
-@test length(c_local) == 4  # 4 vertices
-@test length(d_local) == 4  # 4 vertices **[NEW!]**
-
-total_local_dofs = length(T_local) + length(u_local) + length(p_local) + 
-                   length(φ_local) + length(c_local) + length(d_local)
-@test total_local_dofs == 31  # 4 + 12 + 1 + 6 + 4 + 4 = 31!
-
-println("\n  ✅ Local DOF verification passed!")
-println("     Total local DOFs per element: $total_local_dofs")
-
-# ═══════════════════════════════════════════════════════════════════════
-# STEP 3: Define inline coupling functions (modular physics!)
-# ═══════════════════════════════════════════════════════════════════════
-
-# Original THM-EC coupling functions (from test_thmec_penta_physics.jl)
-# ... (keeping same as before for brevity - see previous file)
-
-# NEW: Damage coupling functions
-
-"""Thermal conductivity change with damage: κ(d) = κ₀·(1 + β_κ·d)"""
-@inline function thermal_damage_coupling(β_κ, κ, ∇N_T, N_d, vol)
-    return β_κ * κ * (∇N_T ⋅ ∇N_T) * N_d * vol
-end
-
-"""Stiffness degradation with damage: σ = (1-d)·C:ε"""
-@inline function mechanical_damage_coupling(C_eff, ∇N_u, N_d, vol)
-    # This coupling is in the diagonal block of K_uu, not a separate coupling
-    # -C:ε(u):ε(δu) contribution from damage
-    return -C_eff * (∇N_u ⋅ ∇N_u) * N_d * vol
-end
-
-"""Permeability change with damage: k(d) = k₀·exp(β_k·d)"""
-@inline function hydraulic_damage_coupling(β_k, k, ∇N_p, N_d, vol)
-    return β_k * k * (∇N_p ⋅ ∇N_p) * N_d * vol
-end
-
-"""Electrical conductivity change: σ_e(d) = σ_e0·(1 - β_σ·d)"""
-@inline function electric_damage_coupling(β_σ, σ_e, ∇N_φ, N_d, vol)
-    return -β_σ * σ_e * (∇N_φ ⋅ ∇N_φ) * N_d * vol
-end
-
-"""Diffusivity change with damage: D(d) = D₀·(1 + β_D·d)"""
-@inline function chemical_damage_coupling(β_D, D, ∇N_c, N_d, vol)
-    return β_D * D * (∇N_c ⋅ ∇N_c) * N_d * vol
-end
-
-"""Elastic energy drives damage: Y = ½ε:(C:ε)"""
-@inline function damage_driving_force(C_eff, ∇N_u, N_d, vol)
-    # ∂Y/∂ε = C:ε, so coupling is (C:ε(u)) · ε(δd)
-    # Simplified as energy density times basis function
-    strain_energy = 0.5 * C_eff * (∇N_u ⋅ ∇N_u)  # Simplified!
-    return strain_energy * N_d * vol
-end
-
-"""Thermal activation of damage: h_T(T) = exp(β_T^d · T)"""
-@inline function thermal_damage_activation(β_T_d, N_T, N_d, vol)
-    return β_T_d * N_T * N_d * vol
-end
-
-"""Chemically-assisted damage (SCC): h_c(c) = 1 + γ_c·c"""
-@inline function chemical_damage_enhancement(γ_c, N_c, N_d, vol)
-    return γ_c * N_c * N_d * vol
-end
-
-"""Pressure-assisted damage: h_p(p) = 1 + γ_p·p"""
-@inline function pressure_damage_enhancement(γ_p, N_p, N_d, vol)
-    return γ_p * N_p * N_d * vol
-end
-
-"""Electric field damage: h_φ(φ) = 1 + γ_φ·|∇φ|"""
-@inline function electric_damage_enhancement(γ_φ, ∇N_φ, N_d, vol)
-    return γ_φ * sqrt(∇N_φ ⋅ ∇N_φ) * N_d * vol
-end
-
-"""Rate-dependent damage evolution: ∂d/∂t = <Y-Y₀>/η · g(d)"""
-@inline function damage_viscosity_regularization(η_inv, N_d, vol)
-    return η_inv * N_d * N_d * vol
-end
-
-# ═══════════════════════════════════════════════════════════════════════
-# STEP 4: Assemble coupled stiffness matrix (30 coupling blocks!)
-# ═══════════════════════════════════════════════════════════════════════
-
-println("\n🔧 Assembling THM-ECD coupled system...")
-
-# Material parameters (original + damage)
-# ... (same THM-EC parameters as before)
-κ = 2.0  # W/(m·K)
-E = 30e9  # Pa
-ν = 0.3
-k_perm = 1e-15  # m²
-σ_e = 0.01  # S/m
-D_chem = 1e-9  # m²/s
-
-# NEW: Damage parameters
-Y₀ = 1e6  # J/m³
-η = 1e12  # J·s/m³
-β_κ = 2.0  # Thermal conductivity increase
-β_k = 10.0  # Permeability increase (exponential!)
-β_σ = 0.8  # Conductivity decrease
-β_D = 3.0  # Diffusivity increase
-β_T_d = 0.001  # 1/K
-γ_c = 1e-3  # m³/mol
-γ_p = 1e-9  # 1/Pa
-γ_φ = 1e-8  # 1/(V/m)
-η_inv = 1.0 / η
-
-# Global system
-n_total = count_total_dofs(dof_mgr)
-K = spzeros(Float64, n_total, n_total)
-F = zeros(Float64, n_total)
-
-println("  Global system size: ($n_total, $n_total)")
-
-for elem in elements
-    T_local = field_dof_range(elem, :T)
-    u_local = field_dof_range(elem, :u)
-    p_local = field_dof_range(elem, :p)
-    φ_local = field_dof_range(elem, :φ)
-    c_local = field_dof_range(elem, :c)
-    d_local = field_dof_range(elem, :d)  # NEW!
-    
-    n_dofs_local = ndofs(elem)
-    K_local = zeros(n_dofs_local, n_dofs_local)
-    
-    # Integration
-    quad = Gauss{4}()  # Order 4 for accuracy
-    ips = integration_points(quad, Tetrahedron{4}())
-    
-    for ip in ips
-        ξ = Vec{3}(ip.ξ)
-        weight = ip.w
-        
-        # Basis functions and derivatives
-        dN_dξ = get_basis_derivatives(Tetrahedron{4}(), Lagrange{1}(), ξ)
-        N = get_basis_functions(Tetrahedron{4}(), Lagrange{1}(), ξ)
-        
-        # Jacobian (simplified for regular tet)
-        J_det = 1.0 / 6.0
-        vol = J_det * weight
-        
-        # Simplified: dN = dN_dξ (for regular element)
-        dN = dN_dξ
-        
-        # ===== Diagonal Blocks (6 blocks) =====
-        
-        # K_TT: Thermal diffusion (damage-dependent!)
-        for i in 1:4, j in 1:4
-            i_T = T_local[i]
-            j_T = T_local[j]
-            K_local[i_T, j_T] += κ * (dN[i] ⋅ dN[j]) * vol
-        end
-        
-        # K_uu: Mechanical stiffness (damage-dependent!)
-        C_eff = E / ((1 + ν) * (1 - 2*ν))  # Simplified
-        for i in 1:4, j in 1:4
-            for comp_i in 1:3, comp_j in 1:3
-                i_u = u_local[(i-1)*3 + comp_i]
-                j_u = u_local[(j-1)*3 + comp_j]
-                if comp_i == comp_j
-                    K_local[i_u, j_u] += C_eff * (dN[i] ⋅ dN[j]) * vol
-                end
-            end
-        end
-        
-        # K_pp: Hydraulic diffusion (damage-dependent!)
-        K_local[p_local[1], p_local[1]] += k_perm * vol
-        
-        # K_φφ: Electric conduction (damage-dependent!)
-        n_φ = length(φ_local)
-        for i in 1:n_φ, j in 1:n_φ
-            K_local[φ_local[i], φ_local[j]] += σ_e * vol / (n_φ * n_φ)
-        end
-        
-        # K_cc: Chemical diffusion (damage-dependent!)
-        for i in 1:4, j in 1:4
-            i_c = c_local[i]
-            j_c = c_local[j]
-            K_local[i_c, j_c] += D_chem * (dN[i] ⋅ dN[j]) * vol
-        end
-        
-        # K_dd: Damage evolution (rate-dependent) **[NEW!]**
-        for i in 1:4, j in 1:4
-            i_d = d_local[i]
-            j_d = d_local[j]
-            val = damage_viscosity_regularization(η_inv, N[j], vol)
-            K_local[i_d, j_d] += val
-        end
-        
-        # ===== Off-Diagonal Coupling Blocks (30 blocks!) =====
-        
-        # Original 20 THM-EC couplings
-        # ... (same as test_thmec_penta_physics.jl - omitted for brevity)
-        
-        # NEW: 10 Damage-related couplings (5 bidirectional pairs)
-        
-        # K_Td / K_dT: Thermal-damage coupling
-        for i_T in 1:4, i_d in 1:4
-            idx_T = T_local[i_T]
-            idx_d = d_local[i_d]
-            
-            # K_Td: Thermal conductivity increase with damage
-            val_Td = thermal_damage_coupling(β_κ, κ, dN[i_T], N[i_d], vol)
-            K_local[idx_T, idx_d] += val_Td
-            
-            # K_dT: Thermal activation of damage (NOT reciprocal!)
-            val_dT = thermal_damage_activation(β_T_d, N[i_T], N[i_d], vol)
-            K_local[idx_d, idx_T] += val_dT
-        end
-        
-        # K_ud / K_du: Mechanical-damage coupling (THE BIG ONE!)
-        for i_u in 1:4, comp in 1:3, i_d in 1:4
-            idx_u = u_local[(i_u-1)*3 + comp]
-            idx_d = d_local[i_d]
-            
-            # K_ud: Stiffness degradation (1-d)·C:ε
-            val_ud = mechanical_damage_coupling(C_eff, dN[i_u], N[i_d], vol)
-            K_local[idx_u, idx_d] += val_ud
-            
-            # K_du: Elastic energy drives damage
-            val_du = damage_driving_force(C_eff, dN[i_u], N[i_d], vol)
-            K_local[idx_d, idx_u] += val_du
-        end
-        
-        # K_pd / K_dp: Hydraulic-damage coupling (EXPONENTIAL effect!)
-        for i_d in 1:4
-            idx_p = p_local[1]
-            idx_d = d_local[i_d]
-            
-            # K_pd: Permeability increase with damage
-            # Simplified: gradient of p is zero for cell DOF
-            val_pd = β_k * k_perm * N[i_d] * vol
-            K_local[idx_p, idx_d] += val_pd
-            
-            # K_dp: Pressure-assisted damage
-            val_dp = pressure_damage_enhancement(γ_p, N[i_d], N[i_d], vol)
-            K_local[idx_d, idx_p] += val_dp
-        end
-        
-        # K_φd / K_dφ: Electric-damage coupling
-        for i_φ in 1:n_φ, i_d in 1:4
-            idx_φ = φ_local[i_φ]
-            idx_d = d_local[i_d]
-            
-            # K_φd: Conductivity decrease with damage
-            val_φd = -β_σ * σ_e * N[i_d] * vol / n_φ
-            K_local[idx_φ, idx_d] += val_φd
-            
-            # K_dφ: Electric field damage (simplified)
-            val_dφ = γ_φ * vol / n_φ
-            K_local[idx_d, idx_φ] += val_dφ
-        end
-        
-        # K_cd / K_dc: Chemical-damage coupling
-        for i_c in 1:4, i_d in 1:4
-            idx_c = c_local[i_c]
-            idx_d = d_local[i_d]
-            
-            # K_cd: Diffusivity increase with damage
-            val_cd = chemical_damage_coupling(β_D, D_chem, dN[i_c], N[i_d], vol)
-            K_local[idx_c, idx_d] += val_cd
-            
-            # K_dc: Chemically-assisted damage (SCC!)
-            val_dc = chemical_damage_enhancement(γ_c, N[i_c], N[i_d], vol)
-            K_local[idx_d, idx_c] += val_dc
-        end
-    end
-    
-    # Scatter to global
-    println("  📤 Scattering coupled local matrix ($(size(K_local,1))×$(size(K_local,2))) to global")
-    dof_global = get_dof_indices(elem)
-    for i in 1:n_dofs_local, j in 1:n_dofs_local
-        K[dof_global[i], dof_global[j]] += K_local[i, j]
-    end
-end
-
-println("\n✓ Assembly complete!")
-println("  ONE coupled system matrix: $(size(K))")
-println("  Total non-zeros: $(nnz(K))")
-@test nnz(K) > 0
-
-# ═══════════════════════════════════════════════════════════════════════
-# STEP 5: Apply boundary conditions and solve
-# ═══════════════════════════════════════════════════════════════════════
-
-println("\n🎯 Applying boundary conditions...")
-
-# Apply simple BCs (at least 8 constraints for 6 fields)
-# Node 1: Fix T, all u components, c, d
-# Node 2: Fix one u component, φ
-# Cell 1: Fix p
-
-n_T_total, _, _ = count_field_dofs(dof_mgr, :T)
-n_u_total, _, _ = count_field_dofs(dof_mgr, :u)
-n_p_total, _, _ = count_field_dofs(dof_mgr, :p)
-n_φ_total, _, _ = count_field_dofs(dof_mgr, :φ)
-n_c_total, _, _ = count_field_dofs(dof_mgr, :c)
-n_d_total, _, _ = count_field_dofs(dof_mgr, :d)
-
-# Calculate offsets
-offset_T = 0
-offset_u = n_T_total
-offset_p = offset_u + n_u_total
-offset_φ = offset_p + n_p_total
-offset_c = offset_φ + n_φ_total
-offset_d = offset_c + n_c_total
-
-bc_dofs = [
-    offset_T + 1,           # T at node 1
-    offset_u + 1,           # ux at node 1
-    offset_u + 2,           # uy at node 1
-    offset_u + 3,           # uz at node 1
-    offset_u + 4,           # ux at node 2
-    offset_p + 1,           # p at cell 1
-    offset_φ + 1,           # φ at edge 1
-    offset_c + 1,           # c at node 1
-    offset_d + 1            # d at node 1 **[NEW!]**
-]
-
-@test all(bc_dofs .<= n_total)
-
-for dof in bc_dofs
-    K[dof, :] .= 0.0
-    K[:, dof] .= 0.0
-    K[dof, dof] = 1.0
-    F[dof] = 0.0
-end
-
-println("  Applied $(length(bc_dofs)) boundary conditions")
-
-# Add small regularization for stability
-ε_reg = 1e-12
-for i in 1:n_total
-    K[i, i] += ε_reg
-end
-
-println("\n🎯 Solving coupled system...")
-println("   Added regularization (ε=$(ε_reg)) for numerical stability")
-
-sol = K \ F
-println("   ✓ Solution converged!")
-
-# ═══════════════════════════════════════════════════════════════════════
-# STEP 6: Extract and display solution fields
-# ═══════════════════════════════════════════════════════════════════════
-
-println("\n📊 Solution extraction:")
-
-# Extract each field
-T_sol = sol[offset_T+1:offset_u]
-u_sol = sol[offset_u+1:offset_p]
-p_sol = sol[offset_p+1:offset_φ]
-φ_sol = sol[offset_φ+1:offset_c]
-c_sol = sol[offset_c+1:offset_d]
-d_sol = sol[offset_d+1:end]  # NEW!
-
-println("  🌡️  Temperature field: $(length(T_sol)) values")
-println("  🏗️  Displacement field: $(length(u_sol)÷3) nodes × 3 components")
-println("  💧 Pore pressure field: $(length(p_sol)) cells")
-println("  ⚡ Electric potential: $(length(φ_sol)) edges")
-println("  🧪 Chemical concentration field: $(length(c_sol)) nodes")
-println("  💥 Damage field (NEW!): $(length(d_sol)) nodes")
-
-if length(d_sol) > 0
-    println("\n  📈 Damage values:")
-    for (i, d_val) in enumerate(d_sol)
-        println("     Node $i: d = $(round(d_val, sigdigits=4))")
-    end
-end
-
-# ═══════════════════════════════════════════════════════════════════════
-# STEP 7: Verification tests
-# ═══════════════════════════════════════════════════════════════════════
-
-@test length(T_sol) == n_T_total
-@test length(u_sol) == n_u_total
-@test length(p_sol) == n_p_total
-@test length(φ_sol) == n_φ_total
-@test length(c_sol) == n_c_total
-@test length(d_sol) == n_d_total  # NEW!
-
-println("\n" * "="^75)
-println("🎉 HEXA-PHYSICS SUCCESS! 6 Fields Fully Coupled!")
-println("="^75)
-println("✅ ONE element type with SIX physics fields!")
-println("✅ ONE local coupled matrix per element (31×31)")
-println("✅ Total: 30 off-diagonal coupling blocks! (HEXA-PHYSICS!)")
-println("✅ Damage-enhanced thermal conductivity: κ(d) = κ₀·(1+β_κ·d)  🆕")
-println("✅ Damage-degraded stiffness: σ = (1-d)·C:ε  🆕")
-println("✅ Damage-enhanced permeability: k(d) = k₀·exp(β_k·d)  🆕")
-println("✅ Damage-reduced conductivity: σ_e(d) = σ_e0·(1-β_σ·d)  🆕")
-println("✅ Damage-enhanced diffusion: D(d) = D₀·(1+β_D·d)  🆕")
-println("✅ Elastic energy drives damage: Y = ½ε:(C:ε)  🆕")
-println("✅ Thermal damage activation: exp(β_T^d·T)  🆕")
-println("✅ Stress corrosion cracking: γ_c·c  🆕")
-println("✅ Pressure damage: γ_p·p  🆕")
-println("✅ Electric field damage: γ_φ·|∇φ|  🆕")
-println("✅ Rate-dependent damage evolution: η  🆕")
-println("✅ Vertex-based damage DOFs (continuous)  🆕")
-println("✅ Type-safe field access: .T, .u, .p, .φ, .c, .d")
-println("="^75)
-
-end  # @testset
-
-# Helper function (same as before)
-function create_simple_tet_mesh()
-    nodes = [
-        Vec{3}((0.0, 0.0, 0.0)),
-        Vec{3}((1.0, 0.0, 0.0)),
-        Vec{3}((0.0, 1.0, 0.0)),
-        Vec{3}((0.0, 0.0, 1.0)),
-        Vec{3}((1.0, 1.0, 1.0))
-    ]
-    elements = [
-        (1, 2, 3, 4),
-        (2, 3, 4, 5)
-    ]
-    return Mesh(nodes, elements, Tetrahedron{4})
-end