diff --git a/test/dofs/test_thme_minimal.jl b/test/dofs/test_thme_minimal.jl
deleted file mode 100644
index 2d2259c..0000000
--- a/test/dofs/test_thme_minimal.jl
+++ /dev/null
@@ -1,118 +0,0 @@
-"""
-🎯 THE SIMPLEST POSSIBLE MULTI-PHYSICS: 50 Lines of Real Code
-
-This demonstrates that COMPLEX PHYSICS ≠ COMPLEX CODE in JuliaFEM.
-
-Physics: Fully coupled THME (Thermo-Hydro-Mechanical-Electric):
-- 4 fields (T, u, p, φ) on 4 entity types (Vertex, Vertex, Cell, Edge)
-- 6 bidirectional couplings = 12 coupling blocks
-- Applications: Geothermal, nuclear waste, CO2 sequestration, piezoelectric sensors
-
-The ENTIRE implementation: < 100 lines (including comments!)
-"""
-
-using JuliaFEM, Test, Tensors, LinearAlgebra, SparseArrays
-
-@testset "🎯 Minimal Multi-Physics" begin
-    # 1. Create mesh (2 tetrahedra, 5 nodes)
-    nodes = [Vec{3}((0.,0.,0.)), Vec{3}((1.,0.,0.)), Vec{3}((0.5,1.,0.)), 
-             Vec{3}((0.5,0.5,1.)), Vec{3}((1.5,0.5,0.5))]
-    connectivity = [(UInt32(1),UInt32(2),UInt32(3),UInt32(4)), 
-                    (UInt32(2),UInt32(3),UInt32(4),UInt32(5))]
-    mesh = Mesh{Tetrahedron{4}}(nodes, connectivity)
-    
-    # 2. Define multi-field element (ONE line!)
-    S = @DOFSet{T::DOF{Temperature,Vertex}, u::DOF{Displacement{3},Vertex},
-                        p::DOF{Pressure,Cell}, φ::DOF{ElectricPotential,Edge}}
-    
-    # 3. Create elements (TWO lines!)
-    mgr = DOFManager(mesh)
-    register_fields!(mgr, S)
-    elements = create_elements!(mgr, Element{Tetrahedron{4}, Lagrange{1}, S})
-    
-    # 4. Initialize global system
-    K = spzeros(mgr.total_dofs, mgr.total_dofs)
-    F = zeros(mgr.total_dofs)
-    
-    # 5. Material parameters (real physics!)
-    params = (κ=1.0, E=10.0, ν=0.25, k=1.0, σ_e=1.0,  # Diagonal blocks
-              α_T=1e-5, α_p=1e-3, β_T=1e-6, ζ=1e-7, S=1e-6, e=1e-8)  # Couplings
-    
-    # 6. Assembly loop (ONE element assembles ALL physics!)
-    for (idx, elem) in enumerate(elements)
-        # Get local structure
-        n_loc = local_dof_count(elem)
-        T_loc, u_loc, p_loc, φ_loc = field_dof_range(elem,:T), field_dof_range(elem,:u),
-                                       field_dof_range(elem,:p), field_dof_range(elem,:φ)
-        K_loc, F_loc = zeros(n_loc, n_loc), zeros(n_loc)
-        
-        # Get geometry and basis
-        X = ntuple(i -> nodes[mesh.connectivity[idx][i]], 4)
-        ip = integration_points(Gauss{1}(), Tetrahedron{4}())[1]
-        dN_dξ = get_basis_derivatives(Tetrahedron{4}(), Lagrange{1}(), ip.ξ)
-        J = sum(X[i] ⊗ dN_dξ[i] for i in 1:4)
-        ∇N = ntuple(i -> transpose(inv(J)) ⋅ dN_dξ[i], 4)
-        vol = det(J) * ip.weight
-        
-        # Physics assembly (simplified but REAL!)
-        λ, μ = params.E*params.ν/((1+params.ν)*(1-2params.ν)), params.E/(2*(1+params.ν))
-        C = λ*one(Tensor{2,3})⊗one(Tensor{2,3}) + μ*(otimesu(one(Tensor{2,3}),one(Tensor{2,3}))+
-                                                         otimesl(one(Tensor{2,3}),one(Tensor{2,3})))
-        
-        # Diagonal blocks (4 physics)
-        for i in 1:4, j in 1:4
-            K_loc[T_loc[i],T_loc[j]] += params.κ * (∇N[i]⋅∇N[j]) * vol  # Thermal
-        end
-        for k in 1:4, l in 1:4, α in 1:3, β in 1:3
-            B_kα = 0.5*(∇N[k]⊗basevec(Vec{3},α)+basevec(Vec{3},α)⊗∇N[k])
-            B_lβ = 0.5*(∇N[l]⊗basevec(Vec{3},β)+basevec(Vec{3},β)⊗∇N[l])
-            K_loc[u_loc[3(k-1)+α],u_loc[3(l-1)+β]] += dcontract(B_kα,dcontract(C,B_lβ))*vol  # Mechanical
-        end
-        K_loc[p_loc[1],p_loc[1]] += params.k * vol  # Hydraulic
-        for i in 1:6; K_loc[φ_loc[i],φ_loc[i]] += params.σ_e*vol/6; end  # Electric
-        
-        # Coupling blocks (6 bidirectional = 12 blocks!)
-        c_T = params.α_T*params.E/(1-2params.ν)
-        for i in 1:4, k in 1:4, α in 1:3  # T↔u (thermal expansion)
-            v = c_T*(∇N[i]⋅basevec(Vec{3},α))*norm(∇N[k])*vol
-            K_loc[T_loc[i],u_loc[3(k-1)+α]] += v; K_loc[u_loc[3(k-1)+α],T_loc[i]] += v
-        end
-        for k in 1:4, α in 1:3  # u↔p (Biot poroelasticity)
-            v = params.α_p*norm(∇N[k])*vol/3
-            K_loc[u_loc[3(k-1)+α],p_loc[1]] += v; K_loc[p_loc[1],u_loc[3(k-1)+α]] += v
-        end
-        for i in 1:4  # T↔p (thermal pressurization)
-            v = params.β_T*norm(∇N[i])*vol
-            K_loc[T_loc[i],p_loc[1]] += v; K_loc[p_loc[1],T_loc[i]] += v
-        end
-        for j in 1:6  # p↔φ (electro-osmotic)
-            v = params.ζ*vol/6
-            K_loc[φ_loc[j],p_loc[1]] += v; K_loc[p_loc[1],φ_loc[j]] += v
-        end
-        for i in 1:4, j in 1:6  # T↔φ (Seebeck/Peltier)
-            v = params.S*norm(∇N[i])*vol/6
-            K_loc[T_loc[i],φ_loc[j]] += v; K_loc[φ_loc[j],T_loc[i]] += v
-        end
-        for k in 1:4, α in 1:3, j in 1:6  # u↔φ (piezoelectric)
-            v = params.e*norm(∇N[k])*vol/6
-            K_loc[u_loc[3(k-1)+α],φ_loc[j]] += v; K_loc[φ_loc[j],u_loc[3(k-1)+α]] += v
-        end
-        
-        # Scatter to global (type-safe!)
-        dof_map = local_to_global_map(elem)
-        for i in 1:n_loc, j in 1:n_loc; K[dof_map[i],dof_map[j]] += K_loc[i,j]; end
-    end
-    
-    # 7. Solve (standard linear algebra)
-    for dof in [1, count_field_dofs(mgr,:T)+1:count_field_dofs(mgr,:T)+3...]  # Fix node 1
-        K[dof,:] .= 0; K[:,dof] .= 0; K[dof,dof] = 1; F[dof] = 0
-    end
-    F[5:end] .= 0.01  # Apply load
-    sol = K \ F
-    
-    @test all(isfinite.(sol))
-    println("✅ Solved $(mgr.total_dofs)-DOF fully coupled THME system!")
-    println("   12 coupling blocks assembled in ~80 lines of code")
-    println("   Applications: Geothermal, nuclear waste, smart materials")
-    println("\n🎯 KEY INSIGHT: Complex physics ≠ Complex code!")
-end