--- a/src/SbpOperators/SbpOperators.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/SbpOperators/SbpOperators.jl	Wed Aug 28 10:50:15 2024 +0200
@@ -43,6 +43,7 @@
 using Sbplib.RegionIndices
 using Sbplib.LazyTensors
 using Sbplib.Grids
+using LinearAlgebra
 
 # Includes
 include("stencil.jl")

--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Wed Aug 28 10:50:15 2024 +0200
@@ -25,3 +25,7 @@
     converted_stencil = convert(Stencil{eltype(g)}, closure_stencil)
     return BoundaryOperator(g, converted_stencil, boundary)
 end
+
+function boundary_restriction(g::MappedGrid, stencil_set::StencilSet, boundary)
+    return boundary_restriction(logicalgrid(g), stencil_set, boundary)
+end

--- a/src/SbpOperators/boundaryops/normal_derivative.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Wed Aug 28 10:50:15 2024 +0200
@@ -28,3 +28,36 @@
     scaled_stencil = scale(closure_stencil,h_inv)
     return BoundaryOperator(g, scaled_stencil, boundary)
 end
+
+function normal_derivative(g::MappedGrid, stencil_set::StencilSet, boundary)
+    k = grid_id(boundary)
+    b_indices = boundary_indices(g, boundary)
+
+    # Compute the weights for the logical derivatives
+    g⁻¹ = metric_tensor_inverse(g)
+    α = map(CartesianIndices(g⁻¹)[b_indices...]) do I # TODO: Fix iterator here
+        gᵏⁱ = g⁻¹[I][k,:]
+        gᵏᵏ = g⁻¹[I][k,k]
+
+        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)
+        end
+
+        αᵢ = componentview(α,i)
+        DiagonalTensor(αᵢ) ∘ ∂_ξᵢ
+    end
+end

--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Wed Aug 28 10:50:15 2024 +0200
@@ -50,3 +50,8 @@
 """
 inner_product(g::ZeroDimGrid, stencil_set::StencilSet) = IdentityTensor{component_type(g)}()
 
+
+function inner_product(g::MappedGrid, stencil_set)
+    J = map(sqrt∘det, metric_tensor(g))
+    DiagonalTensor(J)∘inner_product(logicalgrid(g), stencil_set)
+end

--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Wed Aug 28 10:50:15 2024 +0200
@@ -49,3 +49,8 @@
 Implemented to simplify 1D code for SBP operators.
 """
 inverse_inner_product(g::ZeroDimGrid, stencil_set::StencilSet) = IdentityTensor{component_type(g)}()
+
+function inverse_inner_product(g::MappedGrid, stencil_set)
+    J⁻¹ = map(inv∘sqrt∘det, metric_tensor(g))
+    DiagonalTensor(J⁻¹)∘inner_product(logicalgrid(g), stencil_set)
+end

--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Wed Aug 28 10:50:15 2024 +0200
@@ -51,8 +51,31 @@
     end
     return Δ
 end
+
 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
 
+function laplace(grid::MappedGrid, stencil_set)
+    J = jacobian_determinant(grid)
+    J⁻¹ = DiagonalTensor(map(inv, J))
+
+    Jg = map(*, J, metric_tensor_inverse(grid))
+    lg = logicalgrid(grid)
+
+    return mapreduce(+, CartesianIndices(first(Jg))) do I
+        i, j = I[1], I[2]
+        Jgⁱʲ = componentview(Jg, i, j)
+
+        if i == j
+            J⁻¹∘second_derivative_variable(lg, Jgⁱʲ, stencil_set, i)
+        else
+            Dᵢ = first_derivative(lg, stencil_set, i)
+            Dⱼ = first_derivative(lg, stencil_set, j)
+            J⁻¹∘Dᵢ∘DiagonalTensor(Jgⁱʲ)∘Dⱼ
+        end
+    end
+end
+
+
 """
     sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
 

--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Wed Aug 28 10:50:15 2024 +0200
@@ -6,6 +6,8 @@
 using Sbplib.RegionIndices
 using Sbplib.SbpOperators: BoundaryOperator, Stencil
 
+using StaticArrays
+
 @testset "boundary_restriction" begin
 	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4)
 	e_closure = parse_stencil(stencil_set["e"]["closure"])
@@ -33,7 +35,7 @@
     end
 
     @testset "Application" begin
-        @testset "1D" begin
+        @testset "EquidistantGrid" begin
             e_l, e_r = boundary_restriction.(Ref(g_1D), Ref(stencil_set), boundary_identifiers(g_1D))
             v = eval_on(g_1D,x->1+x^2)
             u = fill(3.124)
@@ -43,7 +45,7 @@
             @test (e_r*v)[1] == v[end]
         end
 
-        @testset "2D" begin
+        @testset "TensorGrid" begin
             e_w, e_e, e_s, e_n = boundary_restriction.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D))
             v = rand(11, 15)
             u = fill(3.124)
@@ -53,5 +55,22 @@
             @test e_s*v == v[:,1]
             @test e_n*v == v[:,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)
+
+            e_w, e_e, e_s, e_n = boundary_restriction.(Ref(mg), Ref(stencil_set), boundary_identifiers(mg))
+            v = rand(10, 13)
+
+            @test e_w*v == v[1,:]
+            @test e_e*v == v[end,:]
+            @test e_s*v == v[:,1]
+            @test e_n*v == v[:,end]
+       end
     end
 end

--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Wed Aug 28 10:50:15 2024 +0200
@@ -6,54 +6,130 @@
 using Sbplib.RegionIndices
 import Sbplib.SbpOperators.BoundaryOperator
 
+using StaticArrays
+using LinearAlgebra
+
 @testset "normal_derivative" begin
-    g_1D = equidistant_grid(0.0, 1.0, 11)
-    g_2D = equidistant_grid((0.0, 0.0), (1.0,1.0), 11, 12)
-    @testset "normal_derivative" begin
-    	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-        @testset "1D" begin
-            d_l = normal_derivative(g_1D, stencil_set, Lower())
-            @test d_l == normal_derivative(g_1D, stencil_set, Lower())
-            @test d_l isa BoundaryOperator{T,Lower} where T
-            @test d_l isa LazyTensor{T,0,1} where T
-        end
-        @testset "2D" begin
-            d_w = normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}())
-            d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,Upper}())
-            Ix = IdentityTensor{Float64}((size(g_2D)[1],))
-            Iy = IdentityTensor{Float64}((size(g_2D)[2],))
-            d_l = normal_derivative(g_2D.grids[1], stencil_set, Lower())
-            d_r = normal_derivative(g_2D.grids[2], stencil_set, Upper())
-            @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}())
-            @test d_w ==  d_l⊗Iy
-            @test d_n ==  Ix⊗d_r
-            @test d_w isa LazyTensor{T,1,2} where T
-            @test d_n isa LazyTensor{T,1,2} where T
+	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+
+    @testset "EquidistantGrid" begin
+        g_1D = equidistant_grid(0.0, 1.0, 11)
+
+        d_l = normal_derivative(g_1D, stencil_set, Lower())
+        @test d_l == normal_derivative(g_1D, stencil_set, Lower())
+        @test d_l isa BoundaryOperator{T,Lower} where T
+        @test d_l isa LazyTensor{T,0,1} where T
+    end
+
+    @testset "TensorGrid" begin
+        g_2D = equidistant_grid((0.0, 0.0), (1.0,1.0), 11, 12)
+        d_w = normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}())
+        d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,Upper}())
+        Ix = IdentityTensor{Float64}((size(g_2D)[1],))
+        Iy = IdentityTensor{Float64}((size(g_2D)[2],))
+        d_l = normal_derivative(g_2D.grids[1], stencil_set, Lower())
+        d_r = normal_derivative(g_2D.grids[2], stencil_set, Upper())
+        @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}())
+        @test d_w ==  d_l⊗Iy
+        @test d_n ==  Ix⊗d_r
+        @test d_w isa LazyTensor{T,1,2} where T
+        @test d_n isa LazyTensor{T,1,2} where T
+
+        @testset "Accuracy" begin
+            v = eval_on(g_2D, (x,y)-> x^2 + (y-1)^2 + x*y)
+            v∂x = eval_on(g_2D, (x,y)-> 2*x + y)
+            v∂y = eval_on(g_2D, (x,y)-> 2*(y-1) + x)
+            # TODO: Test for higher order polynomials?
+            @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(g_2D), Ref(stencil_set), boundary_identifiers(g_2D))
+
+                @test d_w*v ≈ -v∂x[1,:] atol = 1e-13
+                @test d_e*v ≈ v∂x[end,:] atol = 1e-13
+                @test d_s*v ≈ -v∂y[:,1] atol = 1e-13
+                @test d_n*v ≈ v∂y[:,end] 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(g_2D), Ref(stencil_set), boundary_identifiers(g_2D))
+
+                @test d_w*v ≈ -v∂x[1,:] atol = 1e-13
+                @test d_e*v ≈ v∂x[end,:] atol = 1e-13
+                @test d_s*v ≈ -v∂y[:,1] atol = 1e-13
+                @test d_n*v ≈ v∂y[:,end] atol = 1e-13
+            end
         end
     end
-    @testset "Accuracy" begin
-        v = eval_on(g_2D, (x,y)-> x^2 + (y-1)^2 + x*y)
-        v∂x = eval_on(g_2D, (x,y)-> 2*x + y)
-        v∂y = eval_on(g_2D, (x,y)-> 2*(y-1) + x)
-        # TODO: Test for higher order polynomials?
-        @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(g_2D), Ref(stencil_set), boundary_identifiers(g_2D))
+
+    @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)
+
 
-            @test d_w*v ≈ -v∂x[1,:] atol = 1e-13
-            @test d_e*v ≈ v∂x[end,:] atol = 1e-13
-            @test d_s*v ≈ -v∂y[:,1] atol = 1e-13
-            @test d_n*v ≈ v∂y[:,end] atol = 1e-13
+        # 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)
+
+        for bid ∈ boundary_identifiers(mg)
+            @testset let bid=bid
+                @test normal_derivative(mg, stencil_set, bid) isa LazyTensor{<:Any, 1, 2}
+            end
         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(g_2D), Ref(stencil_set), boundary_identifiers(g_2D))
-            
-            @test d_w*v ≈ -v∂x[1,:] atol = 1e-13
-            @test d_e*v ≈ v∂x[end,:] atol = 1e-13
-            @test d_s*v ≈ -v∂y[:,1] atol = 1e-13
-            @test d_n*v ≈ v∂y[:,end] atol = 1e-13
+        @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 = function(x̄)
+                sin(norm(x̄+@SVector[1,1]))
+            end
+            ∇v = function(x̄)
+                ȳ = x̄+@SVector[1,1]
+                cos(norm(ȳ))*(ȳ/norm(ȳ))
+            end
+
+            mg = equidistant_grid(c, 80,80)
+            v̄ = map(v, mg)
+
+            @testset for (order, atol) ∈ [(2,4e-2),(4,2e-3)]
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=order)
+
+                @testset for bId ∈ boundary_identifiers(mg)
+                    ∂ₙv = map(boundary_grid(mg,bId),normal(mg,bId)) do x̄,n̂
+                        n̂⋅∇v(x̄)
+                    end
+
+                    dₙ = normal_derivative(mg, stencil_set, bId)
+                    @test dₙ*v̄ ≈ ∂ₙv atol=atol
+                end
+            end
         end
     end
 end

--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Wed Aug 28 10:50:15 2024 +0200
@@ -6,6 +6,9 @@
 
 import Sbplib.SbpOperators.ConstantInteriorScalingOperator
 
+using StaticArrays
+using LinearAlgebra
+
 @testset "Diagonal-stencil inner_product" begin
     Lx = π/2.
     Ly = Float64(π)
@@ -94,4 +97,43 @@
             end
         end
     end
+
+    @testset "MappedGrid" begin
+        stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        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)
+
+        @test inner_product(mg, stencil_set) isa LazyTensor{<:Any, 2,2}
+
+        @testset "Accuracy" begin
+            v = function(x̄)
+                log(norm(x̄-@SVector[.5, .5]))/2π + log(norm(x̄-@SVector[1.5, 3]))/2π
+            end
+            ∇v = function(x̄)
+                ∇log(ȳ) = ȳ/(ȳ⋅ȳ)
+                ∇log(x̄-@SVector[.5, .5])/2π + ∇log(x̄-@SVector[1.5, 3])/2π
+            end
+
+            mg = equidistant_grid(c, 80,80)
+            v̄ = map(v, mg)
+
+            @testset for (order, atol) ∈ [(2,1e-3),(4,1e-7)]
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=order)
+
+                @test sum(boundary_identifiers(mg)) do bId
+                    ∂ₙv = map(boundary_grid(mg,bId),normal(mg,bId)) do x̄,n̂
+                        n̂⋅∇v(x̄)
+                    end
+                    Hᵧ = inner_product(boundary_grid(mg,bId), stencil_set)
+                    sum(Hᵧ*∂ₙv)
+                end ≈ 2 atol=atol
+
+            end
+        end
+        @test_broken false # Test that it calculates the right thing
+    end
 end

--- a/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Wed Aug 28 10:50:15 2024 +0200
@@ -6,6 +6,8 @@
 
 import Sbplib.SbpOperators.ConstantInteriorScalingOperator
 
+using StaticArrays
+
 @testset "Diagonal-stencil inverse_inner_product" begin
     Lx = π/2.
     Ly = Float64(π)
@@ -82,4 +84,17 @@
             end
         end
     end
+
+    @testset "MappedGrid" begin
+        stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        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)
+
+        @test inverse_inner_product(mg, stencil_set) isa LazyTensor{<:Any, 2,2}
+        @test_broken false # Test that it calculates the right thing
+    end
 end

--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Wed Aug 28 10:35:08 2024 +0200
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Wed Aug 28 10:50:15 2024 +0200
@@ -4,6 +4,8 @@
 using Sbplib.Grids
 using Sbplib.LazyTensors
 
+using StaticArrays
+
 @testset "Laplace" begin
     # Default stencils (4th order)
     operator_path = sbp_operators_path()*"standard_diagonal.toml"
@@ -72,12 +74,12 @@
     g_1D = equidistant_grid(0.0, 1., 101)
     g_3D = equidistant_grid((0.0, -1.0, 0.0), (1., 1., 1.), 51, 101, 52)
 
-    @testset "1D" begin
+    @testset "EquidistantGrid" begin
         Δ = laplace(g_1D, stencil_set)
         @test Δ == second_derivative(g_1D, stencil_set)
         @test Δ isa LazyTensor{Float64,1,1}
     end
-    @testset "3D" begin
+    @testset "TensorGrid" begin
         Δ = laplace(g_3D, stencil_set)
         @test Δ isa LazyTensor{Float64,3,3}
         Dxx = second_derivative(g_3D, stencil_set, 1)
@@ -86,6 +88,26 @@
         @test Δ == Dxx + Dyy + Dzz
         @test Δ isa LazyTensor{Float64,3,3}
     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]
+
+        g = equidistant_grid(c, 30,30)
+
+        @test laplace(g, stencil_set) isa LazyTensor{<:Any,2,2}
+
+        f((x,y)) = sin(4(x + y))
+        Δf((x,y)) = -16sin(4(x + y))
+        gf = map(f,g)
+
+        Δ = laplace(g, stencil_set)
+
+        @test collect(Δ*gf) isa Array{<:Any,2}
+        @test Δ*gf ≈ map(Δf, g)
+    end
 end
 
 @testset "sat_tensors" begin