--- a/src/SbpOperators/boundaryops/normal_derivative.jl	Wed Jun 26 13:42:19 2024 +0200
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Wed Jun 26 14:41:50 2024 +0200
@@ -28,3 +28,14 @@
     scaled_stencil = scale(closure_stencil,h_inv)
     return BoundaryOperator(g, scaled_stencil, boundary)
 end
+
+function normal_derivative(g::MappedGrid, stencil_set::StencilSet, boundary)
+    g⁻¹ = geometric_tensor_inverse(g) # Extract boundary part
+    k = NaN # Dimension of boundary
+    mapreduce(1:ndims(g)) do i
+        gᵏⁱ = componentview(g⁻¹,k,i)
+        gᵏᵏ = componentview(g⁻¹,k,k)
+        # ∂ξᵢ = ...
+        DiagonalTensor(gᵏⁱ./sqrt.(gᵏᵏ)) * ∂ξᵢ # Should the metric expression be mapped lazily?
+    end
+end

--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl	Wed Jun 26 13:42:19 2024 +0200
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Wed Jun 26 14:41:50 2024 +0200
@@ -6,6 +6,8 @@
 using Sbplib.RegionIndices
 import Sbplib.SbpOperators.BoundaryOperator
 
+using StaticArrays
+
 @testset "normal_derivative" begin
 	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
 
@@ -58,4 +60,45 @@
             end
         end
     end
+
+    @testset "MappedGrid" begin
+        c = Chart(unitsquare()) do (ξ,η)
+            @SVector[2ξ + η*(1-η), 3η+(1+η/2)*ξ^2]
+        end
+        Grids.jacobian(c::typeof(c), (ξ,η)) = @SMatrix[2 1-2η; (2+η)*ξ 3+ξ^2/2]
+
+        mg = equidistant_grid(c, 10,13)
+
+        b_w, b_e, b_s, b_n = boundary_identifiers(mg)
+
+        @test_broken normal_derivative(mg, stencil_set, b_w) isa LazyTensor{<:Any, 1, 2}
+        @test_broken normal_derivative(mg, stencil_set, b_e) isa LazyTensor{<:Any, 1, 2}
+        @test_broken normal_derivative(mg, stencil_set, b_s) isa LazyTensor{<:Any, 1, 2}
+        @test_broken normal_derivative(mg, stencil_set, b_n) isa LazyTensor{<:Any, 1, 2}
+
+        @testset "Accuracy" begin
+            v = map(x̄ -> NaN, mg)
+            dₙv = map(x̄ -> NaN, mg)
+
+            @testset "2nd order" begin
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+                d_w, d_e, d_s, d_n = normal_derivative.(Ref(mg), Ref(stencil_set), boundary_identifiers(mg))
+
+                @test_broken d_w*v ≈ dₙv atol = 1e-13
+                @test_broken d_e*v ≈ dₙv atol = 1e-13
+                @test_broken d_s*v ≈ dₙv atol = 1e-13
+                @test_broken d_n*v ≈ dₙv atol = 1e-13
+            end
+
+            @testset "4th order" begin
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+                d_w, d_e, d_s, d_n = normal_derivative.(Ref(mg), Ref(stencil_set), boundary_identifiers(mg))
+
+                @test_broken d_w*v ≈ dₙv atol = 1e-13
+                @test_broken d_e*v ≈ dₙv atol = 1e-13
+                @test_broken d_s*v ≈ dₙv atol = 1e-13
+                @test_broken d_n*v ≈ dₙv atol = 1e-13
+            end
+        end
+    end
 end