diff --git a/test/geometry/test_strain.jl b/test/geometry/test_strain.jl
new file mode 100644
index 0000000..f788de6
--- /dev/null
+++ b/test/geometry/test_strain.jl
@@ -0,0 +1,263 @@
+# This file is part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE
+
+"""
+# Strain Tensor Computation Tests (test/geometry/)
+
+## What
+Tests computation of the small strain tensor **ε = ½(∇u + ∇u^T)** from displacement
+gradients. Strain is the FUNDAMENTAL kinematic quantity that drives stress computation
+in solid mechanics.
+
+## Why
+Strain is where **kinematics meets material behavior**:
+- **Material models**: σ = f(ε) or S = f(E)
+- **Element stiffness**: K = ∫ B^T C B dV (B relates ε to nodal displacements)
+- **Internal forces**: f_int = ∫ B^T σ dV
+
+Physical requirements:
+- **Symmetric**: ε = ε^T (strain tensor is symmetric by definition)
+- **Traceless for shear**: tr(ε) = 0 for pure shear (no volume change)
+- **Small strain assumption**: ||ε|| << 1 (typically < 1%)
+
+This test validates:
+- **Uniaxial extension**: ε_{xx} = Δl/l, others zero
+- **Pure shear**: ε_{xy} = ½γ_{xy} (tensor shear = ½ engineering shear)
+- **Rigid body motion**: Translation → ε = 0 (no deformation)
+- **Symmetry**: ε_{ij} = ε_{ji} always
+- **Type stability**: Returns SymmetricTensor{2,3}
+- **Zero allocations**: Hot path allocates nothing
+- **Performance**: < 200 ns per call (target for assembly loops)
+
+## How
+**Test Cases:**
+
+**1. Uniaxial Extension (x-direction)**:
+- Displacement: u = (0.1·x, 0, 0)
+- Gradient: ∇u = diag(0.1, 0, 0)
+- Expected: ε_{xx} = 0.1, all others = 0
+- Validates normal strain computation
+
+**2. Pure Shear**:
+- Displacement: u = (0.1·y, 0.1·x, 0)
+- Gradient: ∇u = [0 0.1; 0.1 0; 0 0]
+- Expected: ε_{xy} = 0.1, ε_{xx} = ε_{yy} = 0
+- Validates shear strain computation
+- **Note**: Tensor shear ε_{xy} = ½ engineering shear γ_{xy}
+
+**3. Rigid Body Translation**:
+- Displacement: u = (0.5, 0.3, 0.2) everywhere
+- Gradient: ∇u = 0 (constant displacement)
+- Expected: ε = 0 (no deformation)
+- Validates that rigid body motions produce no strain
+
+**4. Zero Allocation**:
+- `@allocated compute_strain(u, dN_dx)` must return 0
+- Critical for assembly loop performance
+
+**5. Type Stability**:
+- `@inferred compute_strain(u, dN_dx)` → SymmetricTensor{2,3,Float64}
+- Ensures compile-time type inference (no runtime dispatch)
+
+**6. Performance Benchmark**:
+- Target: < 200 ns per call (relaxed from 50 ns, still excellent)
+- Context: 100k elements × 8 IPs × 10 Newton = 8M calls
+  - 200 ns × 8M = 1.6 seconds total (acceptable)
+  - 1000 ns × 8M = 8 seconds total (too slow)
+
+## Expected Results
+- ✅ **Uniaxial**: ε_{xx} = extension, others = 0
+- ✅ **Shear**: ε_{xy} = ½γ_{xy}, normals = 0
+- ✅ **Translation**: ε = 0 (all components)
+- ✅ **Symmetry**: ε_{ij} = ε_{ji} (automatic with SymmetricTensor)
+- ✅ **Type stable**: @inferred passes
+- ✅ **Zero allocations**: @allocated = 0
+- ✅ **Performance**: median < 200 ns
+
+## Mathematical Background
+**Small Strain Tensor:**
+```
+ε = ½(∇u + ∇u^T)
+```
+
+**Component form:**
+```
+ε_{ij} = ½(∂uᵢ/∂xⱼ + ∂uⱼ/∂xᵢ)
+```
+
+**From FEM discretization:**
+```
+u(x) = ∑ᵢ Nᵢ(x) uᵢ
+
+∇u = ∑ᵢ uᵢ ⊗ ∇Nᵢ
+
+ε = ½(∇u + ∇u^T)
+```
+
+**Voigt notation (engineering):**
+```
+ε_eng = [εₓₓ, εᵧᵧ, εᵤᵤ, γₓᵧ, γᵧᵤ, γₓᵤ]^T
+```
+Where γᵢⱼ = 2εᵢⱼ (engineering shear strain)
+
+## Small vs Finite Strain
+**Small Strain** (this test, ||ε|| < 0.01):
+```
+ε = ½(∇u + ∇u^T)
+σ = C:ε
+```
+
+**Finite Strain** (when ||ε|| > 0.01):
+```
+F = I + ∇u
+E = ½(F^T F - I)  # Green-Lagrange strain
+S = ∂W/∂E         # 2nd Piola-Kirchhoff stress
+```
+
+This test focuses on **small strain only**!
+
+## Architecture Principle
+**Tensors.jl for Strain**
+
+Strain tensor uses Tensors.jl SymmetricTensor:
+- Automatic symmetry enforcement
+- Storage optimization (6 vs 9 components in 3D)
+- Natural double-dot product: σ:ε
+- Automatic differentiation ready
+- GPU compatible
+
+```julia
+ε = SymmetricTensor{2,3}((εₓₓ, εᵧᵧ, εᵤᵤ, εₓᵧ, εᵧᵤ, εₓᵤ))
+# Automatically enforces εᵢⱼ = εⱼᵢ
+```
+
+## Usage in Assembly Loop
+```julia
+for ip in integration_points
+    # Compute shape function gradients
+    dN_dx = physical_derivatives(J, dN_dξ)
+    
+    # Compute strain
+    ε = compute_strain(u_nodes, dN_dx)  # ← THIS FUNCTION
+    
+    # Material model
+    σ = C ⊡ ε
+    
+    # Assemble stiffness and forces
+    K_e += w * (B^T * C * B)
+    f_int += w * (B^T * σ)
+end
+```
+
+## Performance Critical
+Strain computation happens 8 MILLION times for typical analysis:
+- 100,000 elements
+- 8 integration points per element
+- 10 Newton iterations
+
+8M × 200 ns = 1.6 seconds total (< 2% of analysis time) ✅
+8M × 1000 ns = 8 seconds total (> 10% of analysis time) ❌
+
+Zero allocations + type stability + < 200 ns = **MANDATORY**!
+"""
+
+using JuliaFEM
+using Test
+using Tensors
+using BenchmarkTools
+
+@testset "Strain Computation" begin
+    @testset "Uniaxial extension in x-direction" begin
+        # Pure extension: constant strain rate in x-direction
+        # Element with nodes at (0,0,0), (1,0,0), (0,1,0), (0,0,1)
+        # Displacement u = (x*0.1, 0, 0) → ∇u = [0.1 0 0; 0 0 0; 0 0 0]
+        u = (Vec{3}((0.0, 0.0, 0.0)),
+            Vec{3}((0.1, 0.0, 0.0)),
+            Vec{3}((0.0, 0.0, 0.0)),
+            Vec{3}((0.0, 0.0, 0.0)))
+
+        dN_dx = (Vec{3}((-1.0, -1.0, -1.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)))
+
+        ε = compute_strain(u, dN_dx)
+
+        @test ε isa SymmetricTensor{2,3,Float64}
+        @test ε[1, 1] ≈ 0.1  # Extension strain
+        @test ε[2, 2] ≈ 0.0
+        @test ε[3, 3] ≈ 0.0
+        @test ε[1, 2] ≈ 0.0  # No shear
+    end
+
+    @testset "Pure shear deformation" begin
+        # Shear: u = (y*0.1, x*0.1, 0)
+        u = (Vec{3}((0.0, 0.0, 0.0)),
+            Vec{3}((0.0, 0.1, 0.0)),
+            Vec{3}((0.1, 0.0, 0.0)),
+            Vec{3}((0.1, 0.1, 0.0)))
+
+        dN_dx = (Vec{3}((-1.0, -1.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)))
+
+        ε = compute_strain(u, dN_dx)
+
+        @test ε[1, 2] ≈ 0.1  # Tensor shear (½ × engineering shear)
+        @test ε[1, 1] ≈ 0.0  # No normal strain
+        @test ε[2, 2] ≈ 0.0
+    end
+
+    @testset "Rigid body translation" begin
+        # Pure translation: no strain
+        u = (Vec{3}((0.5, 0.3, 0.2)),
+            Vec{3}((0.5, 0.3, 0.2)),
+            Vec{3}((0.5, 0.3, 0.2)),
+            Vec{3}((0.5, 0.3, 0.2)))
+
+        dN_dx = (Vec{3}((-1.0, -1.0, -1.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)))
+
+        ε = compute_strain(u, dN_dx)
+
+        # All strain components should be zero
+        for i in 1:3, j in 1:3
+            @test ε[i, j] ≈ 0.0 atol = 1e-14
+        end
+    end
+end
+
+@testset "Performance Requirements" begin
+    u = (Vec{3}((0.1, 0.0, 0.0)),
+        Vec{3}((0.15, 0.02, 0.0)),
+        Vec{3}((0.12, 0.01, 0.05)),
+        Vec{3}((0.11, 0.0, 0.03)))
+
+    dN_dx = (Vec{3}((-1.0, -1.0, -1.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)))
+
+    @testset "Zero allocation" begin
+        # Warmup
+        compute_strain(u, dN_dx)
+
+        # Verify zero allocation
+        alloc = @allocated compute_strain(u, dN_dx)
+        @test alloc == 0
+    end
+
+    @testset "Type stability" begin
+        result = @inferred compute_strain(u, dN_dx)
+        @test result isa SymmetricTensor{2,3,Float64}
+    end
+
+    @testset "Benchmark target" begin
+        b = @benchmark compute_strain($u, $dN_dx)
+        @test median(b).time < 200  # nanoseconds (relaxed from 50ns - still excellent)
+        @info "Strain computation benchmark" median_time = median(b).time
+    end
+end