diff --git a/test/basis/test_gradient_properties.jl b/test/basis/test_gradient_properties.jl
new file mode 100644
index 0000000..8289cc1
--- /dev/null
+++ b/test/basis/test_gradient_properties.jl
@@ -0,0 +1,90 @@
+# Test gradient properties for all elements
+using Test
+using StaticArrays, Tensors
+
+include("../../src/topology/api.jl")
+include("../../src/topology/segments.jl")
+include("../../src/topology/triangles.jl")
+include("../../src/topology/quadrilaterals.jl")
+include("../../src/topology/tetrahedra.jl")
+include("../../src/topology/hexahedra.jl")
+include("../../src/topology/pyramids.jl")
+include("../../src/topology/wedges.jl")
+include("../../src/basis/api.jl")
+include("../../src/basis/basis_generated.jl")
+
+@testset "Gradient Properties" begin
+    @testset "Linear elements - constant gradients" begin
+        # Tri3 gradients should be constant everywhere
+        dN1 = get_basis_derivatives(Tri3(), Lagrange{1}(), Vec{2}((0.2, 0.3)))
+        dN2 = get_basis_derivatives(Tri3(), Lagrange{1}(), Vec{2}((0.5, 0.4)))
+
+        for i in 1:3
+            @test dN1[i] ≈ dN2[i] atol = 1e-14
+        end
+
+        # Tet4 gradients should be constant everywhere
+        dN1 = get_basis_derivatives(Tet4(), Lagrange{1}(), Vec{3}((0.1, 0.2, 0.3)))
+        dN2 = get_basis_derivatives(Tet4(), Lagrange{1}(), Vec{3}((0.2, 0.1, 0.4)))
+
+        for i in 1:4
+            @test dN1[i] ≈ dN2[i] atol = 1e-14
+        end
+    end
+
+    @testset "Sum of gradients = 0 (rigid body motion)" begin
+        # For any basis, sum of all gradients should be zero vector
+        # (ensures constant field has zero gradient)
+
+        # 1D
+        dN = get_basis_derivatives(Seg2(), Lagrange{1}(), Vec{1}((0.3,)))
+        @test sum(dN)[1] ≈ 0.0 atol = 1e-14
+
+        # 2D
+        dN = get_basis_derivatives(Tri6(), Lagrange{2}(), Vec{2}((0.3, 0.4)))
+        sum_grad = sum(dN)
+        @test sum_grad[1] ≈ 0.0 atol = 1e-14
+        @test sum_grad[2] ≈ 0.0 atol = 1e-14
+
+        dN = get_basis_derivatives(Quad9(), Lagrange{2}(), Vec{2}((0.5, -0.5)))
+        sum_grad = sum(dN)
+        @test sum_grad[1] ≈ 0.0 atol = 1e-14
+        @test sum_grad[2] ≈ 0.0 atol = 1e-14
+
+        # 3D
+        dN = get_basis_derivatives(Hex27(), Lagrange{2}(), Vec{3}((0.2, 0.3, -0.1)))
+        sum_grad = sum(dN)
+        @test sum_grad[1] ≈ 0.0 atol = 1e-14
+        @test sum_grad[2] ≈ 0.0 atol = 1e-14
+        @test sum_grad[3] ≈ 0.0 atol = 1e-14
+    end
+
+    @testset "Gradient return types" begin
+        # Test that gradients have correct Vec type
+        dN = get_basis_derivatives(Tri3(), Lagrange{1}(), Vec{2}((0.5, 0.25)))
+        @test dN isa SVector{3,Vec{2,Float64}}
+        @test dN[1] isa Vec{2,Float64}
+
+        dN = get_basis_derivatives(Hex8(), Lagrange{1}(), Vec{3}((0.0, 0.0, 0.0)))
+        @test dN isa SVector{8,Vec{3,Float64}}
+        @test dN[1] isa Vec{3,Float64}
+    end
+
+    @testset "Numerical gradient check (finite differences)" begin
+        # Simple finite difference check for Quad4
+        h = 1e-8
+        xi = Vec{2}((0.3, 0.4))
+
+        # Analytical gradient
+        dN = get_basis_derivatives(Quad4(), Lagrange{1}(), xi)
+
+        # Finite difference in u-direction
+        N_plus = get_basis_functions(Quad4(), Lagrange{1}(), Vec{2}((xi[1] + h, xi[2])))
+        N_minus = get_basis_functions(Quad4(), Lagrange{1}(), Vec{2}((xi[1] - h, xi[2])))
+        dN_u_fd = (N_plus - N_minus) / (2h)
+
+        for i in 1:4
+            @test dN[i][1] ≈ dN_u_fd[i] atol = 1e-6
+        end
+    end
+end