--- a/Notes.md	Mon Feb 21 10:38:19 2022 +0100
+++ b/Notes.md	Mon Feb 21 23:36:41 2022 +0100
@@ -147,6 +147,7 @@
  - [ ] How do we handle mixes of periodic and non-periodic grids? Seems it should be supported on the grid level and on the 1d operator level. Between there it should be transparent.
  - [ ] Can we have a trait to tell if a TensorMapping is transposable?
  - [ ] Is it ok to have "Constructors" for abstract types which create subtypes? For example a Grids() functions that gives different kind of grids based on input?
+ - [ ] Figure out how to treat the borrowing parameters of operators. Include in into the struct? Expose via function dispatched on the operator type and grid?
 
 ## Regions and tensormappings
 - [ ] Use a trait to indicate if a TensorMapping uses indices with regions.

--- a/src/Grids/AbstractGrid.jl	Mon Feb 21 10:38:19 2022 +0100
+++ b/src/Grids/AbstractGrid.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -7,7 +7,7 @@
 
 """
 abstract type AbstractGrid end
-
+export AbstractGrid
 function dimension end
 function points end
 export dimension, points

--- a/src/SbpOperators/SbpOperators.jl	Mon Feb 21 10:38:19 2022 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -21,4 +21,14 @@
 include("boundaryops/boundary_restriction.jl")
 include("boundaryops/normal_derivative.jl")
 
+
+export boundary_quadrature
+export boundary_restriction
+export inner_product
+export inverse_inner_product
+export Laplace
+export laplace
+export normal_derivative
+export second_derivative
+
 end # module

--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Mon Feb 21 10:38:19 2022 +0100
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -1,18 +1,22 @@
 """
-    boundary_restriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
-    boundary_restriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
+    boundary_restriction(grid, closure_stencil::Stencil, boundary)
 
-Creates the boundary restriction operator `e` as a `TensorMapping`
+Creates boundary restriction operators `e` as `TensorMapping`s on `boundary`
 
-`e` is the restriction of a grid function to the boundary specified by `boundary` or `region` using some `closure_stencil`.
-`e'` is the prolongation of a grid function on the boundary to the whole grid using the same `closure_stencil`.
+`e` is the restriction of a grid function to `boundary` using a `Stencil` `closure_stencil`.
+`e'` is the prolongation of a grid function on `boundary`` to the whole grid using the same `closure_stencil`.
 On a one-dimensional `grid`, `e` is a `BoundaryOperator`. On a multi-dimensional `grid`, `e` is the inflation of
-a `BoundaryOperator`. Also see the documentation of `SbpOperators.boundary_operator(...)` for more details.
+a `BoundaryOperator`. See also [`SbpOperators.boundary_operator`](@ref).
 """
-function boundary_restriction(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary)
+function boundary_restriction(grid, closure_stencil::Stencil, boundary)
     converted_stencil = convert(Stencil{eltype(grid)}, closure_stencil)
     return SbpOperators.boundary_operator(grid, converted_stencil, boundary)
 end
-boundary_restriction(grid::EquidistantGrid{1}, closure_stencil, region::Region) = boundary_restriction(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
+
+"""
+    boundary_restriction(grid, stencil_set, boundary)
 
-export boundary_restriction
+Creates a `boundary_restriction` operator on `grid` given a parsed TOML
+`stencil_set`.
+"""
+boundary_restriction(grid, stencil_set, boundary) = boundary_restriction(grid, parse_stencil(stencil_set["e"]["closure"]), boundary)

--- a/src/SbpOperators/boundaryops/normal_derivative.jl	Mon Feb 21 10:38:19 2022 +0100
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -1,18 +1,23 @@
 """
-    normal_derivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
-    normal_derivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
+    normal_derivative(grid, closure_stencil::Stencil, boundary)
 
 Creates the normal derivative boundary operator `d` as a `TensorMapping`
 
-`d` is the normal derivative of a grid function at the boundary specified by `boundary` or `region` using some `closure_stencil`.
+`d` computes the normal derivative of a grid function  on `boundary` a `Stencil` `closure_stencil`.
 `d'` is the prolongation of the normal derivative of a grid function to the whole grid using the same `closure_stencil`.
 On a one-dimensional `grid`, `d` is a `BoundaryOperator`. On a multi-dimensional `grid`, `d` is the inflation of
-a `BoundaryOperator`. Also see the documentation of `SbpOperators.boundary_operator(...)` for more details.
+a `BoundaryOperator`. See also [`SbpOperators.boundary_operator`](@ref).
 """
-function normal_derivative(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary)
+function normal_derivative(grid, closure_stencil::Stencil, boundary)
     direction = dim(boundary)
     h_inv = inverse_spacing(grid)[direction]
     return SbpOperators.boundary_operator(grid, scale(closure_stencil,h_inv), boundary)
 end
-normal_derivative(grid::EquidistantGrid{1}, closure_stencil, region::Region) = normal_derivative(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
-export normal_derivative
+
+"""
+    normal_derivative(grid, stencil_set, boundary)
+
+Creates a `normal_derivative` operator on `grid` given a parsed TOML
+`stencil_set`.
+"""
+normal_derivative(grid, stencil_set, boundary) = normal_derivative(grid, parse_stencil(stencil_set["e"]["closure"]), boundary)

--- a/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Mon Feb 21 10:38:19 2022 +0100
+++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -16,4 +16,3 @@
     return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even, direction)
 end
 second_derivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils) = second_derivative(grid,inner_stencil,closure_stencils,1)
-export second_derivative

--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl	Mon Feb 21 10:38:19 2022 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -10,8 +10,8 @@
 
 On a 1-dimensional grid, `H` is a `ConstantInteriorScalingOperator`. On a
 N-dimensional grid, `H` is the outer product of the 1-dimensional inner
-product operators for each coordinate direction. Also see the documentation of
-On a 0-dimensional grid, `H` is a 0-dimensional `IdentityMapping`.
+product operators for each coordinate direction. On a 0-dimensional grid,
+`H` is a 0-dimensional `IdentityMapping`.
 """
 function inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
     Hs = ()
@@ -22,7 +22,6 @@
 
     return foldl(⊗, Hs)
 end
-export inner_product
 
 function inner_product(grid::EquidistantGrid{1}, interior_weight, closure_weights)
     h = spacing(grid)[1]
@@ -32,3 +31,15 @@
 end
 
 inner_product(grid::EquidistantGrid{0}, interior_weight, closure_weights) = IdentityMapping{eltype(grid)}()
+
+"""
+    inner_product(grid, stencil_set)
+
+Creates a `inner_product` operator on `grid` given a parsed TOML
+`stencil_set`.
+"""
+function inner_product(grid, stencil_set)
+    inner_stencil = parse_scalar(stencil_set["H"]["inner"])
+    closure_stencils = parse_tuple(stencil_set["H"]["closure"])
+    return inner_product(grid, inner_stencil, closure_stencils)
+end
\ No newline at end of file

--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Mon Feb 21 10:38:19 2022 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -25,6 +25,17 @@
     H⁻¹ = SbpOperators.ConstantInteriorScalingOperator(grid, h⁻¹*1/interior_weight, h⁻¹./closure_weights)
     return H⁻¹
 end
-export inverse_inner_product
 
 inverse_inner_product(grid::EquidistantGrid{0}, interior_weight, closure_weights) = IdentityMapping{eltype(grid)}()
+
+"""
+    inverse_inner_product(grid, stencil_set)
+
+Creates a `inverse_inner_product` operator on `grid` given a parsed TOML
+`stencil_set`.
+"""
+function inverse_inner_product(grid, stencil_set)
+    inner_stencil = parse_scalar(stencil_set["H"]["inner"])
+    closure_stencils = parse_tuple(stencil_set["H"]["closure"])
+    return inverse_inner_product(grid, inner_stencil, closure_stencils)
+end 
\ No newline at end of file

--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Mon Feb 21 10:38:19 2022 +0100
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -1,5 +1,37 @@
 """
-    laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils)
+    Laplace{T, Dim, DiffOp} <: TensorMapping{T, Dim, Dim}
+
+Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a
+`TensorMapping`. Additionally `Laplace` stores the stencil set (parsed from TOML) 
+used to construct the `TensorMapping`.
+"""
+struct Laplace{T, Dim, DiffOp<:TensorMapping{T, Dim, Dim}} <: TensorMapping{T, Dim, Dim}
+    D::DiffOp# Differential operator
+    stencil_set # Stencil set of the operator
+end
+
+"""
+    `Laplace(grid::Equidistant, stencil_set)`
+
+Creates the `Laplace`` operator `Δ` on `grid` given a parsed TOML
+`stencil_set`. See also [`laplace`](@ref).
+"""
+function Laplace(grid::EquidistantGrid, stencil_set)
+    inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+    closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+    Δ = laplace(grid, inner_stencil,closure_stencils)
+    return Laplace(Δ,stencil_set)
+end
+
+LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D)
+LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D)
+LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...)
+
+# TODO: Implement pretty printing of Laplace once pretty printing of TensorMappings is implemented. 
+# Base.show(io::IO, L::Laplace) = ...
+
+"""
+    laplace(grid::EquidistantGrid, inner_stencil, closure_stencils)
 
 Creates the Laplace operator operator `Δ` as a `TensorMapping`
 
@@ -9,7 +41,7 @@
 
 On a one-dimensional `grid`, `Δ` is equivalent to `second_derivative`. On a
 multi-dimensional `grid`, `Δ` is the sum of multi-dimensional `second_derivative`s
-where the sum is carried out lazily.
+where the sum is carried out lazily.  See also [`second_derivative`](@ref).
 """
 function laplace(grid::EquidistantGrid, inner_stencil, closure_stencils)
     Δ = second_derivative(grid, inner_stencil, closure_stencils, 1)
@@ -17,5 +49,4 @@
         Δ += second_derivative(grid, inner_stencil, closure_stencils, d)
     end
     return Δ
-end
-export laplace
+end
\ No newline at end of file

--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl	Mon Feb 21 10:38:19 2022 +0100
+++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -32,14 +32,8 @@
             @test e_w isa TensorMapping{T,1,2} where T
         end
     end
-
-    op_l = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Lower}())
-    op_r = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Upper}())
-
-    op_w = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Lower}())
-    op_e = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Upper}())
-    op_s = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Lower}())
-    op_n = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Upper}())
+    (op_l, op_r) = map(id -> boundary_operator(g_1D, closure_stencil, id), boundary_identifiers(g_1D))
+    (op_w, op_e, op_s, op_n) = map(id -> boundary_operator(g_2D, closure_stencil, id), boundary_identifiers(g_2D))
 
     @testset "Sizes" begin
         @testset "1D" begin

--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Mon Feb 21 10:38:19 2022 +0100
+++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -2,7 +2,6 @@
 
 using Sbplib.SbpOperators
 using Sbplib.Grids
-using Sbplib.RegionIndices
 using Sbplib.LazyTensors
 
 import Sbplib.SbpOperators.BoundaryOperator
@@ -15,14 +14,12 @@
 
     @testset "boundary_restriction" begin
         @testset "1D" begin
-            e_l = boundary_restriction(g_1D,e_closure,Lower())
-            @test e_l == boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Lower}())
+            e_l = boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Lower}())
             @test e_l == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),Lower())
             @test e_l isa BoundaryOperator{T,Lower} where T
             @test e_l isa TensorMapping{T,0,1} where T
 
-            e_r = boundary_restriction(g_1D,e_closure,Upper())
-            @test e_r == boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Upper}())
+            e_r = boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Upper}())
             @test e_r == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),Upper())
             @test e_r isa BoundaryOperator{T,Upper} where T
             @test e_r isa TensorMapping{T,0,1} where T
@@ -37,8 +34,8 @@
 
     @testset "Application" begin
         @testset "1D" begin
-            e_l = boundary_restriction(g_1D, e_closure, CartesianBoundary{1,Lower}())
-            e_r = boundary_restriction(g_1D, e_closure, CartesianBoundary{1,Upper}())
+            (e_l, e_r) = 
+                map(id -> boundary_restriction(g_1D, e_closure, id), boundary_identifiers(g_1D))
 
             v = evalOn(g_1D,x->1+x^2)
             u = fill(3.124)
@@ -49,10 +46,8 @@
         end
 
         @testset "2D" begin
-            e_w = boundary_restriction(g_2D, e_closure, CartesianBoundary{1,Lower}())
-            e_e = boundary_restriction(g_2D, e_closure, CartesianBoundary{1,Upper}())
-            e_s = boundary_restriction(g_2D, e_closure, CartesianBoundary{2,Lower}())
-            e_n = boundary_restriction(g_2D, e_closure, CartesianBoundary{2,Upper}())
+            (e_w, e_e, e_s, e_n) = 
+                map(id -> boundary_restriction(g_2D, e_closure, id), boundary_identifiers(g_2D))
 
             v = rand(11, 15)
             u = fill(3.124)

--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl	Mon Feb 21 10:38:19 2022 +0100
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -2,7 +2,6 @@
 
 using Sbplib.SbpOperators
 using Sbplib.Grids
-using Sbplib.RegionIndices
 using Sbplib.LazyTensors
 
 import Sbplib.SbpOperators.BoundaryOperator
@@ -14,8 +13,7 @@
     	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
     	d_closure = parse_stencil(stencil_set["d1"]["closure"])
         @testset "1D" begin
-            d_l = normal_derivative(g_1D, d_closure, Lower())
-            @test d_l == normal_derivative(g_1D, d_closure, CartesianBoundary{1,Lower}())
+            d_l = normal_derivative(g_1D, d_closure, CartesianBoundary{1,Lower}())
             @test d_l isa BoundaryOperator{T,Lower} where T
             @test d_l isa TensorMapping{T,0,1} where T
         end
@@ -24,8 +22,8 @@
             d_n = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Upper}())
             Ix = IdentityMapping{Float64}((size(g_2D)[1],))
             Iy = IdentityMapping{Float64}((size(g_2D)[2],))
-            d_l = normal_derivative(restrict(g_2D,1),d_closure,Lower())
-            d_r = normal_derivative(restrict(g_2D,2),d_closure,Upper())
+            d_l = normal_derivative(restrict(g_2D,1),d_closure,CartesianBoundary{1,Lower}())
+            d_r = normal_derivative(restrict(g_2D,2),d_closure,CartesianBoundary{1,Upper}())
             @test d_w ==  d_l⊗Iy
             @test d_n ==  Ix⊗d_r
             @test d_w isa TensorMapping{T,1,2} where T
@@ -40,10 +38,8 @@
         @testset "2nd order" begin
         	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
         	d_closure = parse_stencil(stencil_set["d1"]["closure"])
-            d_w = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Lower}())
-            d_e = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Upper}())
-            d_s = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Lower}())
-            d_n = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Upper}())
+            (d_w, d_e, d_s, d_n) = 
+                map(id -> normal_derivative(g_2D, d_closure, id), boundary_identifiers(g_2D))
 
             @test d_w*v ≈ -v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ v∂x[end,:] atol = 1e-13
@@ -54,10 +50,8 @@
         @testset "4th order" begin
             stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         	d_closure = parse_stencil(stencil_set["d1"]["closure"])
-            d_w = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Lower}())
-            d_e = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Upper}())
-            d_s = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Lower}())
-            d_n = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Upper}())
+            (d_w, d_e, d_s, d_n) = 
+                map(id -> normal_derivative(g_2D, d_closure, id), boundary_identifiers(g_2D))
 
             @test d_w*v ≈ -v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ v∂x[end,:] atol = 1e-13

--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Mon Feb 21 10:38:19 2022 +0100
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -4,25 +4,25 @@
 using Sbplib.Grids
 using Sbplib.LazyTensors
 
+# Default stencils (4th order)
+operator_path = sbp_operators_path()*"standard_diagonal.toml"
+stencil_set = read_stencil_set(operator_path; order=4)
+inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
+closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
+g_1D = EquidistantGrid(101, 0.0, 1.)
+g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
+
 @testset "Laplace" begin
-    g_1D = EquidistantGrid(101, 0.0, 1.)
-    g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
     @testset "Constructors" begin
-        stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-        inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
-        closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
         @testset "1D" begin
-            L = laplace(g_1D, inner_stencil, closure_stencils)
-            @test L == second_derivative(g_1D, inner_stencil, closure_stencils)
-            @test L isa TensorMapping{T,1,1}  where T
+            Δ = laplace(g_1D, inner_stencil, closure_stencils)
+            @test Laplace(g_1D, stencil_set) == Laplace(Δ, stencil_set)
+            @test Laplace(g_1D, stencil_set) isa TensorMapping{T,1,1}  where T
         end
         @testset "3D" begin
-            L = laplace(g_3D, inner_stencil, closure_stencils)
-            @test L isa TensorMapping{T,3,3} where T
-            Dxx = second_derivative(g_3D, inner_stencil, closure_stencils, 1)
-            Dyy = second_derivative(g_3D, inner_stencil, closure_stencils, 2)
-            Dzz = second_derivative(g_3D, inner_stencil, closure_stencils, 3)
-            @test L == Dxx + Dyy + Dzz
+            Δ = laplace(g_3D, inner_stencil, closure_stencils)
+            @test Laplace(g_3D, stencil_set) == Laplace(Δ,stencil_set)
+            @test Laplace(g_3D, stencil_set) isa TensorMapping{T,3,3} where T
         end
     end
 
@@ -42,30 +42,44 @@
         # 2nd order interior stencil, 1st order boundary stencil,
         # implies that L*v should be exact for binomials up to order 2.
         @testset "2nd order" begin
-            stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-            inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
-            closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
-            L = laplace(g_3D, inner_stencil, closure_stencils)
-            @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
-            @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
-            @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
-            @test L*v ≈ Δv rtol = 5e-2 norm = l2
+            stencil_set = read_stencil_set(operator_path; order=2)
+            Δ = Laplace(g_3D, stencil_set)
+            @test Δ*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
+            @test Δ*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
+            @test Δ*polynomials[3] ≈ polynomials[1] atol = 5e-9
+            @test Δ*v ≈ Δv rtol = 5e-2 norm = l2
         end
 
         # 4th order interior stencil, 2nd order boundary stencil,
         # implies that L*v should be exact for binomials up to order 3.
         @testset "4th order" begin
-            stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-            inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
-            closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
-            L = laplace(g_3D, inner_stencil, closure_stencils)
+            stencil_set = read_stencil_set(operator_path; order=4)
+            Δ = Laplace(g_3D, stencil_set)
             # NOTE: high tolerances for checking the "exact" differentiation
             # due to accumulation of round-off errors/cancellation errors?
-            @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
-            @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
-            @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
-            @test L*polynomials[4] ≈ polynomials[2] atol = 5e-9
-            @test L*v ≈ Δv rtol = 5e-4 norm = l2
+            @test Δ*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
+            @test Δ*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
+            @test Δ*polynomials[3] ≈ polynomials[1] atol = 5e-9
+            @test Δ*polynomials[4] ≈ polynomials[2] atol = 5e-9
+            @test Δ*v ≈ Δv rtol = 5e-4 norm = l2
         end
     end
 end
+
+@testset "laplace" begin
+    @testset "1D" begin
+        Δ = laplace(g_1D, inner_stencil, closure_stencils)
+        @test Δ == second_derivative(g_1D, inner_stencil, closure_stencils)
+        @test Δ isa TensorMapping{T,1,1}  where T
+    end
+    @testset "3D" begin
+        Δ = laplace(g_3D, inner_stencil, closure_stencils)
+        @test Δ isa TensorMapping{T,3,3} where T
+        Dxx = second_derivative(g_3D, inner_stencil, closure_stencils, 1)
+        Dyy = second_derivative(g_3D, inner_stencil, closure_stencils, 2)
+        Dzz = second_derivative(g_3D, inner_stencil, closure_stencils, 3)
+        @test Δ == Dxx + Dyy + Dzz
+        @test Δ isa TensorMapping{T,3,3} where T
+    end
+end
+