--- a/test/testSbpOperators.jl	Sat Feb 06 15:26:14 2021 +0100
+++ b/test/testSbpOperators.jl	Mon Feb 15 11:30:34 2021 +0100
@@ -24,9 +24,12 @@
     @test eltype(s) == Float64
     @test SbpOperators.scale(s, 2) == Stencil((-2,2), (2.,4.,4.,6.,8.))
 
-    @test Stencil((1,2,3,4), center=1) == Stencil((0, 3),(1,2,3,4))
-    @test Stencil((1,2,3,4), center=2) == Stencil((-1, 2),(1,2,3,4))
-    @test Stencil((1,2,3,4), center=4) == Stencil((-3, 0),(1,2,3,4))
+    @test Stencil(1,2,3,4; center=1) == Stencil((0, 3),(1,2,3,4))
+    @test Stencil(1,2,3,4; center=2) == Stencil((-1, 2),(1,2,3,4))
+    @test Stencil(1,2,3,4; center=4) == Stencil((-3, 0),(1,2,3,4))
+
+    @test CenteredStencil(1,2,3,4,5) == Stencil((-2, 2), (1,2,3,4,5))
+    @test_throws ArgumentError CenteredStencil(1,2,3,4)
 end
 
 @testset "parse_rational" begin
@@ -67,40 +70,40 @@
 
     parsed_toml = TOML.parse(toml_str)
     @testset "get_stencil" begin
-        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil") == Stencil((-1/2, 0., 1/2), center=2)
-        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=1) == Stencil((-1/2, 0., 1/2); center=1)
-        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=3) == Stencil((-1/2, 0., 1/2); center=3)
+        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil") == Stencil(-1/2, 0., 1/2, center=2)
+        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=1) == Stencil(-1/2, 0., 1/2; center=1)
+        @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=3) == Stencil(-1/2, 0., 1/2; center=3)
 
-        @test get_stencil(parsed_toml, "order2", "H", "inner") == Stencil((1.,), center=1)
+        @test get_stencil(parsed_toml, "order2", "H", "inner") == Stencil(1.; center=1)
 
         @test_throws AssertionError get_stencil(parsed_toml, "meta", "type")
         @test_throws AssertionError get_stencil(parsed_toml, "order2", "D1", "closure_stencils")
     end
 
     @testset "get_stencils" begin
-        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(1,)) == (Stencil((-1., 1.), center=1),)
-        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(2,)) == (Stencil((-1., 1.), center=2),)
-        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=[2]) == (Stencil((-1., 1.), center=2),)
+        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(1,)) == (Stencil(-1., 1., center=1),)
+        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(2,)) == (Stencil(-1., 1., center=2),)
+        @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=[2]) == (Stencil(-1., 1., center=2),)
 
         @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=[1,1,1,1]) == (
-            Stencil((    2.,    -5.,      4.,     -1.,    0.,    0.), center=1),
-            Stencil((    1.,    -2.,      1.,      0.,    0.,    0.), center=1),
-            Stencil(( -4/43,  59/43, -110/43,   59/43, -4/43,    0.), center=1),
-            Stencil(( -1/49,     0.,   59/49, -118/49, 64/49, -4/49), center=1),
+            Stencil(    2.,    -5.,      4.,     -1.,    0.,    0., center=1),
+            Stencil(    1.,    -2.,      1.,      0.,    0.,    0., center=1),
+            Stencil( -4/43,  59/43, -110/43,   59/43, -4/43,    0., center=1),
+            Stencil( -1/49,     0.,   59/49, -118/49, 64/49, -4/49, center=1),
         )
 
         @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(4,2,3,1)) == (
-            Stencil((    2.,    -5.,      4.,     -1.,    0.,    0.), center=4),
-            Stencil((    1.,    -2.,      1.,      0.,    0.,    0.), center=2),
-            Stencil(( -4/43,  59/43, -110/43,   59/43, -4/43,    0.), center=3),
-            Stencil(( -1/49,     0.,   59/49, -118/49, 64/49, -4/49), center=1),
+            Stencil(    2.,    -5.,      4.,     -1.,    0.,    0., center=4),
+            Stencil(    1.,    -2.,      1.,      0.,    0.,    0., center=2),
+            Stencil( -4/43,  59/43, -110/43,   59/43, -4/43,    0., center=3),
+            Stencil( -1/49,     0.,   59/49, -118/49, 64/49, -4/49, center=1),
         )
 
         @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=1:4) == (
-            Stencil((    2.,    -5.,      4.,     -1.,    0.,    0.), center=1),
-            Stencil((    1.,    -2.,      1.,      0.,    0.,    0.), center=2),
-            Stencil(( -4/43,  59/43, -110/43,   59/43, -4/43,    0.), center=3),
-            Stencil(( -1/49,     0.,   59/49, -118/49, 64/49, -4/49), center=4),
+            Stencil(    2.,    -5.,      4.,     -1.,    0.,    0., center=1),
+            Stencil(    1.,    -2.,      1.,      0.,    0.,    0., center=2),
+            Stencil( -4/43,  59/43, -110/43,   59/43, -4/43,    0., center=3),
+            Stencil( -1/49,     0.,   59/49, -118/49, 64/49, -4/49, center=4),
         )
 
         @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3))
@@ -116,8 +119,8 @@
 end
 
 @testset "VolumeOperator" begin
-    inner_stencil = Stencil(1/4 .* (1.,2.,1.),center=2)
-    closure_stencils = (Stencil(1/2 .* (1.,1.),center=1),Stencil((0.,1.),center=2))
+    inner_stencil = CenteredStencil(1/4, 2/4, 1/4)
+    closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(0.,1.; center=2))
     g_1D = EquidistantGrid(11,0.,1.)
     g_2D = EquidistantGrid((11,12),(0.,0.),(1.,1.))
     g_3D = EquidistantGrid((11,12,10),(0.,0.,0.),(1.,1.,1.))
@@ -238,13 +241,13 @@
 
     @testset "Constructors" begin
         @testset "1D" begin
-            Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
-            @test Dₓₓ == SecondDerivative(g_1D,op.innerStencil,op.closureStencils,1)
+            Dₓₓ = second_derivative(g_1D,op.innerStencil,op.closureStencils)
+            @test Dₓₓ == second_derivative(g_1D,op.innerStencil,op.closureStencils,1)
             @test Dₓₓ isa VolumeOperator
         end
         @testset "2D" begin
-            Dₓₓ = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,1)
-            D2 = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
+            Dₓₓ = second_derivative(g_2D,op.innerStencil,op.closureStencils,1)
+            D2 = second_derivative(g_1D,op.innerStencil,op.closureStencils)
             I = IdentityMapping{Float64}(size(g_2D)[2])
             @test Dₓₓ == D2⊗I
             @test Dₓₓ isa TensorMapping{T,2,2} where T
@@ -269,7 +272,7 @@
             # implies that L*v should be exact for monomials up to order 2.
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
+                Dₓₓ = second_derivative(g_1D,op.innerStencil,op.closureStencils)
                 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
                 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
                 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
@@ -280,7 +283,7 @@
             # implies that L*v should be exact for monomials up to order 3.
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
+                Dₓₓ = second_derivative(g_1D,op.innerStencil,op.closureStencils)
                 # NOTE: high tolerances for checking the "exact" differentiation
                 # due to accumulation of round-off errors/cancellation errors?
                 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
@@ -306,7 +309,7 @@
             # implies that L*v should be exact for binomials up to order 2.
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2)
+                Dyy = second_derivative(g_2D,op.innerStencil,op.closureStencils,2)
                 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
                 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
                 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
@@ -317,7 +320,7 @@
             # implies that L*v should be exact for binomials up to order 3.
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2)
+                Dyy = second_derivative(g_2D,op.innerStencil,op.closureStencils,2)
                 # NOTE: high tolerances for checking the "exact" differentiation
                 # due to accumulation of round-off errors/cancellation errors?
                 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
@@ -338,21 +341,21 @@
         @testset "1D" begin
             # Create all tensor mappings included in Laplace
             Δ = laplace(g_1D, op.innerStencil, op.closureStencils)
-            H = quadrature(g_1D, op.quadratureClosure)
-            Hi = InverseDiagonalQuadrature(g_1D, op.quadratureClosure)
+            H = inner_product(g_1D, op.quadratureClosure)
+            Hi = inverse_inner_product(g_1D, op.quadratureClosure)
 
             (id_l, id_r) = boundary_identifiers(g_1D)
 
-            e_l = BoundaryRestriction(g_1D,op.eClosure,id_l)
-            e_r = BoundaryRestriction(g_1D,op.eClosure,id_r)
+            e_l = boundary_restriction(g_1D,op.eClosure,id_l)
+            e_r = boundary_restriction(g_1D,op.eClosure,id_r)
             e_dict = Dict(Pair(id_l,e_l),Pair(id_r,e_r))
 
-            d_l = NormalDerivative(g_1D,op.dClosure,id_l)
-            d_r = NormalDerivative(g_1D,op.dClosure,id_r)
+            d_l = normal_derivative(g_1D,op.dClosure,id_l)
+            d_r = normal_derivative(g_1D,op.dClosure,id_r)
             d_dict = Dict(Pair(id_l,d_l),Pair(id_r,d_r))
 
-            H_l = boundary_quadrature(g_1D,op.quadratureClosure,id_r)
-            H_r = boundary_quadrature(g_1D,op.quadratureClosure,id_r)
+            H_l = inner_product(boundary_grid(g_1D,id_l),op.quadratureClosure)
+            H_r = inner_product(boundary_grid(g_1D,id_r),op.quadratureClosure)
             Hb_dict = Dict(Pair(id_l,H_l),Pair(id_r,H_r))
 
             # TODO: Not sure why this doesnt work? Comparing the fields of
@@ -372,37 +375,37 @@
         @testset "3D" begin
             # Create all tensor mappings included in Laplace
             Δ = laplace(g_3D, op.innerStencil, op.closureStencils)
-            H = quadrature(g_3D, op.quadratureClosure)
-            Hi = InverseDiagonalQuadrature(g_3D, op.quadratureClosure)
+            H = inner_product(g_3D, op.quadratureClosure)
+            Hi = inverse_inner_product(g_3D, op.quadratureClosure)
 
             (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D)
 
-            e_l = BoundaryRestriction(g_3D,op.eClosure,id_l)
-            e_r = BoundaryRestriction(g_3D,op.eClosure,id_r)
-            e_s = BoundaryRestriction(g_3D,op.eClosure,id_s)
-            e_n = BoundaryRestriction(g_3D,op.eClosure,id_n)
-            e_b = BoundaryRestriction(g_3D,op.eClosure,id_b)
-            e_t = BoundaryRestriction(g_3D,op.eClosure,id_t)
+            e_l = boundary_restriction(g_3D,op.eClosure,id_l)
+            e_r = boundary_restriction(g_3D,op.eClosure,id_r)
+            e_s = boundary_restriction(g_3D,op.eClosure,id_s)
+            e_n = boundary_restriction(g_3D,op.eClosure,id_n)
+            e_b = boundary_restriction(g_3D,op.eClosure,id_b)
+            e_t = boundary_restriction(g_3D,op.eClosure,id_t)
             e_dict = Dict(Pair(id_l,e_l),Pair(id_r,e_r),
                           Pair(id_s,e_s),Pair(id_n,e_n),
                           Pair(id_b,e_b),Pair(id_t,e_t))
 
-            d_l = NormalDerivative(g_3D,op.dClosure,id_l)
-            d_r = NormalDerivative(g_3D,op.dClosure,id_r)
-            d_s = NormalDerivative(g_3D,op.dClosure,id_s)
-            d_n = NormalDerivative(g_3D,op.dClosure,id_n)
-            d_b = NormalDerivative(g_3D,op.dClosure,id_b)
-            d_t = NormalDerivative(g_3D,op.dClosure,id_t)
+            d_l = normal_derivative(g_3D,op.dClosure,id_l)
+            d_r = normal_derivative(g_3D,op.dClosure,id_r)
+            d_s = normal_derivative(g_3D,op.dClosure,id_s)
+            d_n = normal_derivative(g_3D,op.dClosure,id_n)
+            d_b = normal_derivative(g_3D,op.dClosure,id_b)
+            d_t = normal_derivative(g_3D,op.dClosure,id_t)
             d_dict = Dict(Pair(id_l,d_l),Pair(id_r,d_r),
                           Pair(id_s,d_s),Pair(id_n,d_n),
                           Pair(id_b,d_b),Pair(id_t,d_t))
 
-            H_l = boundary_quadrature(g_3D,op.quadratureClosure,id_r)
-            H_r = boundary_quadrature(g_3D,op.quadratureClosure,id_r)
-            H_s = boundary_quadrature(g_3D,op.quadratureClosure,id_s)
-            H_n = boundary_quadrature(g_3D,op.quadratureClosure,id_n)
-            H_b = boundary_quadrature(g_3D,op.quadratureClosure,id_b)
-            H_t = boundary_quadrature(g_3D,op.quadratureClosure,id_t)
+            H_l = inner_product(boundary_grid(g_3D,id_l),op.quadratureClosure)
+            H_r = inner_product(boundary_grid(g_3D,id_r),op.quadratureClosure)
+            H_s = inner_product(boundary_grid(g_3D,id_s),op.quadratureClosure)
+            H_n = inner_product(boundary_grid(g_3D,id_n),op.quadratureClosure)
+            H_b = inner_product(boundary_grid(g_3D,id_b),op.quadratureClosure)
+            H_t = inner_product(boundary_grid(g_3D,id_t),op.quadratureClosure)
             Hb_dict = Dict(Pair(id_l,H_l),Pair(id_r,H_r),
                           Pair(id_s,H_s),Pair(id_n,H_n),
                           Pair(id_b,H_b),Pair(id_t,H_t))
@@ -426,14 +429,15 @@
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
             L = laplace(g_1D, op.innerStencil, op.closureStencils)
-            @test L == SecondDerivative(g_1D, op.innerStencil, op.closureStencils)
+            @test L == second_derivative(g_1D, op.innerStencil, op.closureStencils)
             @test L isa TensorMapping{T,1,1}  where T
         end
         @testset "3D" begin
             L = laplace(g_3D, op.innerStencil, op.closureStencils)
-            Dxx = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,1)
-            Dyy = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,2)
-            Dzz = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,3)
+            @test L isa TensorMapping{T,3,3} where T
+            Dxx = second_derivative(g_3D, op.innerStencil, op.closureStencils,1)
+            Dyy = second_derivative(g_3D, op.innerStencil, op.closureStencils,2)
+            Dzz = second_derivative(g_3D, op.innerStencil, op.closureStencils,3)
             @test L == Dxx + Dyy + Dzz
             @test L isa TensorMapping{T,3,3} where T
         end
@@ -470,7 +474,8 @@
         # 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
-            L = Laplace(g_3D,sbp_operators_path()*"standard_diagonal.toml"; order=2)
+            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+            L = laplace(g_3D,op.innerStencil,op.closureStencils)
             @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
@@ -480,7 +485,8 @@
         # 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
-            L = Laplace(g_3D,sbp_operators_path()*"standard_diagonal.toml"; order=4)
+            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+            L = laplace(g_3D,op.innerStencil,op.closureStencils)
             # 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
@@ -492,7 +498,7 @@
     end
 end
 
-@testset "Quadrature diagonal" begin
+@testset "Diagonal-stencil inner_product" begin
     Lx = π/2.
     Ly = Float64(π)
     Lz = 1.
@@ -500,65 +506,37 @@
     g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
     g_3D = EquidistantGrid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz))
     integral(H,v) = sum(H*v)
-    @testset "quadrature" begin
+    @testset "inner_product" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        @testset "0D" begin
+            H = inner_product(EquidistantGrid{Float64}(),op.quadratureClosure)
+            @test H == IdentityMapping{Float64}()
+            @test H isa TensorMapping{T,0,0} where T
+        end
         @testset "1D" begin
-            H = quadrature(g_1D,op.quadratureClosure)
-            inner_stencil = Stencil((1.,),center=1)
-            @test H == quadrature(g_1D,inner_stencil,op.quadratureClosure)
+            H = inner_product(g_1D,op.quadratureClosure)
+            inner_stencil = CenteredStencil(1.)
+            @test H == inner_product(g_1D,op.quadratureClosure,inner_stencil)
             @test H isa TensorMapping{T,1,1} where T
         end
         @testset "2D" begin
-            H = quadrature(g_2D,op.quadratureClosure)
-            H_x = quadrature(restrict(g_2D,1),op.quadratureClosure)
-            H_y = quadrature(restrict(g_2D,2),op.quadratureClosure)
+            H = inner_product(g_2D,op.quadratureClosure)
+            H_x = inner_product(restrict(g_2D,1),op.quadratureClosure)
+            H_y = inner_product(restrict(g_2D,2),op.quadratureClosure)
             @test H == H_x⊗H_y
             @test H isa TensorMapping{T,2,2} where T
         end
     end
 
-    @testset "boundary_quadrature" begin
-        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-        @testset "1D" begin
-            (id_l, id_r) = boundary_identifiers(g_1D)
-            @test boundary_quadrature(g_1D,op.quadratureClosure,id_l) == IdentityMapping{Float64}()
-            @test boundary_quadrature(g_1D,op.quadratureClosure,id_r) == IdentityMapping{Float64}()
-
-        end
-        @testset "2D" begin
-            (id_w, id_e, id_s, id_n) = boundary_identifiers(g_2D)
-            H_x = quadrature(restrict(g_2D,1),op.quadratureClosure)
-            H_y = quadrature(restrict(g_2D,2),op.quadratureClosure)
-            @test boundary_quadrature(g_2D,op.quadratureClosure,id_w) == H_y
-            @test boundary_quadrature(g_2D,op.quadratureClosure,id_e) == H_y
-            @test boundary_quadrature(g_2D,op.quadratureClosure,id_s) == H_x
-            @test boundary_quadrature(g_2D,op.quadratureClosure,id_n) == H_x
-        end
-        @testset "3D" begin
-            (id_w, id_e,
-             id_s, id_n,
-             id_t, id_b) = boundary_identifiers(g_3D)
-            H_xy = quadrature(restrict(g_3D,[1,2]),op.quadratureClosure)
-            H_xz = quadrature(restrict(g_3D,[1,3]),op.quadratureClosure)
-            H_yz = quadrature(restrict(g_3D,[2,3]),op.quadratureClosure)
-            @test boundary_quadrature(g_3D,op.quadratureClosure,id_w) == H_yz
-            @test boundary_quadrature(g_3D,op.quadratureClosure,id_e) == H_yz
-            @test boundary_quadrature(g_3D,op.quadratureClosure,id_s) == H_xz
-            @test boundary_quadrature(g_3D,op.quadratureClosure,id_n) == H_xz
-            @test boundary_quadrature(g_3D,op.quadratureClosure,id_t) == H_xy
-            @test boundary_quadrature(g_3D,op.quadratureClosure,id_b) == H_xy
-        end
-    end
-
     @testset "Sizes" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
-            H = quadrature(g_1D,op.quadratureClosure)
+            H = inner_product(g_1D,op.quadratureClosure)
             @test domain_size(H) == size(g_1D)
             @test range_size(H) == size(g_1D)
         end
         @testset "2D" begin
-            H = quadrature(g_2D,op.quadratureClosure)
+            H = inner_product(g_2D,op.quadratureClosure)
             @test domain_size(H) == size(g_2D)
             @test range_size(H) == size(g_2D)
         end
@@ -575,7 +553,7 @@
 
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = quadrature(g_1D,op.quadratureClosure)
+                H = inner_product(g_1D,op.quadratureClosure)
                 for i = 1:2
                     @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
                 end
@@ -584,7 +562,7 @@
 
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = quadrature(g_1D,op.quadratureClosure)
+                H = inner_product(g_1D,op.quadratureClosure)
                 for i = 1:4
                     @test integral(H,v[i]) ≈ v[i+1][end] -  v[i+1][1] rtol = 1e-14
                 end
@@ -598,13 +576,13 @@
             u = evalOn(g_2D,(x,y)->sin(x)+cos(y))
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = quadrature(g_2D,op.quadratureClosure)
+                H = inner_product(g_2D,op.quadratureClosure)
                 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
                 @test integral(H,u) ≈ π rtol = 1e-4
             end
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = quadrature(g_2D,op.quadratureClosure)
+                H = inner_product(g_2D,op.quadratureClosure)
                 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
                 @test integral(H,u) ≈ π rtol = 1e-8
             end
@@ -612,27 +590,32 @@
     end
 end
 
-@testset "InverseDiagonalQuadrature" begin
+@testset "Diagonal-stencil inverse_inner_product" begin
     Lx = π/2.
     Ly = Float64(π)
     g_1D = EquidistantGrid(77, 0.0, Lx)
     g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
-    @testset "Constructors" begin
+    @testset "inverse_inner_product" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        @testset "0D" begin
+            Hi = inverse_inner_product(EquidistantGrid{Float64}(),op.quadratureClosure)
+            @test Hi == IdentityMapping{Float64}()
+            @test Hi isa TensorMapping{T,0,0} where T
+        end
         @testset "1D" begin
-            Hi = InverseDiagonalQuadrature(g_1D, op.quadratureClosure);
-            inner_stencil = Stencil((1.,),center=1)
+            Hi = inverse_inner_product(g_1D, op.quadratureClosure);
+            inner_stencil = CenteredStencil(1.)
             closures = ()
             for i = 1:length(op.quadratureClosure)
                 closures = (closures...,Stencil(op.quadratureClosure[i].range,1.0./op.quadratureClosure[i].weights))
             end
-            @test Hi == InverseQuadrature(g_1D,inner_stencil,closures)
+            @test Hi == inverse_inner_product(g_1D,closures,inner_stencil)
             @test Hi isa TensorMapping{T,1,1} where T
         end
         @testset "2D" begin
-            Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
-            Hi_x = InverseDiagonalQuadrature(restrict(g_2D,1),op.quadratureClosure)
-            Hi_y = InverseDiagonalQuadrature(restrict(g_2D,2),op.quadratureClosure)
+            Hi = inverse_inner_product(g_2D,op.quadratureClosure)
+            Hi_x = inverse_inner_product(restrict(g_2D,1),op.quadratureClosure)
+            Hi_y = inverse_inner_product(restrict(g_2D,2),op.quadratureClosure)
             @test Hi == Hi_x⊗Hi_y
             @test Hi isa TensorMapping{T,2,2} where T
         end
@@ -641,12 +624,12 @@
     @testset "Sizes" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
-            Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
+            Hi = inverse_inner_product(g_1D,op.quadratureClosure)
             @test domain_size(Hi) == size(g_1D)
             @test range_size(Hi) == size(g_1D)
         end
         @testset "2D" begin
-            Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
+            Hi = inverse_inner_product(g_2D,op.quadratureClosure)
             @test domain_size(Hi) == size(g_2D)
             @test range_size(Hi) == size(g_2D)
         end
@@ -658,15 +641,15 @@
             u = evalOn(g_1D,x->x^3-x^2+1)
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = quadrature(g_1D,op.quadratureClosure)
-                Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
+                H = inner_product(g_1D,op.quadratureClosure)
+                Hi = inverse_inner_product(g_1D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = quadrature(g_1D,op.quadratureClosure)
-                Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
+                H = inner_product(g_1D,op.quadratureClosure)
+                Hi = inverse_inner_product(g_1D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
@@ -676,15 +659,15 @@
             u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y))
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = quadrature(g_2D,op.quadratureClosure)
-                Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
+                H = inner_product(g_2D,op.quadratureClosure)
+                Hi = inverse_inner_product(g_2D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = quadrature(g_2D,op.quadratureClosure)
-                Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
+                H = inner_product(g_2D,op.quadratureClosure)
+                Hi = inverse_inner_product(g_2D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
@@ -705,7 +688,7 @@
             @test op_l isa TensorMapping{T,0,1} where T
 
             op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1])
-            @test op_r == BoundaryRestriction(g_1D,closure_stencil,Upper())
+            @test op_r == BoundaryOperator(g_1D,closure_stencil,Upper())
             @test op_r == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Upper}())
             @test op_r isa TensorMapping{T,0,1} where T
         end
@@ -836,28 +819,28 @@
 
 end
 
-@testset "BoundaryRestriction" begin
+@testset "boundary_restriction" begin
     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
     g_1D = EquidistantGrid(11, 0.0, 1.0)
     g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
 
-    @testset "Constructors" begin
+    @testset "boundary_restriction" begin
         @testset "1D" begin
-            e_l = BoundaryRestriction(g_1D,op.eClosure,Lower())
-            @test e_l == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
+            e_l = boundary_restriction(g_1D,op.eClosure,Lower())
+            @test e_l == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
             @test e_l == BoundaryOperator(g_1D,op.eClosure,Lower())
             @test e_l isa BoundaryOperator{T,Lower} where T
             @test e_l isa TensorMapping{T,0,1} where T
 
-            e_r = BoundaryRestriction(g_1D,op.eClosure,Upper())
-            @test e_r == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}())
+            e_r = boundary_restriction(g_1D,op.eClosure,Upper())
+            @test e_r == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}())
             @test e_r == BoundaryOperator(g_1D,op.eClosure,Upper())
             @test e_r isa BoundaryOperator{T,Upper} where T
             @test e_r isa TensorMapping{T,0,1} where T
         end
 
         @testset "2D" begin
-            e_w = BoundaryRestriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
+            e_w = boundary_restriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
             @test e_w isa InflatedTensorMapping
             @test e_w isa TensorMapping{T,1,2} where T
         end
@@ -865,8 +848,8 @@
 
     @testset "Application" begin
         @testset "1D" begin
-            e_l = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
-            e_r = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}())
+            e_l = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
+            e_r = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}())
 
             v = evalOn(g_1D,x->1+x^2)
             u = fill(3.124)
@@ -877,10 +860,10 @@
         end
 
         @testset "2D" begin
-            e_w = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}())
-            e_e = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}())
-            e_s = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}())
-            e_n = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}())
+            e_w = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}())
+            e_e = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}())
+            e_s = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}())
+            e_n = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}())
 
             v = rand(11, 15)
             u = fill(3.124)
@@ -893,25 +876,25 @@
     end
 end
 
-@testset "NormalDerivative" begin
+@testset "normal_derivative" begin
     g_1D = EquidistantGrid(11, 0.0, 1.0)
     g_2D = EquidistantGrid((11,12), (0.0, 0.0), (1.0,1.0))
-    @testset "Constructors" begin
+    @testset "normal_derivative" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
-            d_l = NormalDerivative(g_1D, op.dClosure, Lower())
-            @test d_l == NormalDerivative(g_1D, op.dClosure, CartesianBoundary{1,Lower}())
+            d_l = normal_derivative(g_1D, op.dClosure, Lower())
+            @test d_l == normal_derivative(g_1D, op.dClosure, CartesianBoundary{1,Lower}())
             @test d_l isa BoundaryOperator{T,Lower} where T
             @test d_l isa TensorMapping{T,0,1} where T
         end
         @testset "2D" begin
             op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-            d_w = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
-            d_n = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
+            d_w = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
+            d_n = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
             Ix = IdentityMapping{Float64}((size(g_2D)[1],))
             Iy = IdentityMapping{Float64}((size(g_2D)[2],))
-            d_l = NormalDerivative(restrict(g_2D,1),op.dClosure,Lower())
-            d_r = NormalDerivative(restrict(g_2D,2),op.dClosure,Upper())
+            d_l = normal_derivative(restrict(g_2D,1),op.dClosure,Lower())
+            d_r = normal_derivative(restrict(g_2D,2),op.dClosure,Upper())
             @test d_w ==  d_l⊗Iy
             @test d_n ==  Ix⊗d_r
             @test d_w isa TensorMapping{T,1,2} where T
@@ -925,10 +908,10 @@
         # TODO: Test for higher order polynomials?
         @testset "2nd order" begin
             op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-            d_w = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
-            d_e = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
-            d_s = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
-            d_n = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
+            d_w = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
+            d_e = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
+            d_s = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
+            d_n = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
 
             @test d_w*v ≈ v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ -v∂x[end,:] atol = 1e-13
@@ -938,10 +921,10 @@
 
         @testset "4th order" begin
             op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-            d_w = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
-            d_e = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
-            d_s = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
-            d_n = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
+            d_w = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
+            d_e = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
+            d_s = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
+            d_n = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
 
             @test d_w*v ≈ v∂x[1,:] atol = 1e-13
             @test d_e*v ≈ -v∂x[end,:] atol = 1e-13