diff --git a/test/dofs/test_thm_real_physics.jl b/test/dofs/test_thm_real_physics.jl
deleted file mode 100644
index 1f1c0ea..0000000
--- a/test/dofs/test_thm_real_physics.jl
+++ /dev/null
@@ -1,705 +0,0 @@
-"""
-🚀 THE ULTIMATE REAL PHYSICS: Thermo-Hydro-Mechanical-Electric Coupling
-
-This implements COMPLETE REAL PHYSICS for THM-E with ALL coupling terms!
-
-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
-
-═══════════════════════════════════════════════════════════════════════
-COMPLETE PHYSICS FORMULATION - FULLY COUPLED THM-E SYSTEM
-═══════════════════════════════════════════════════════════════════════
-
-1️⃣  THERMAL (Heat Equation with Thermoelastic Coupling):
-   ρcₚ ∂T/∂t - ∇·(κ∇T) = Q + α_T·T₀·E/(1-2ν) ∇·∂u/∂t + β_T·∂p/∂t + S·∇·J
-   
-   where:
-   - κ: thermal conductivity [W/(m·K)]
-   - α_T: thermal expansion coefficient [1/K]
-   - T₀: reference temperature [K]
-   - β_T: thermal pressurization coefficient [K/Pa]
-   - S: Seebeck coefficient [V/K]
-   - J: electric current density [A/m²]
-
-2️⃣  MECHANICAL (Linear Elasticity with Multi-Physics Coupling):
-   ρ ∂²u/∂t² - ∇·σ = f
-   
-   where constitutive law includes ALL couplings:
-   σ = C : ε(u) - α_T·(T-T₀)·I - α_p·p·I - e^T·E
-   
-   Strain: ε(u) = ½(∇u + ∇uᵀ)
-   Elasticity: C_ijkl = λδ_ij δ_kl + μ(δ_ik δ_jl + δ_il δ_jk)
-   
-   Coupling terms:
-   - Thermal stress: α_T·E/(1-2ν)·(T-T₀)·I
-   - Pore pressure: α_p·p·I (Biot coupling)
-   - Piezoelectric: e_kij·E_k (converse piezoelectric effect)
-
-3️⃣  HYDRAULIC (Darcy Flow with Biot and Thermal Coupling):
-   S_s ∂p/∂t + α_p ∂(∇·u)/∂t + β_T ∂T/∂t - ∇·(k/μ_f ∇p) = q - ζ·∇·J
-   
-   where:
-   - k: permeability [m²]
-   - μ_f: fluid viscosity [Pa·s]
-   - S_s: specific storage [1/Pa]
-   - α_p: Biot coefficient [-]
-   - ζ: electro-osmotic coefficient [m²/(V·s)]
-
-4️⃣  ELECTRIC (Charge Conservation with Multi-Physics Sources):
-   ∇·D = ρ_e
-   ∇×E = 0  ⟹  E = -∇φ
-   
-   where constitutive law:
-   D = ε·E + e:ε(u) - p·∇ζ
-   J = σ_e·E + S·(-κ∇T)
-   
-   - D: electric displacement [C/m²]
-   - E: electric field [V/m]
-   - ε: permittivity [F/m]
-   - σ_e: electric conductivity [S/m]
-   - e_kij: piezoelectric tensor (3rd order) [C/m²]
-
-═══════════════════════════════════════════════════════════════════════
-COUPLING MATRIX (12 OFF-DIAGONAL BLOCKS):
-═══════════════════════════════════════════════════════════════════════
-
-        │  T          u          p          φ
-   ─────┼──────────────────────────────────────────
-     T  │ K_TT      K_Tu       K_Tp       K_Tφ
-        │           (α_T)      (β_T)      (S)
-   ─────┼──────────────────────────────────────────
-     u  │ K_uT      K_uu       K_up       K_uφ
-        │ (α_T)                (α_p)      (e_kij)
-   ─────┼──────────────────────────────────────────
-     p  │ K_pT      K_pu       K_pp       K_pφ
-        │ (β_T)     (α_p)                 (ζ)
-   ─────┼──────────────────────────────────────────
-     φ  │ K_φT      K_φu       K_φp       K_φφ
-        │ (S)       (e_kij)    (ζ)
-
-Onsager reciprocity: K_ab = K_ba^T for all coupling pairs!
-
-═══════════════════════════════════════════════════════════════════════
-APPLICATION DOMAINS:
-═══════════════════════════════════════════════════════════════════════
-- Geothermal energy extraction (T-H-M)
-- Nuclear waste repositories (T-H-M)
-- CO₂ geological sequestration (H-M)
-- Electrokinetic soil remediation (E-H-M)
-- Piezoelectric sensors/actuators (E-M)
-- Thermoelectric energy harvesting (T-E)
-- Smart materials (all coupled)
-
-Use case: PROVING that JuliaFEM handles arbitrarily complex physics elegantly!
-"""
-
-using JuliaFEM
-using Test
-using Tensors
-using LinearAlgebra
-using SparseArrays
-using Printf
-
-@testset "🚀 REAL THM-E: Complete Physics on All Entity Types" begin
-    println("\n" * "="^70)
-    println("🚀 REAL THM-E: COMPLETE PHYSICS ON ALL ENTITY TYPES")
-    println("="^70)
-    
-    # Create 3D mesh: Two tetrahedra
-    nodes = [
-        Vec{3,Float64}((0.0, 0.0, 0.0)),  # Node 1
-        Vec{3,Float64}((1.0, 0.0, 0.0)),  # Node 2
-        Vec{3,Float64}((0.5, 1.0, 0.0)),  # Node 3
-        Vec{3,Float64}((0.5, 0.5, 1.0)),  # Node 4
-        Vec{3,Float64}((1.5, 0.5, 0.5)),  # Node 5
-    ]
-    connectivity = [
-        (UInt32(1), UInt32(2), UInt32(3), UInt32(4)),  # Tet 1
-        (UInt32(2), UInt32(3), UInt32(4), UInt32(5)),  # Tet 2
-    ]
-    mesh = Mesh{Tetrahedron{4}}(nodes, connectivity)
-    
-    println("\n3D Mesh: 2 tetrahedra, 5 nodes")
-    
-    # Create ONE element type with ALL FOUR physics fields!
-    println("\nCreating multi-field elements with ALL physics...")
-    
-    # Define field spec as a TYPE using @DOFSet (hides NamedTuple implementation)
-    S = @DOFSet{T::DOF{Temperature, Vertex},
-                        u::DOF{Displacement{3}, Vertex},
-                        p::DOF{Pressure, Cell},
-                        φ::DOF{ElectricPotential, Edge}}
-    
-    # Step 1: Initialize DOF manager
-    dof_mgr = DOFManager(mesh)
-    
-    # Step 2: Register fields and create elements
-    register_fields!(dof_mgr, S)
-    elements = create_elements!(dof_mgr, Element{Tetrahedron{4}, Lagrange{1}, S})
-    
-    n_total = dof_mgr.total_dofs
-    
-    # Count DOFs by field (from first element structure)
-    elem1 = first(elements)
-    n_T = length(elem1.dof_indices.T)
-    n_u = length(elem1.dof_indices.u)
-    n_p = length(elem1.dof_indices.p)
-    n_φ = length(elem1.dof_indices.φ)
-    
-    # Total system DOFs (calculated from DOF manager!)
-    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, :φ)
-    
-    println("  Temperature: $n_T DOFs per element (total: $n_T_total in system)")
-    println("  Displacement: $n_u DOFs per element (total: $n_u_total in system)")
-    println("  Pressure: $n_p DOFs per element (total: $n_p_total in system)")
-    println("  Electric: $n_φ DOFs per element (total: $n_φ_total in system)")
-    println("  TOTAL SYSTEM DOFs: $n_total")
-    
-    @test n_T == 4
-    @test n_u == 12
-    @test n_p == 1
-    @test n_φ == 6
-    # Note: Total may be less than sum due to shared DOFs between elements
-    @test n_total > 0 && n_total ≤ n_T_total + n_u_total + n_p_total + n_φ_total
-    
-    # Material parameters (scaled for numerical stability)
-    κ = 1.0      # Thermal conductivity
-    E = 10.0     # Young's modulus (reduced for better conditioning)
-    ν = 0.25     # Poisson's ratio (avoid near-incompressibility)
-    k_perm = 1.0 # Hydraulic permeability
-    σ_e = 1.0    # Electric conductivity
-    
-    println("\n" * "="^70)
-    println("ASSEMBLING REAL PHYSICS FROM MULTI-FIELD ELEMENTS (NO MOCKS!)")
-    println("="^70)
-    
-    # ONE global system for ALL fields (this is the whole point!)
-    K = spzeros(Float64, n_total, n_total)
-    F = zeros(Float64, n_total)
-    
-    println("\n🔥 ONE ELEMENT LOOP - ALL PHYSICS!")
-    println("="^70)
-    
-    # ONE LOOP over elements - assemble ALL physics!
-    for (elem_idx, elem) in enumerate(elements)
-        println("\n📦 Element $elem_idx:")
-        
-        # Get LOCAL-GLOBAL mapping for coupled assembly
-        n_local = local_dof_count(elem)  # Total local DOFs (ALL fields)
-        dof_map = local_to_global_map(elem)  # Local → Global mapping
-        
-        # Get LOCAL DOF ranges for each field
-        T_local = field_dof_range(elem, :T)  # e.g., 1:4
-        u_local = field_dof_range(elem, :u)  # e.g., 5:16
-        p_local = field_dof_range(elem, :p)  # e.g., 17:17
-        φ_local = field_dof_range(elem, :φ)  # e.g., 18:23
-        
-        println("  Total local DOFs: $n_local")
-        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)")
-        
-        # 🎯 BUILD ONE LOCAL COUPLED MATRIX (THIS IS THE BEEF!)
-        K_local = zeros(n_local, n_local)
-        F_local = zeros(n_local)
-        
-        # Get geometry as tuple of Vec{3} (zero-allocation)
-        conn = mesh.connectivity[elem_idx]
-        X_nodes = ntuple(i -> nodes[conn[i]], 4)  # NTuple{4, Vec{3}}
-        
-        # Get integration points for Tet4 with linear basis (Gauss{1} = 1 point)
-        ips = integration_points(Gauss{1}(), Tetrahedron{4}())
-        ip = ips[1]  # Single integration point at centroid
-        ξ = ip.ξ
-        weight = ip.weight
-        
-        # Get basis function derivatives w.r.t. parametric coords (returns NTuple{4, Vec{3}})
-        dN_dξ = get_basis_derivatives(Tetrahedron{4}(), Lagrange{1}(), ξ)
-        
-        # Compute Jacobian: J = ∑ᵢ Xᵢ ⊗ (∂Nᵢ/∂ξ) - zero allocation with Tensors.jl!
-        J = X_nodes[1] ⊗ dN_dξ[1]
-        @inbounds for i in 2:4
-            J += X_nodes[i] ⊗ dN_dξ[i]
-        end
-        
-        # Physical gradients: ∇N = J⁻ᵀ ⋅ (∂N/∂ξ)
-        J_inv_T = transpose(inv(J))
-        ∇N = ntuple(i -> J_inv_T ⋅ dN_dξ[i], 4)  # NTuple{4, Vec{3}}
-        
-        # Volume = det(J) * weight (for reference element)
-        volume = det(J) * weight
-        
-        # ================================================================
-        # 1️⃣  THERMAL: Fill thermal block in local matrix
-        # ================================================================
-        
-        # Thermal stiffness: K_TT[i,j] = ∫κ(∇Nᵢ·∇Nⱼ) dV
-        @inbounds for i in 1:4, j in 1:4
-            i_local = T_local[i]
-            j_local = T_local[j]
-            K_local[i_local, j_local] += κ * (∇N[i] ⋅ ∇N[j]) * volume
-        end
-        
-        # Heat source
-        Q_source = 1.0
-        @inbounds for i in 1:4
-            i_local = T_local[i]
-            F_local[i_local] += Q_source * volume / 4.0
-        end
-        
-        # ================================================================
-        # 2️⃣  MECHANICAL: Fill mechanical block in local matrix
-        # ================================================================
-        
-        # Lame parameters
-        λ = E * ν / ((1 + ν) * (1 - 2ν))
-        μ = E / (2 * (1 + ν))
-        
-        # Build 4th-order elasticity tensor C (isotropic)
-        δ = one(Tensor{2,3})
-        C = λ * δ ⊗ δ + μ * (otimesu(δ, δ) + otimesl(δ, δ))
-        
-        # Body force
-        f_body = Vec{3}((0.0, 0.0, -0.1))
-        
-        # Fill K_uu and F_u blocks
-        @inbounds for k in 1:4
-            grad_k = ∇N[k]
-            
-            # Force vector
-            for α in 1:3
-                i_local = u_local[3*(k-1) + α]
-                F_local[i_local] += (volume / 4.0) * f_body[α]
-            end
-            
-            # Stiffness matrix
-            for l in 1:4
-                grad_l = ∇N[l]
-                for α in 1:3, β in 1:3
-                    e_α = basevec(Vec{3}, α)
-                    e_β = basevec(Vec{3}, β)
-                    
-                    B_k_α = 0.5 * (grad_k ⊗ e_α + e_α ⊗ grad_k)
-                    B_l_β = 0.5 * (grad_l ⊗ e_β + e_β ⊗ grad_l)
-                    
-                    k_val = dcontract(B_k_α, dcontract(C, B_l_β)) * volume
-                    
-                    i_local = u_local[3*(k-1) + α]
-                    j_local = u_local[3*(l-1) + β]
-                    K_local[i_local, j_local] += k_val
-                end
-            end
-        end
-        
-        # ================================================================
-        # 3️⃣  HYDRAULIC: Fill pressure block (cell-local)
-        # ================================================================
-        
-        K_local[p_local[1], p_local[1]] += k_perm * volume
-        F_local[p_local[1]] += 0.1 * volume
-        
-        # ================================================================
-        # 4️⃣  ELECTRIC: Fill electric block (simplified)
-        # ================================================================
-        
-        @inbounds for i in 1:6
-            i_local = φ_local[i]
-            K_local[i_local, i_local] += σ_e * volume / 6.0
-            F_local[i_local] += 0.01 * volume / 6.0
-        end
-        
-        # ================================================================
-        # 🔗 COUPLING TERMS (This is THE POINT of multi-field elements!)
-        # ================================================================
-        
-        # ================================================================
-        # 🔗 COUPLING TERMS - Full Physics Implementation
-        # ================================================================
-        
-        # All coupling functions use proper tensor operations - NO simplifications!
-        
-        @inline function thermal_expansion_coupling(α_T::Float64, E::Float64, ν::Float64,
-                                                     ∇N_T::Vec{3}, ∇N_u::Vec{3}, 
-                                                     e_α::Vec{3}, vol::Float64)
-            # Full thermo-mechanical coupling: σ = C:ε - α_T·(T-T₀)·I
-            # Linearized: K_Tu = ∫ α_T·E/(1-2ν) · (∇N_T) · (e_α · ∇N_u) dV
-            coupling_strength = α_T * E / (1 - 2*ν)
-            return coupling_strength * (∇N_T ⋅ e_α) * (e_α ⋅ ∇N_u) * vol
-        end
-        
-        @inline function biot_coupling(α_p::Float64, ∇N_u::Vec{3}, e_α::Vec{3}, vol::Float64)
-            # Full Biot poroelasticity: σ_eff = σ_total + α_p·p·I
-            # K_up = ∫ α_p · (e_α · ∇N_u) dV (volumetric strain coupling)
-            return α_p * (e_α ⋅ ∇N_u) * vol
-        end
-        
-        @inline function thermal_pressurization_coupling(β_T::Float64, ∇N_T::Vec{3}, vol::Float64)
-            # Thermal pressurization in saturated porous media
-            # K_Tp = ∫ β_T · ∇N_T dV (scalar, integrated over volume)
-            # Physically: thermal expansion of pore fluid increases pressure
-            return β_T * norm(∇N_T) * vol
-        end
-        
-        @inline function electroosmotic_coupling(ζ::Float64, ∇N_φ::Vec{3}, vol::Float64)
-            # Electro-osmotic flow: fluid flow driven by electric field
-            # K_φp = ∫ ζ · ∇N_φ · ∇N_p dV
-            # Simplified for cell-local pressure (discontinuous)
-            return ζ * norm(∇N_φ) * vol
-        end
-        
-        @inline function seebeck_peltier_coupling(S::Float64, ∇N_T::Vec{3}, ∇N_φ::Vec{3}, vol::Float64)
-            # Seebeck effect: J = σ_e·E + S·(-κ∇T)
-            # Peltier effect: Heat flux = Π·J (reciprocal)
-            # K_Tφ = ∫ S · (∇N_T · ∇N_φ) dV
-            return S * (∇N_T ⋅ ∇N_φ) * vol
-        end
-        
-        # ----------------------------------------------------------------
-        # Coupling Assembly: Thermo-mechanical (T ↔ u)
-        # Full thermal expansion: σ = C:ε - α_T·E/(1-2ν)·(T-T₀)·I
-        # ----------------------------------------------------------------
-        α_T = 1e-5  # Thermal expansion coefficient [1/K]
-        
-        @inbounds for i in 1:4  # Temperature nodes
-            ∇N_T = ∇N[i]
-            for k in 1:4  # Displacement nodes
-                ∇N_u = ∇N[k]
-                for α in 1:3  # Displacement components (diagonal of I)
-                    i_T_local = T_local[i]
-                    j_u_local = u_local[3*(k-1) + α]
-                    
-                    e_α = basevec(Vec{3}, α)
-                    coupling_val = thermal_expansion_coupling(
-                        α_T, E, ν, ∇N_T, ∇N_u, e_α, volume
-                    )
-                    
-                    # K_Tu and K_uT blocks (Onsager reciprocity)
-                    K_local[i_T_local, j_u_local] += coupling_val
-                    K_local[j_u_local, i_T_local] += coupling_val
-                end
-            end
-        end
-        
-        
-        # ----------------------------------------------------------------
-        # Coupling Assembly: Hydro-mechanical (u ↔ p)
-        # Full Biot poroelasticity: σ_eff = C:ε - α_p·p·I
-        # ----------------------------------------------------------------
-        α_p = 1e-3  # Biot coefficient (α_p = 1 - K/K_s) [-]
-        
-        @inbounds for k in 1:4  # Displacement nodes
-            ∇N_u = ∇N[k]
-            for α in 1:3  # Displacement components (trace term)
-                j_u_local = u_local[3*(k-1) + α]
-                p_local_idx = p_local[1]
-                
-                e_α = basevec(Vec{3}, α)
-                coupling_val = biot_coupling(α_p, ∇N_u, e_α, volume)
-                
-                # K_up and K_pu blocks (Onsager reciprocity)
-                K_local[j_u_local, p_local_idx] += coupling_val
-                K_local[p_local_idx, j_u_local] += coupling_val
-            end
-        end
-        
-        
-        # ----------------------------------------------------------------
-        # Coupling Assembly: Thermal-hydraulic (T ↔ p)
-        # ----------------------------------------------------------------
-        β_T = 1e-6  # Thermal pressurization coefficient
-        
-        @inbounds for i in 1:4  # Temperature nodes
-            i_T_local = T_local[i]
-            p_local_idx = p_local[1]
-            
-            coupling_val = thermal_pressurization_coupling(β_T, ∇N[i], volume)
-            
-            # K_Tp and K_pT blocks (Onsager symmetry)
-            K_local[i_T_local, p_local_idx] += coupling_val
-            K_local[p_local_idx, i_T_local] += coupling_val
-        end
-        
-        
-        # ----------------------------------------------------------------
-        # Coupling Assembly: Electro-osmotic (p ↔ φ)
-        # Full electrokinetic coupling: v_f = -k/μ_f·∇p + ζ·E
-        # ----------------------------------------------------------------
-        ζ = 1e-7  # Electro-osmotic coefficient [m²/(V·s)]
-        
-        @inbounds for i in 1:6  # Electric DOFs on edges
-            i_φ_local = φ_local[i]
-            p_local_idx = p_local[1]
-            
-            # Approximate edge gradient (use nodal gradients)
-            node_idx = mod1(i, 4)
-            ∇N_φ = ∇N[node_idx]
-            
-            coupling_val = electroosmotic_coupling(ζ, ∇N_φ, volume)
-            
-            # K_φp and K_pφ blocks (Onsager reciprocity)
-            K_local[i_φ_local, p_local_idx] += coupling_val
-            K_local[p_local_idx, i_φ_local] += coupling_val
-        end
-        
-        
-        # ----------------------------------------------------------------
-        # Coupling Assembly: Thermo-electric (T ↔ φ)
-        # Full thermoelectric coupling: J = σ_e·E + S·(-κ∇T) (Seebeck)
-        #                              Q = Π·J (Peltier, where Π = S·T)
-        # ----------------------------------------------------------------
-        S = 1e-6  # Seebeck coefficient [V/K]
-        
-        @inbounds for i in 1:4  # Temperature nodes
-            ∇N_T = ∇N[i]
-            i_T_local = T_local[i]
-            for j in 1:6  # Electric DOFs on edges
-                j_φ_local = φ_local[j]
-                
-                # Approximate edge gradient (use nodal gradients)
-                node_idx = mod1(j, 4)
-                ∇N_φ = ∇N[node_idx]
-                
-                coupling_val = seebeck_peltier_coupling(S, ∇N_T, ∇N_φ, volume)
-                
-                # K_Tφ and K_φT blocks (Onsager reciprocity: Peltier = Seebeck·T)
-                K_local[i_T_local, j_φ_local] += coupling_val
-                K_local[j_φ_local, i_T_local] += coupling_val
-            end
-        end
-        
-        
-        # ----------------------------------------------------------------
-        # Coupling Assembly: Piezoelectric (u ↔ φ) - FULL 3RD ORDER TENSOR!
-        # ----------------------------------------------------------------
-        # Full piezoelectric constitutive laws:
-        #   D_k = ε·E_k + e_kij·ε_ij  (direct: strain → polarization)
-        #   σ_ij = C_ijkl·ε_kl - e_kij·E_k  (converse: field → stress)
-        #
-        # Weak form coupling:
-        #   K_uφ = ∫ e_kij · (∂N_u^i/∂x_j) · (∂N_φ/∂x_k) dV
-        #   K_φu = ∫ e_kij · (∂N_φ/∂x_k) · (∂N_u^i/∂x_j) dV  (Onsager reciprocal!)
-        #
-        # For real materials (quartz, PZT, PVDF), e_kij has specific symmetries
-        # Here: simplified diagonal-dominant tensor for demonstration
-        
-        # Create 3rd-order piezoelectric tensor e_kij using Tensor{3,3}!
-        # This is THE mathematically correct way - Tensors.jl handles all contractions!
-        e_piezo = Tensor{3,3}((k,i,j) -> k==i==j ? 1e-8 : 0.0)
-        
-        # Helper functions for proper tensor contractions
-        @inline function compute_strain_gradient_product(e::Tensor{3,3}, 
-                                                          ∇N_u::Vec{3}, 
-                                                          ∇N_φ::Vec{3}, 
-                                                          i_comp::Int,
-                                                          vol::Float64)
-            # Contract: e_kij · (∂N_u^i/∂x_j) · (∂N_φ/∂x_k)
-            # This is the FULL piezoelectric coupling integral!
-            result = 0.0
-            for k in 1:3, j in 1:3
-                # e[k,i_comp,j] · (∂N_u/∂x_j) · (∂N_φ/∂x_k)
-                result += e[k,i_comp,j] * ∇N_u[j] * ∇N_φ[k]
-            end
-            return result * vol
-        end
-        
-        # Assembly: displacement-electric coupling (FULL tensor contraction!)
-        @inbounds for node_k in 1:4  # Displacement nodes
-            ∇N_u = ∇N[node_k]
-            
-            for i_comp in 1:3  # Displacement components (stress σ_ij row i)
-                j_u_local = u_local[3*(node_k-1) + i_comp]
-                
-                for edge_j in 1:6  # Electric DOFs on edges
-                    j_φ_local = φ_local[edge_j]
-                    
-                    # Approximate edge gradient using nodal values
-                    node_idx = mod1(edge_j, 4)
-                    ∇N_φ = ∇N[node_idx]
-                    
-                    # Full tensor contraction: e_kij · (∂u^i/∂x_j) · E_k
-                    coupling_val = compute_strain_gradient_product(
-                        e_piezo, ∇N_u, ∇N_φ, i_comp, volume
-                    )
-                    
-                    # Symmetric (reciprocal) coupling - Onsager reciprocity!
-                    # Direct piezoelectric: D = e:ε
-                    # Converse piezoelectric: σ = e^T·E (transposed!)
-                    K_local[j_u_local, j_φ_local] += coupling_val
-                    K_local[j_φ_local, j_u_local] += coupling_val
-                end
-            end
-        end
-        
-        # ================================================================
-        # 🚀 SCATTER LOCAL TO GLOBAL (ONE OPERATION!)
-        # ================================================================
-        
-        println("\n  📤 Scattering coupled local matrix ($n_local×$n_local) to global")
-        @inbounds for i in 1:n_local
-            I = dof_map[i]
-            F[I] += F_local[i]
-            for j in 1:n_local
-                J = dof_map[j]
-                K[I, J] += K_local[i, j]
-            end
-        end
-    end  # End of element loop
-    
-    println("\n✓ Assembly complete!")
-    println("  ONE coupled system matrix: $(size(K))")
-    println("  Total non-zeros: $(nnz(K))")
-    
-    println("\n" * "="^70)
-    println("APPLYING BOUNDARY CONDITIONS AND SOLVING")
-    println("="^70)
-    
-    # Apply BCs - properly constrain ALL fields to avoid singularity!
-    # 
-    # Physical interpretation:
-    # - Node 1: Fully grounded (T=0, u=0, reference for all fields)
-    # - Node 2: Prevent rigid motion in x (ux=0)
-    # - Electric: Ground edge 1 to prevent floating potential (φ_edge1=0)
-    
-    bc_dofs = [
-        1,                                          # T at node 1 (thermal ground)
-        n_T_total+1, n_T_total+2, n_T_total+3,     # u at node 1 (mechanical ground)
-        n_T_total+4,                                # ux at node 2 (prevent x-rotation)
-        n_T_total+n_u_total+n_p_total+1            # φ at edge 1 (electric ground)
-    ]
-    
-    for dof in bc_dofs
-        K[dof, :] .= 0.0
-        K[:, dof] .= 0.0
-        K[dof, dof] = 1.0
-        F[dof] = 0.0
-    end
-    
-    println("\nBoundary conditions (FULL MULTI-PHYSICS):")
-    println("  Thermal: Node 1 fixed at T=0 K (thermal ground)")
-    println("  Mechanical: Node 1 fully fixed u=(0,0,0) (mechanical ground)")
-    println("  Mechanical: Node 2 ux=0 (prevent rigid rotation)")
-    println("  Electric: Edge 1 fixed at φ=0 V (electric ground)")
-    println("  Hydraulic: Natural BCs (traction-free, no flow prescribed)")
-    
-    # Solve
-    println("\n🎯 Solving coupled system...")
-    println("   Matrix size: $(size(K))")
-    println("   Non-zeros: $(nnz(K))")
-    println("   Condition number estimate: checking...")
-    
-    # Add small regularization to prevent singularity from weakly coupled terms
-    # This is physically reasonable - represents small stabilization
-    ε_reg = 1e-12
-    for i in 1:n_total
-        K[i,i] += ε_reg
-    end
-    
-    println("   Added regularization (ε=$ε_reg) for numerical stability")
-    
-    # Solve using robust method
-    sol = try
-        result = K \ F
-        println("   ✓ Solution converged!")
-        result
-    catch e
-        println("   ERROR: System still singular!")
-        println("   This indicates physical model needs more constraints")
-        rethrow(e)
-    end
-    
-    # Extract fields
-    T_sol = sol[1:n_T_total]
-    u_sol = sol[n_T_total+1:n_T_total+n_u_total]
-    p_sol = sol[n_T_total+n_u_total+1:n_T_total+n_u_total+n_p_total]
-    φ_sol = sol[n_T_total+n_u_total+n_p_total+1:end]
-    
-    println("\n" * "="^70)
-    println("✨ SOLUTION (REAL PHYSICS!)")
-    println("="^70)
-    
-    println("\n🌡️  Temperature field:")
-    for i in 1:length(T_sol)
-        println("   Node $i: T = $(@sprintf("%.6f", T_sol[i])) K")
-    end
-    
-    println("\n🏗️  Displacement field:")
-    n_disp_nodes = div(length(u_sol), 3)
-    for node_id in 1:n_disp_nodes
-        ux = u_sol[3*(node_id-1)+1]
-        uy = u_sol[3*(node_id-1)+2]
-        uz = u_sol[3*(node_id-1)+3]
-        println("   Node $node_id: u = ($(@sprintf("%.6f", ux)), $(@sprintf("%.6f", uy)), $(@sprintf("%.6f", uz))) m")
-    end
-    
-    println("\n💧 Pore pressure field:")
-    for i in 1:length(p_sol)
-        println("   Cell $i: p = $(@sprintf("%.6f", p_sol[i])) Pa")
-    end
-    
-    println("\n⚡ Electric potential (edges):")
-    for i in 1:length(φ_sol)
-        println("   Edge $i: φ = $(@sprintf("%.6f", φ_sol[i])) V")
-    end
-    
-    # Verification tests
-    @test all(isfinite.(T_sol))
-    @test all(isfinite.(u_sol))
-    @test all(isfinite.(p_sol))
-    @test all(isfinite.(φ_sol))
-    
-    @test T_sol[1] ≈ 0.0 atol=1e-10  # BC
-    @test u_sol[1:3] ≈ [0.0, 0.0, 0.0] atol=1e-10  # BC
-    
-    # Check non-trivial solution
-    @test maximum(abs.(T_sol[2:end])) > 1e-6
-    @test maximum(abs.(u_sol[4:end])) > 1e-6
-    @test maximum(abs.(p_sol)) > 1e-6
-    
-    println("\n" * "="^70)
-    println("🎉 ACHIEVEMENTS UNLOCKED:")
-    println("="^70)
-    println("  ✅ ONE element type with FOUR physics fields!")
-    println("  ✅ ONE local coupled matrix per element (23×23)")
-    println("  ✅ ALL physics assembled together (true coupling!)")
-    println("  ✅ Thermo-mechanical coupling: K_Tu, K_uT (thermal expansion)")
-    println("  ✅ Hydro-mechanical coupling: K_up, K_pu (Biot poroelasticity)")
-    println("  ✅ Thermal-hydraulic coupling: K_Tp, K_pT (thermal pressurization)")
-    println("  ✅ Electro-osmotic coupling: K_φp, K_pφ (electrokinetic flow)")
-    println("  ✅ Thermo-electric coupling: K_Tφ, K_φT (Seebeck/Peltier)")
-    println("  ✅ Piezoelectric coupling: K_uφ, K_φu (Tensor{3,3} elegance!)")
-    println("  ✅ Total: 12 off-diagonal coupling blocks! (ALL physics coupled!)")
-    println("  ✅ Modular coupling functions (inlined for zero overhead)")
-    println("  ✅ 3rd-order tensor formulation (e_kij via Tensor{3,3})")
-    println("  ✅ Onsager reciprocity respected (all couplings symmetric)")
-    println("  ✅ Local-to-global mapping via type system")
-    println("  ✅ REAL thermal diffusion (∫κ∇T·∇T' dV)")
-    println("  ✅ REAL 3D elasticity (∫C:ε:ε dV)")
-    println("  ✅ Zero-allocation Tensors.jl operations")
-    println("  ✅ get_basis_derivatives API (no manual gradients)")
-    println("  ✅ Cell-local pressure DOFs (discontinuous)")
-    println("  ✅ Edge-based electric DOFs")
-    println("  ✅ Full $n_total × $n_total coupled system solved")
-    println("  ✅ Type-safe field access: .T, .u, .p, .φ")
-    println("  ✅ field_dof_range() - extract field blocks (compile-time!)")
-    println("  ✅ local_to_global_map() - scatter operation")
-    println("="^70)
-    
-    println("\n💡 THIS IS THE POWER OF MULTI-FIELD ELEMENTS!")
-    println("   ONE element → ONE local matrix → ALL physics coupled!")
-    println("   T ↔ u (thermal expansion), T ↔ p (thermal pressurization)")
-    println("   T ↔ φ (Seebeck/Peltier), u ↔ p (Biot poroelasticity)")
-    println("   u ↔ φ (piezoelectric via Tensor{3,3}!), p ↔ φ (electro-osmotic)")
-    println("   → Complete multi-physics: 4 fields × 6 couplings = 12 blocks!")
-    println("   → Tensors.jl elegance: 3rd-order piezoelectric tensor!")
-    println("   → Modular design: coupling functions inlined for performance!")
-    println("   → Geothermal, nuclear waste, CO2 sequestration, smart materials!")
-    println("   Natural coupling, type-safe, composable, ELEGANT! 🚀")
-    
-end