diff -r 529b533a1dbb -r b30db2ea34ed src/SbpOperators/boundaryops/normal_derivative.jl
--- a/src/SbpOperators/boundaryops/normal_derivative.jl	Sun Jun 30 10:57:05 2024 +0200
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Sun Jun 30 15:57:22 2024 +0200
@@ -42,13 +42,19 @@
         gᵏⁱ./sqrt(gᵏᵏ)
     end
 
+    σ = ScalingTensor(
+        Grids._boundary_sign(component_type(g), boundary),
+        size(boundary_grid(g,boundary)),
+    )
+
+
     # Assemble difference operator
     mapreduce(+,1:ndims(g)) do i
         if i == k
             ∂_ξᵢ = normal_derivative(logicalgrid(g), stencil_set, boundary)
         else
             e = boundary_restriction(logicalgrid(g), stencil_set, boundary)
-            ∂_ξᵢ = e ∘ first_derivative(logicalgrid(g), stencil_set, i)
+            ∂_ξᵢ = σ ∘ e ∘ first_derivative(logicalgrid(g), stencil_set, i)
         end
 
         αᵢ = componentview(α,i)

diff -r 529b533a1dbb -r b30db2ea34ed test/SbpOperators/boundaryops/normal_derivative_test.jl
--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl	Sun Jun 30 10:57:05 2024 +0200
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Sun Jun 30 15:57:22 2024 +0200
@@ -66,15 +66,44 @@
             @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)
+
+        # x̄((ξ, η)) = @SVector[ξ, η*(1+ξ*(ξ-1))]
+        # J((ξ, η)) = @SMatrix[
+        #     1         0;
+        #     η*(2ξ-1)  1+ξ*(ξ-1);
+        # ]
+        # mg = mapped_grid(x̄, J, 20, 21)
+
+
+        # x̄((ξ, η)) = @SVector[ξ,η]
+        # J((ξ, η)) = @SMatrix[
+        #     1  0;
+        #     0  1;
+        # ]
+        # mg = mapped_grid(identity, J, 10, 11)
 
-        @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}
+        for bid ∈ boundary_identifiers(mg)
+            @testset let bid=bid
+                @test normal_derivative(mg, stencil_set, bid) isa LazyTensor{<:Any, 1, 2}
+            end
+        end
+
+        @testset "Consistency" begin
+            v = map(identity, mg)
+
+             @testset "4nd order" begin
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+
+                for bid ∈ boundary_identifiers(mg)
+                    @testset let bid=bid
+                        d = normal_derivative(mg, stencil_set, bid)
+                        @test d*v ≈ normal(mg, bid) rtol=1e-13
+                    end
+                end
+             end
+        end
 
         @testset "Accuracy" begin
             v = map(x̄ -> NaN, mg)