diff --git a/test/elements/test_interpolate.jl b/test/elements/test_interpolate.jl
deleted file mode 100644
index bde74e1..0000000
--- a/test/elements/test_interpolate.jl
+++ /dev/null
@@ -1,234 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-using Test
-using JuliaFEM
-using Tensors
-
-@testset "Interpolation: Single-field Scalar" begin
-    # Thermal element: 1 scalar per node
-    S = @DOFSet{T::DOF{Temperature, Vertex}}
-    elem = Element{Triangle{3}, Lagrange{1}, S}(
-        UInt(1),
-        (UInt64(1), UInt64(2), UInt64(3))  # Flat tuple of DOF indices
-    )
-    
-    # Global solution: T values at nodes = [10.0, 20.0, 30.0]
-    u_global = zeros(100)
-    u_global[1] = 10.0
-    u_global[2] = 20.0
-    u_global[3] = 30.0
-    
-    # Interpolate at reference center (1/3, 1/3) in barycentric
-    # For linear triangle in reference coords (-1,-1) to (1,1), center is (0,0)
-    ξ = Vec((0.0, 0.0))
-    
-    vals = interpolate_fields(elem, u_global, ξ)
-    @test vals isa NamedTuple
-    @test haskey(vals, :T)
-    @test haskey(vals, :∇T)
-    
-    # At center of linear triangle, value should be average
-    # For Tri3 with vertices at (-1,-1), (1,-1), (-1,1),
-    # N1(0,0) = N2(0,0) = N3(0,0) for equilateral triangle
-    @test vals.T isa Float64
-    @test vals.∇T isa Vec{2}
-    
-    # Test single field extraction
-    T_val, T_grad = interpolate_field(elem, u_global, :T, ξ)
-    @test T_val ≈ vals.T
-    @test T_grad ≈ vals.∇T
-    
-    # Test value-only extraction
-    T_only = interpolate_field_value(elem, u_global, :T, ξ)
-    @test T_only ≈ vals.T
-end
-
-@testset "Interpolation: Single-field Vector 2D" begin
-    # 2D elasticity on triangle
-    S = @DOFSet{u::DOF{Displacement{2}, Vertex}}
-    elem = Element{Triangle{3}, Lagrange{1}, S}(
-        UInt(1),
-        tuple(UInt64.(1:6)...)  # ux1, uy1, ux2, uy2, ux3, uy3
-    )
-    
-    # Global solution: displacement vectors at nodes
-    u_global = zeros(100)
-    u_global[1:6] = [1.0, 0.0, 0.0, 1.0, 1.0, 1.0]  # Node 1: (1,0), Node 2: (0,1), Node 3: (1,1)
-    
-    ξ = Vec((0.0, 0.0))
-    
-    vals = interpolate_fields(elem, u_global, ξ)
-    @test vals isa NamedTuple
-    @test haskey(vals, :u)
-    @test haskey(vals, :∇u)
-    @test vals.u isa Vec{2}
-    @test vals.∇u isa Tensor{2,2}  # ∇u is 2x2 gradient tensor
-    
-    # Test single field
-    u_val, u_grad = interpolate_field(elem, u_global, :u, ξ)
-    @test u_val isa Vec{2}
-    @test u_grad isa Tensor{2,2}
-    @test u_val ≈ vals.u
-    @test u_grad ≈ vals.∇u
-    
-    # Test value-only
-    u_only = interpolate_field_value(elem, u_global, :u, ξ)
-    @test u_only isa Vec{2}
-    @test u_only ≈ vals.u
-end
-
-@testset "Interpolation: Single-field Vector 3D" begin
-    # 3D elasticity on tetrahedron
-    S = @DOFSet{u::DOF{Displacement{3}, Vertex}}
-    elem = Element{Tetrahedron{4}, Lagrange{1}, S}(
-        UInt(1),
-        tuple(UInt64.(1:12)...)  # 4 nodes × 3 components
-    )
-    
-    # Set up a simple displacement field
-    u_global = zeros(100)
-    u_global[1:12] = [
-        1.0, 0.0, 0.0,  # Node 1: (1,0,0)
-        0.0, 1.0, 0.0,  # Node 2: (0,1,0)
-        0.0, 0.0, 1.0,  # Node 3: (0,0,1)
-        1.0, 1.0, 1.0   # Node 4: (1,1,1)
-    ]
-    
-    # Interpolate at reference center
-    ξ = Vec((0.25, 0.25, 0.25))
-    
-    vals = interpolate_fields(elem, u_global, ξ)
-    @test vals isa NamedTuple
-    @test vals.u isa Vec{3}
-    @test vals.∇u isa Tensor{2,3}  # ∇u is 3x3 gradient tensor
-    
-    # Test consistency
-    u_val, u_grad = interpolate_field(elem, u_global, :u, ξ)
-    @test u_val ≈ vals.u
-    @test u_grad ≈ vals.∇u
-end
-
-@testset "Interpolation: Multi-field (Temperature + Displacement)" begin
-    # Thermo-mechanical element
-    S = @DOFSet{T::DOF{Temperature,Vertex}, u::DOF{Displacement{3},Vertex}}
-    elem = Element{Tetrahedron{4}, Lagrange{1}, S}(
-        UInt(1),
-        tuple(UInt64.([1, 2, 3, 4, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21])...)
-    )
-    
-    # Global solution with both fields
-    u_global = zeros(100)
-    u_global[1:4] = [300.0, 350.0, 400.0, 375.0]  # Temperature
-    u_global[10:21] = [
-        0.1, 0.0, 0.0,  # Node 1 displacement
-        0.0, 0.1, 0.0,  # Node 2 displacement
-        0.0, 0.0, 0.1,  # Node 3 displacement
-        0.1, 0.1, 0.1   # Node 4 displacement
-    ]
-    
-    ξ = Vec((0.25, 0.25, 0.25))
-    
-    # Test interpolate_fields (all fields at once)
-    vals = interpolate_fields(elem, u_global, ξ)
-    @test vals isa NamedTuple
-    @test haskey(vals, :T)
-    @test haskey(vals, :∇T)
-    @test haskey(vals, :u)
-    @test haskey(vals, :∇u)
-    @test vals.T isa Float64
-    @test vals.∇T isa Vec{3}
-    @test vals.u isa Vec{3}
-    @test vals.∇u isa Tensor{2,3}
-    
-    # Temperature should be in reasonable range
-    @test 300.0 ≤ vals.T ≤ 400.0
-    
-    # Test individual field extraction
-    T_val, T_grad = interpolate_field(elem, u_global, :T, ξ)
-    @test T_val ≈ vals.T
-    @test T_grad ≈ vals.∇T
-    
-    u_val, u_grad = interpolate_field(elem, u_global, :u, ξ)
-    @test u_val ≈ vals.u
-    @test u_grad ≈ vals.∇u
-    
-    # Test value-only extraction
-    T_only = interpolate_field_value(elem, u_global, :T, ξ)
-    @test T_only ≈ vals.T
-    
-    u_only = interpolate_field_value(elem, u_global, :u, ξ)
-    @test u_only ≈ vals.u
-end
-
-@testset "Zero Allocation: Interpolation functions" begin
-    # Single-field scalar
-    S_scalar = @DOFSet{T::DOF{Temperature, Vertex}}
-    elem_scalar = Element{Triangle{3}, Lagrange{1}, S_scalar}(
-        UInt(1),
-        (UInt64(1), UInt64(2), UInt64(3))
-    )
-    u_global = rand(100)
-    ξ2 = Vec((0.0, 0.0))
-    
-    # Warm up
-    interpolate_fields(elem_scalar, u_global, ξ2)
-    interpolate_field(elem_scalar, u_global, :T, ξ2)
-    interpolate_field_value(elem_scalar, u_global, :T, ξ2)
-    
-    # Check zero allocation
-    alloc1 = @allocated interpolate_fields(elem_scalar, u_global, ξ2)
-    alloc2 = @allocated interpolate_field(elem_scalar, u_global, :T, ξ2)
-    alloc3 = @allocated interpolate_field_value(elem_scalar, u_global, :T, ξ2)
-    @test alloc1 == 0
-    @test alloc2 == 0
-    @test alloc3 == 0
-    
-    # Single-field vector 3D
-    S_vec3d = @DOFSet{u::DOF{Displacement{3}, Vertex}}
-    elem_vec3d = Element{Tetrahedron{4}, Lagrange{1}, S_vec3d}(
-        UInt(1),
-        tuple(UInt64.(1:12)...)
-    )
-    ξ3 = Vec((0.25, 0.25, 0.25))
-    
-    # Warm up
-    interpolate_fields(elem_vec3d, u_global, ξ3)
-    interpolate_field(elem_vec3d, u_global, :u, ξ3)
-    interpolate_field_value(elem_vec3d, u_global, :u, ξ3)
-    
-    # Check zero allocation
-    alloc1 = @allocated interpolate_fields(elem_vec3d, u_global, ξ3)
-    alloc2 = @allocated interpolate_field(elem_vec3d, u_global, :u, ξ3)
-    alloc3 = @allocated interpolate_field_value(elem_vec3d, u_global, :u, ξ3)
-    @test alloc1 == 0
-    @test alloc2 == 0
-    @test alloc3 == 0
-    
-    # Multi-field
-    S_multi = @DOFSet{T::DOF{Temperature,Vertex}, u::DOF{Displacement{3},Vertex}}
-    elem_multi = Element{Tetrahedron{4}, Lagrange{1}, S_multi}(
-        UInt(1),
-        tuple(UInt64.([1, 2, 3, 4, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21])...)
-    )
-    
-    # Warm up
-    interpolate_fields(elem_multi, u_global, ξ3)
-    interpolate_field(elem_multi, u_global, :T, ξ3)
-    interpolate_field(elem_multi, u_global, :u, ξ3)
-    
-    # Check zero allocation
-    alloc1 = @allocated interpolate_fields(elem_multi, u_global, ξ3)
-    alloc2 = @allocated interpolate_field(elem_multi, u_global, :T, ξ3)
-    alloc3 = @allocated interpolate_field(elem_multi, u_global, :u, ξ3)
-    @test alloc1 == 0
-    @test alloc2 == 0
-    @test alloc3 == 0
-end
-
-println("✓ All interpolation tests passed!")
-println("  - Single-field scalar interpolation works")
-println("  - Single-field vector interpolation works (2D and 3D)")
-println("  - Multi-field interpolation works")
-println("  - Field values and gradients computed correctly")
-println("  - Zero allocations verified for all cases")