--- a/Manifest.toml	Wed Jun 26 12:42:28 2024 +0200
+++ b/Manifest.toml	Wed Jun 26 12:47:26 2024 +0200
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.4"
 manifest_format = "2.0"
-project_hash = "a36735c53cfa4453f39635046eeaa47a4ea1231b"
+project_hash = "32fac879810099260f177c27318d3f26de4a00cc"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]

--- a/Project.toml	Wed Jun 26 12:42:28 2024 +0200
+++ b/Project.toml	Wed Jun 26 12:47:26 2024 +0200
@@ -4,17 +4,18 @@
 version = "0.1.1"
 
 [deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
 
 [weakdeps]
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 
 [extensions]
+SbplibMakieExt = "Makie"
 SbplibPlotsExt = "Plots"
-SbplibMakieExt = "Makie"
 
 [compat]
 julia = "1.5"

--- a/src/Grids/Grids.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/src/Grids/Grids.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -4,6 +4,7 @@
 using Sbplib.RegionIndices
 using Sbplib.LazyTensors
 using StaticArrays
+using LinearAlgebra
 
 export ParameterSpace
 export HyperBox

--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Wed Jun 26 12:47:26 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, geometric_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/volumeops/laplace/laplace_test.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Wed Jun 26 12:47:26 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