diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index ca28e2f..a890604 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -190,6 +190,9 @@ export get_gauss_points!  # NEW: Zero-allocation integration point API
 include("geometry/jacobian.jl")
 export compute_jacobian, physical_derivatives
 
+include("geometry/strain.jl")
+export compute_strain
+
 # ============================================================================
 # BASIS: Interpolation schemes (consolidated from FEMBasis.jl)
 # ============================================================================
diff --git a/src/geometry/strain.jl b/src/geometry/strain.jl
new file mode 100644
index 0000000..3389d3c
--- /dev/null
+++ b/src/geometry/strain.jl
@@ -0,0 +1,45 @@
+# This file is part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE
+
+"""
+Strain computation utilities for finite element analysis.
+
+This module provides functions for computing strain tensors from nodal
+displacements and basis function derivatives.
+"""
+
+using Tensors
+
+"""
+    compute_strain(u_elem, dN_dx) -> SymmetricTensor{2,3}
+
+Compute small strain tensor from nodal displacements.
+
+Small strain: ε = ½(∇u + ∇u^T) where ∇u = ∑ᵢ uᵢ ⊗ (dNᵢ/dx)
+
+# Arguments
+- `u_elem::NTuple{N, Vec{3}}`: Nodal displacement vectors
+- `dN_dx::NTuple{N, Vec{3}}`: Basis derivatives in physical coordinates
+
+# Returns
+- `ε::SymmetricTensor{2,3}`: Small strain tensor
+
+# Example
+```julia
+u = (Vec{3}((0.1, 0.0, 0.0)), Vec{3}((0.15, 0.02, 0.0)))
+dN_dx = (Vec{3}((-1.0, -1.0, -1.0)), Vec{3}((1.0, 0.0, 0.0)))
+ε = compute_strain(u, dN_dx)
+```
+"""
+function compute_strain(
+    u_elem::NTuple{N,Vec{3}},
+    dN_dx::NTuple{N,Vec{3}}
+) where N
+    # Displacement gradient: ∇u = ∑ᵢ uᵢ ⊗ (dNᵢ/dx)
+    ∇u = sum(u_i ⊗ dN_i for (u_i, dN_i) in zip(u_elem, dN_dx))
+
+    # Small strain (symmetric part): ε = ½(∇u + ∇u^T)
+    ε = symmetric(∇u)
+
+    return ε
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index eb7a173..15bfb27 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,6 +7,9 @@ using JuliaFEM, Test
     @testset "test_dirichlet.jl" begin
         include("test_dirichlet.jl")
     end
+    @testset "test_strain.jl" begin
+        include("test_strain.jl")
+    end
     @testset "test_elasticity_1d.jl" begin
         include("test_elasticity_1d.jl")
     end
diff --git a/test/test_strain.jl b/test/test_strain.jl
new file mode 100644
index 0000000..c55ca27
--- /dev/null
+++ b/test/test_strain.jl
@@ -0,0 +1,103 @@
+# This file is part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE
+
+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