diff --git a/test/geometry/test_deformation_gradient.jl b/test/geometry/test_deformation_gradient.jl
deleted file mode 100644
index 879e359..0000000
--- a/test/geometry/test_deformation_gradient.jl
+++ /dev/null
@@ -1,561 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-"""
-# Deformation Gradient Tests (test/geometry/)
-
-## What
-Tests computation of the deformation gradient **F = ∂x/∂X = I + ∂u/∂X**, which maps
-from reference (undeformed) to current (deformed) configuration in finite strain mechanics.
-
-## Why
-The deformation gradient is **FUNDAMENTAL** to finite strain analysis:
-- **Green-Lagrange strain**: E = ½(F^T F - I)
-- **Cauchy-Green deformation tensor**: C = F^T F
-- **Volume ratio**: J = det(F) (incompressibility requires J = 1)
-- **Polar decomposition**: F = RU (rotation × stretch)
-- **Material frame**: All stress/strain in reference configuration
-
-Physical requirements:
-- **det(F) > 0**: No material inversion (orientation preserved)
-- **det(F) ≈ 1**: Nearly incompressible materials (rubber, metal plasticity)
-- **F = I**: Undeformed configuration (u = 0)
-
-This test validates:
-- **Identity case**: u = 0 → F = I, det(F) = 1
-- **Pure translation**: ∇u = 0 → F = I (rigid body motion)
-- **Pure stretch**: Diagonal F with stretch ratios
-- **Simple shear**: Off-diagonal F terms
-- **Small vs finite strain**: Difference in formulations
-- **Physical constraints**: det(F) > 0 always
-- **Zero allocations**: Hot path allocates nothing
-
-## How
-**Test Cases:**
-
-**1. Identity (u = 0)**:
-- Unit cube Hex8 element, zero displacement
-- Expected: F = I, det(F) = 1.0
-- Both finite strain and small strain give same result
-
-**2. Pure Translation**:
-- Uniform displacement u = (0.5, 0.5, 0.5) at all nodes
-- Expected: ∇u = 0 → F = I + 0 = I
-- Validates that rigid body motion doesn't deform material
-
-**3. Pure Stretch (x-direction)**:
-- Displacement u_x = 0.1·X (10% engineering strain)
-- Expected: F = diag(1.1, 1.0, 1.0), det(F) = 1.1
-- Validates uniaxial extension
-
-**4. Simple Shear**:
-- Displacement u_x = 0.1·y (shear deformation)
-- Expected: F_{12} = 0.1, det(F) = 1.0 (volume preserving)
-- Validates shear kinematics
-
-**5. Small vs Finite Strain Difference**:
-- 20% stretch: significant displacement gradient
-- **Finite strain**: F = I + ∇u (includes gradient)
-- **Small strain**: F = I (ignores gradient, approximation)
-- Validates that formulations differ for large deformations
-
-**6. Physical Constraint**:
-- All physical deformations must have det(F) > 0
-- Negative det(F) → inverted element (unphysical)
-
-**7. Tet10 (Quadratic) Elements**:
-- Higher-order elements with 10 nodes
-- Same tests as Hex8 but with quadratic basis
-- Validates that API works for all element types
-
-**8. Zero Allocations**:
-- `@allocated compute_deformation_gradient(...)` must return 0
-- Critical for performance in assembly loops
-
-## Expected Results
-- ✅ **Identity**: F = I, det(F) = 1.0
-- ✅ **Translation**: F = I (∇u = 0)
-- ✅ **Stretch**: F_{ii} = 1 + ε_{ii}, det(F) = product of stretches
-- ✅ **Shear**: Off-diagonal terms non-zero, det(F) = 1.0
-- ✅ **Finite ≠ Small**: Different F for large deformations (20%+)
-- ✅ **Physical**: det(F) > 0 always
-- ✅ **Tet10**: Works with quadratic elements
-- ✅ **Zero allocations**: @allocated = 0
-
-## Mathematical Background
-**Deformation Gradient (Finite Strain):**
-```
-F = ∂x/∂X = I + ∂u/∂X
-```
-
-Where:
-- x = X + u (current position = reference + displacement)
-- X = reference coordinates
-- u = displacement vector
-- ∂u/∂X = displacement gradient
-
-**Small Strain Approximation:**
-```
-F ≈ I   (ignores ∂u/∂X, valid for ||∇u|| << 1)
-```
-
-**Computation:**
-```julia
-# Displacement gradient
-∇u = ∑ᵢ uᵢ ⊗ (∂Nᵢ/∂X)
-
-# Physical derivatives
-∂Nᵢ/∂X = J⁻¹ · ∂Nᵢ/∂ξ
-
-# Deformation gradient
-F = I + ∇u
-```
-
-## Formulation Comparison
-**Finite Strain** (use when ||∇u|| > 0.01):
-- F = I + ∇u (full nonlinear kinematics)
-- E = ½(F^T F - I) (Green-Lagrange strain)
-- S = ∂W/∂E (2nd Piola-Kirchhoff stress)
-- Required for: Rubber, large rotations, metal forming
-
-**Small Strain** (use when ||∇u|| < 0.01):
-- F ≈ I (linear approximation)
-- ε = ½(∇u + ∇u^T) (engineering strain)
-- σ = C:ε (Cauchy stress)
-- Valid for: Linear elasticity, small vibrations
-
-## Architecture Principle
-**Tensors.jl for Kinematics**
-
-All kinematic quantities use Tensors.jl:
-- Displacement: Vec{3,Float64}
-- Gradient: Tensor{2,3} (⊗ outer product)
-- Deformation gradient: Tensor{2,3}
-- Identity: one(Tensor{2,3})
-
-This provides:
-- Natural tensor notation (F[i,j])
-- Automatic differentiation ready
-- Zero-allocation operations
-- GPU compatible
-
-## Usage Pattern
-```julia
-# In assembly loop:
-for ip in integration_points
-    # Get basis derivatives
-    dN_dξ = get_basis_derivatives(topology, basis, ip.ξ)
-    
-    # Compute Jacobian
-    J = compute_jacobian(X_nodes, dN_dξ)
-    
-    # Compute deformation gradient
-    F = compute_deformation_gradient(X_nodes, u_nodes, dN_dξ, J, FiniteStrain())
-    
-    # Compute strain
-    E = 0.5 * (F' ⊡ F - I)
-    
-    # Material model
-    S, 𝔻, state = compute_stress(material, E, state_old, Δt)
-    
-    # ... assemble
-end
-```
-
-## Critical Performance Path
-Deformation gradient is computed at EVERY integration point, EVERY Newton iteration.
-For 100k elements × 8 integration points × 10 Newton iterations = 8M calls!
-
-Zero allocations are MANDATORY!
-"""
-
-using Test
-using Tensors
-using LinearAlgebra
-
-# Use JuliaFEM for basis functions
-using JuliaFEM
-
-# Load our new deformation gradient code
-include("../src/physics/deformation_gradient.jl")
-
-@testset "Deformation Gradient - Low Level API" begin
-
-    @testset "Identity case (u = 0)" begin
-        # Unit cube element, no displacement
-        X_nodes = (
-            Vec(0.0, 0.0, 0.0),
-            Vec(1.0, 0.0, 0.0),
-            Vec(1.0, 1.0, 0.0),
-            Vec(0.0, 1.0, 0.0),
-            Vec(0.0, 0.0, 1.0),
-            Vec(1.0, 0.0, 1.0),
-            Vec(1.0, 1.0, 1.0),
-            Vec(0.0, 1.0, 1.0)
-        )
-
-        # Zero displacement
-        u_nodes = tuple([zero(Vec{3,Float64}) for _ in 1:8]...)
-
-        # At element center ξ = (0, 0, 0)
-        ξ = Vec(0.0, 0.0, 0.0)
-
-        # Hex8 basis function derivatives at center (new API)
-        dN_dξ = get_basis_derivatives(Hexahedron(), Lagrange{Hexahedron,1}(), ξ)
-
-        # Compute Jacobian
-        J = zero(Tensor{2,3,Float64,9})
-        for i in 1:8
-            J += X_nodes[i] ⊗ dN_dξ[i]
-        end
-
-        # Finite strain: Should give F = I + 0 = I
-        F_finite = compute_deformation_gradient(X_nodes, u_nodes, dN_dξ, J, FiniteStrain())
-        @test F_finite ≈ one(Tensor{2,3})
-        @test det(F_finite) ≈ 1.0
-
-        # Small strain: Should also give F = I
-        F_small = compute_deformation_gradient(X_nodes, u_nodes, dN_dξ, J, SmallStrain())
-        @test F_small ≈ one(Tensor{2,3})
-        @test det(F_small) ≈ 1.0
-    end
-
-    @testset "Pure translation" begin
-        # Unit cube
-        X_nodes = (
-            Vec(0.0, 0.0, 0.0),
-            Vec(1.0, 0.0, 0.0),
-            Vec(1.0, 1.0, 0.0),
-            Vec(0.0, 1.0, 0.0),
-            Vec(0.0, 0.0, 1.0),
-            Vec(1.0, 0.0, 1.0),
-            Vec(1.0, 1.0, 1.0),
-            Vec(0.0, 1.0, 1.0)
-        )
-
-        # Uniform translation: u = (0.5, 0.5, 0.5) everywhere
-        u_const = Vec(0.5, 0.5, 0.5)
-        u_nodes = tuple([u_const for _ in 1:8]...)
-
-        ξ = Vec(0.0, 0.0, 0.0)
-        dN_dξ = get_basis_derivatives(Hexahedron(), Lagrange{Hexahedron,1}(), ξ)
-
-        J = zero(Tensor{2,3,Float64,9})
-        for i in 1:8
-            J += X_nodes[i] ⊗ dN_dξ[i]
-        end
-
-        # Pure translation ⇒ ∇u = 0 ⇒ F = I
-        F = compute_deformation_gradient(X_nodes, u_nodes, dN_dξ, J, FiniteStrain())
-        @test F ≈ one(Tensor{2,3})
-        @test det(F) ≈ 1.0
-    end
-
-    @testset "Pure stretch in x-direction" begin
-        # Unit cube
-        X_nodes = (
-            Vec(0.0, 0.0, 0.0),
-            Vec(1.0, 0.0, 0.0),
-            Vec(1.0, 1.0, 0.0),
-            Vec(0.0, 1.0, 0.0),
-            Vec(0.0, 0.0, 1.0),
-            Vec(1.0, 0.0, 1.0),
-            Vec(1.0, 1.0, 1.0),
-            Vec(0.0, 1.0, 1.0)
-        )
-
-        # Stretch: u_x = 0.1 * X (10% stretch in x)
-        u_nodes = (
-            Vec(0.0, 0.0, 0.0),  # u = 0.1 * 0 = 0
-            Vec(0.1, 0.0, 0.0),  # u = 0.1 * 1 = 0.1
-            Vec(0.1, 0.0, 0.0),  # u = 0.1 * 1 = 0.1
-            Vec(0.0, 0.0, 0.0),  # u = 0.1 * 0 = 0
-            Vec(0.0, 0.0, 0.0),  # u = 0.1 * 0 = 0
-            Vec(0.1, 0.0, 0.0),  # u = 0.1 * 1 = 0.1
-            Vec(0.1, 0.0, 0.0),  # u = 0.1 * 1 = 0.1
-            Vec(0.0, 0.0, 0.0)   # u = 0.1 * 0 = 0
-        )
-
-        ξ = Vec(0.0, 0.0, 0.0)
-        dN_dξ = get_basis_derivatives(Hexahedron(), Lagrange{Hexahedron,1}(), ξ)
-
-        J = zero(Tensor{2,3,Float64,9})
-        for i in 1:8
-            J += X_nodes[i] ⊗ dN_dξ[i]
-        end
-
-        F = compute_deformation_gradient(X_nodes, u_nodes, dN_dξ, J, FiniteStrain())
-
-        # Expected: F = [1.1  0  0]
-        #               [0    1  0]
-        #               [0    0  1]
-        @test F[1, 1] ≈ 1.1 atol = 1e-10
-        @test F[2, 2] ≈ 1.0 atol = 1e-10
-        @test F[3, 3] ≈ 1.0 atol = 1e-10
-        @test F[1, 2] ≈ 0.0 atol = 1e-10
-        @test F[1, 3] ≈ 0.0 atol = 1e-10
-        @test F[2, 3] ≈ 0.0 atol = 1e-10
-        @test det(F) ≈ 1.1 atol = 1e-10
-    end
-
-    @testset "Simple shear" begin
-        # Unit cube
-        X_nodes = (
-            Vec(0.0, 0.0, 0.0),
-            Vec(1.0, 0.0, 0.0),
-            Vec(1.0, 1.0, 0.0),
-            Vec(0.0, 1.0, 0.0),
-            Vec(0.0, 0.0, 1.0),
-            Vec(1.0, 0.0, 1.0),
-            Vec(1.0, 1.0, 1.0),
-            Vec(0.0, 1.0, 1.0)
-        )
-
-        # Shear: u_x = 0.1 * y
-        u_nodes = (
-            Vec(0.0, 0.0, 0.0),   # y=0
-            Vec(0.0, 0.0, 0.0),   # y=0
-            Vec(0.1, 0.0, 0.0),   # y=1
-            Vec(0.1, 0.0, 0.0),   # y=1
-            Vec(0.0, 0.0, 0.0),   # y=0
-            Vec(0.0, 0.0, 0.0),   # y=0
-            Vec(0.1, 0.0, 0.0),   # y=1
-            Vec(0.1, 0.0, 0.0)    # y=1
-        )
-
-        ξ = Vec(0.0, 0.0, 0.0)
-        dN_dξ = get_basis_derivatives(Hexahedron(), Lagrange{Hexahedron,1}(), ξ)
-
-        J = zero(Tensor{2,3,Float64,9})
-        for i in 1:8
-            J += X_nodes[i] ⊗ dN_dξ[i]
-        end
-
-        F = compute_deformation_gradient(X_nodes, u_nodes, dN_dξ, J, FiniteStrain())
-
-        # Expected: F = [1    0.1  0]
-        #               [0    1    0]
-        #               [0    0    1]
-        @test F[1, 1] ≈ 1.0 atol = 1e-10
-        @test F[1, 2] ≈ 0.1 atol = 1e-10
-        @test F[2, 2] ≈ 1.0 atol = 1e-10
-        @test F[3, 3] ≈ 1.0 atol = 1e-10
-        @test det(F) ≈ 1.0 atol = 1e-10
-    end
-
-    @testset "Small vs Finite strain difference" begin
-        # Setup with significant displacement gradient
-        X_nodes = (
-            Vec(0.0, 0.0, 0.0),
-            Vec(1.0, 0.0, 0.0),
-            Vec(1.0, 1.0, 0.0),
-            Vec(0.0, 1.0, 0.0),
-            Vec(0.0, 0.0, 1.0),
-            Vec(1.0, 0.0, 1.0),
-            Vec(1.0, 1.0, 1.0),
-            Vec(0.0, 1.0, 1.0)
-        )
-
-        # 20% stretch in x
-        u_nodes = (
-            Vec(0.0, 0.0, 0.0),
-            Vec(0.2, 0.0, 0.0),
-            Vec(0.2, 0.0, 0.0),
-            Vec(0.0, 0.0, 0.0),
-            Vec(0.0, 0.0, 0.0),
-            Vec(0.2, 0.0, 0.0),
-            Vec(0.2, 0.0, 0.0),
-            Vec(0.0, 0.0, 0.0)
-        )
-
-        ξ = Vec(0.0, 0.0, 0.0)
-        dN_dξ = get_basis_derivatives(Hexahedron(), Lagrange{Hexahedron,1}(), ξ)
-
-        J = zero(Tensor{2,3,Float64,9})
-        for i in 1:8
-            J += X_nodes[i] ⊗ dN_dξ[i]
-        end
-
-        F_finite = compute_deformation_gradient(X_nodes, u_nodes, dN_dξ, J, FiniteStrain())
-        F_small = compute_deformation_gradient(X_nodes, u_nodes, dN_dξ, J, SmallStrain())
-
-        # Finite strain includes gradient
-        @test F_finite[1, 1] ≈ 1.2 atol = 1e-10
-
-        # Small strain ignores gradient
-        @test F_small[1, 1] ≈ 1.0 atol = 1e-10
-
-        # They should be different!
-        @test !(F_finite ≈ F_small)
-    end
-
-    @testset "Physical constraint: det(F) > 0" begin
-        # Physical deformation must preserve orientation
-        X_nodes = (
-            Vec(0.0, 0.0, 0.0),
-            Vec(1.0, 0.0, 0.0),
-            Vec(1.0, 1.0, 0.0),
-            Vec(0.0, 1.0, 0.0),
-            Vec(0.0, 0.0, 1.0),
-            Vec(1.0, 0.0, 1.0),
-            Vec(1.0, 1.0, 1.0),
-            Vec(0.0, 1.0, 1.0)
-        )
-
-        # Small positive stretch
-        u_nodes = (
-            Vec(0.0, 0.0, 0.0),
-            Vec(0.05, 0.0, 0.0),
-            Vec(0.05, 0.0, 0.0),
-            Vec(0.0, 0.0, 0.0),
-            Vec(0.0, 0.0, 0.0),
-            Vec(0.05, 0.0, 0.0),
-            Vec(0.05, 0.0, 0.0),
-            Vec(0.0, 0.0, 0.0)
-        )
-
-        ξ = Vec(0.0, 0.0, 0.0)
-        dN_dξ = get_basis_derivatives(Hexahedron(), Lagrange{Hexahedron,1}(), ξ)
-
-        J = zero(Tensor{2,3,Float64,9})
-        for i in 1:8
-            J += X_nodes[i] ⊗ dN_dξ[i]
-        end
-
-        F = compute_deformation_gradient(X_nodes, u_nodes, dN_dξ, J, FiniteStrain())
-
-        @test det(F) > 0  # Physical requirement
-    end
-end
-
-@testset "Deformation Gradient - Tet10 Element" begin
-
-    @testset "Tet10: Identity case" begin
-        # Regular tetrahedron nodes (4 corners + 6 edge midpoints)
-        X_nodes = (
-            Vec(0.0, 0.0, 0.0),              # 1: corner
-            Vec(1.0, 0.0, 0.0),              # 2: corner
-            Vec(0.0, 1.0, 0.0),              # 3: corner
-            Vec(0.0, 0.0, 1.0),              # 4: corner
-            Vec(0.5, 0.0, 0.0),              # 5: edge 1-2
-            Vec(0.5, 0.5, 0.0),              # 6: edge 2-3
-            Vec(0.0, 0.5, 0.0),              # 7: edge 3-1
-            Vec(0.0, 0.0, 0.5),              # 8: edge 1-4
-            Vec(0.5, 0.0, 0.5),              # 9: edge 2-4
-            Vec(0.0, 0.5, 0.5)               # 10: edge 3-4
-        )
-
-        # Zero displacement
-        u_nodes = tuple([zero(Vec{3,Float64}) for _ in 1:10]...)
-
-        # At element centroid ξ = (1/4, 1/4, 1/4)
-        ξ = Vec(0.25, 0.25, 0.25)
-
-        # Tet10 basis function derivatives (new API)
-        dN_dξ = get_basis_derivatives(Tetrahedron(), Lagrange{Tetrahedron,2}(), ξ)
-
-        # Compute Jacobian
-        J = zero(Tensor{2,3,Float64,9})
-        for i in 1:10
-            J += X_nodes[i] ⊗ dN_dξ[i]
-        end
-
-        F = compute_deformation_gradient(X_nodes, u_nodes, dN_dξ, J, FiniteStrain())
-
-        @test F ≈ one(Tensor{2,3}) atol = 1e-10
-        @test det(F) ≈ 1.0 atol = 1e-10
-    end
-
-    @testset "Tet10: Uniform stretch" begin
-        # Regular tetrahedron
-        X_nodes = (
-            Vec(0.0, 0.0, 0.0),
-            Vec(1.0, 0.0, 0.0),
-            Vec(0.0, 1.0, 0.0),
-            Vec(0.0, 0.0, 1.0),
-            Vec(0.5, 0.0, 0.0),
-            Vec(0.5, 0.5, 0.0),
-            Vec(0.0, 0.5, 0.0),
-            Vec(0.0, 0.0, 0.5),
-            Vec(0.5, 0.0, 0.5),
-            Vec(0.0, 0.5, 0.5)
-        )
-
-        # Isotropic expansion: u = 0.1 * X
-        u_nodes = (
-            Vec(0.0, 0.0, 0.0),
-            Vec(0.1, 0.0, 0.0),
-            Vec(0.0, 0.1, 0.0),
-            Vec(0.0, 0.0, 0.1),
-            Vec(0.05, 0.0, 0.0),
-            Vec(0.05, 0.05, 0.0),
-            Vec(0.0, 0.05, 0.0),
-            Vec(0.0, 0.0, 0.05),
-            Vec(0.05, 0.0, 0.05),
-            Vec(0.0, 0.05, 0.05)
-        )
-
-        ξ = Vec(0.25, 0.25, 0.25)
-        dN_dξ = get_basis_derivatives(Tetrahedron(), Lagrange{Tetrahedron,2}(), ξ)
-
-        J = zero(Tensor{2,3,Float64,9})
-        for i in 1:10
-            J += X_nodes[i] ⊗ dN_dξ[i]
-        end
-
-        F = compute_deformation_gradient(X_nodes, u_nodes, dN_dξ, J, FiniteStrain())
-
-        # Expected: F ≈ 1.1 * I
-        @test F[1, 1] ≈ 1.1 atol = 1e-10
-        @test F[2, 2] ≈ 1.1 atol = 1e-10
-        @test F[3, 3] ≈ 1.1 atol = 1e-10
-        @test abs(F[1, 2]) < 1e-10
-        @test abs(F[1, 3]) < 1e-10
-        @test abs(F[2, 3]) < 1e-10
-        @test det(F) ≈ 1.1^3 atol = 1e-10
-    end
-end
-
-@testset "Deformation Gradient - Zero Allocation" begin
-
-    @testset "Verify zero allocations" begin
-        # Setup
-        X_nodes = (
-            Vec(0.0, 0.0, 0.0),
-            Vec(1.0, 0.0, 0.0),
-            Vec(1.0, 1.0, 0.0),
-            Vec(0.0, 1.0, 0.0),
-            Vec(0.0, 0.0, 1.0),
-            Vec(1.0, 0.0, 1.0),
-            Vec(1.0, 1.0, 1.0),
-            Vec(0.0, 1.0, 1.0)
-        )
-
-        u_nodes = (
-            Vec(0.0, 0.0, 0.0),
-            Vec(0.1, 0.0, 0.0),
-            Vec(0.1, 0.0, 0.0),
-            Vec(0.0, 0.0, 0.0),
-            Vec(0.0, 0.0, 0.0),
-            Vec(0.1, 0.0, 0.0),
-            Vec(0.1, 0.0, 0.0),
-            Vec(0.0, 0.0, 0.0)
-        )
-
-        ξ = Vec(0.0, 0.0, 0.0)
-        dN_dξ = get_basis_derivatives(Hexahedron(), Lagrange{Hexahedron,1}(), ξ)
-
-        J = zero(Tensor{2,3,Float64,9})
-        for i in 1:8
-            J += X_nodes[i] ⊗ dN_dξ[i]
-        end
-
-        # Warm up (compile)
-        F = compute_deformation_gradient(X_nodes, u_nodes, dN_dξ, J, FiniteStrain())
-
-        # Measure allocations
-        allocs = @allocated compute_deformation_gradient(X_nodes, u_nodes, dN_dξ, J, FiniteStrain())
-
-        @test allocs == 0  # Zero allocations!
-    end
-end