--- a/src/Grids/geometry.jl	Tue Feb 25 22:48:42 2025 +0100
+++ b/src/Grids/geometry.jl	Tue Feb 25 23:33:56 2025 +0100
@@ -118,5 +118,31 @@
     return true
 end
 
+function Grids.jacobian(s::TransfiniteInterpolationSurface, ξ̄)
+    u, v = ξ̄
+
+    c₁, c₂, c₃, c₄ = s.c₁, s.c₂, s.c₃, s.c₄
+    P₀₀ = c₁(0)
+    P₁₀ = c₂(0)
+    P₁₁ = c₃(0)
+    P₀₁ = c₄(0)
+
+    ∂x̄∂ξ₁ = (1-v)*jacobian(c₁,u) + c₂(v) - v*jacobian(c₃,1-u) -c₄(1-v) - (
+        -(1-v)*P₀₀ + (1-v)*P₁₀ + v*P₁₁ - v*P₀₁
+    )
+
+
+    (1-v)*c₁(u) + u*c₂(v) + v*c₃(1-u) + (1-u)*c₄(1-v) - (
+        (1-u)*(1-v)*P₀₀ + u*(1-v)*P₁₀ + u*v*P₁₁ + (1-u)*v*P₀₁
+    )
+
+
+    ∂x̄∂ξ₂ = -c₁(u) + u*jacobian(c₂,v) + c₃(1-u) - (1-u)*jacobian(c₄,1-v) - (
+        -(1-u)*P₀₀ - u*P₁₀ + u*P₁₁ + (1-u)*P₀₁
+    )
+
+    return [∂x̄∂ξ₁ ∂x̄∂ξ₂]
+end
+
 # TODO: Implement jacobian() for the different mapping helpers
 # TODO: Add doc strings

--- a/test/Grids/geometry_test.jl	Tue Feb 25 22:48:42 2025 +0100
+++ b/test/Grids/geometry_test.jl	Tue Feb 25 23:33:56 2025 +0100
@@ -171,7 +171,6 @@
     end
 
     @testset "Evaluation" begin
-        a, b, c, d = [1,0],[2,1/4],[2.5,1],[-1/3,1]
         cs = polygon_edges([0,0],[1,0],[1,1],[0,1])
         ti = TransfiniteInterpolationSurface(cs...)
 
@@ -236,4 +235,40 @@
         @test_throws Exception check_transfiniteinterpolation(s3)
         @test_throws Exception check_transfiniteinterpolation(s4)
     end
+
+    @testset "Grids.jacobian" begin
+        cs = polygon_edges([0,0],[1,0],[1,1],[0,1])
+        ti = TransfiniteInterpolationSurface(cs...)
+
+        @test Grids.jacobian(ti, [0,0]) isa SMatrix
+
+        @test Grids.jacobian(ti, [0,0]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [1,0]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [1,1]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [0,1]) == [1 0; 0 1]
+
+        @test Grids.jacobian(ti, [1/2, 0]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [1/2, 1]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [0,1/2]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [1,1/2]) == [1 0; 0 1]
+
+
+        a, b, c, d = [1,0],[2,1/4],[2.5,1],[-1/3,1]
+        cs = polygon_edges(a,b,c,d)
+        ti = TransfiniteInterpolationSurface(cs...)
+
+        @test Grids.jacobian(ti, [0,0]) == [b-a d-a]
+        @test Grids.jacobian(ti, [1,0]) == [b-a c-b]
+        @test Grids.jacobian(ti, [1,1]) == [c-d c-b]
+        @test Grids.jacobian(ti, [0,1]) == [c-d d-a]
+
+
+        mid(x,y) = (x+y)/2
+        @test Grids.jacobian(ti, [1/2, 0]) ≈ [b-a mid(c,d)-mid(a,b)]
+        @test Grids.jacobian(ti, [1/2, 1]) ≈ [c-d mid(c,d)-mid(a,b)]
+        @test Grids.jacobian(ti, [0, 1/2]) ≈ [mid(b,c)-mid(a,d) d-a]
+        @test Grids.jacobian(ti, [1, 1/2]) ≈ [mid(b,c)-mid(a,d) c-b]
+
+        # TODO: Some test with curved edges?
+    end
 end