--- a/src/SbpOperators/SbpOperators.jl	Sat Feb 13 16:05:02 2021 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Mon Feb 15 11:13:12 2021 +0100
@@ -10,8 +10,8 @@
 include("volumeops/volume_operator.jl")
 include("volumeops/derivatives/secondderivative.jl")
 include("volumeops/laplace/laplace.jl")
-include("volumeops/quadratures/quadrature.jl")
-include("volumeops/quadratures/inverse_quadrature.jl")
+include("volumeops/inner_products/inner_product.jl")
+include("volumeops/inner_products/inverse_inner_product.jl")
 include("boundaryops/boundary_operator.jl")
 include("boundaryops/boundary_restriction.jl")
 include("boundaryops/normal_derivative.jl")

--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Sat Feb 13 16:05:02 2021 +0100
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Mon Feb 15 11:13:12 2021 +0100
@@ -1,6 +1,6 @@
 """
-    BoundaryRestriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
-    BoundaryRestriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
+    boundary_restriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
+    boundary_restriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
 
 Creates the boundary restriction operator `e` as a `TensorMapping`
 
@@ -9,7 +9,7 @@
 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.
 """
-BoundaryRestriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary) = SbpOperators.boundary_operator(grid, closure_stencil, boundary)
-BoundaryRestriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = BoundaryRestriction(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
+boundary_restriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary) = SbpOperators.boundary_operator(grid, closure_stencil, boundary)
+boundary_restriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = boundary_restriction(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
 
-export BoundaryRestriction
+export boundary_restriction

--- a/src/SbpOperators/boundaryops/normal_derivative.jl	Sat Feb 13 16:05:02 2021 +0100
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Mon Feb 15 11:13:12 2021 +0100
@@ -1,6 +1,6 @@
 """
-    NormalDerivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
-    NormalDerivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
+    normal_derivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
+    normal_derivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
 
 Creates the normal derivative boundary operator `d` as a `TensorMapping`
 
@@ -9,10 +9,10 @@
 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.
 """
-function NormalDerivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
+function normal_derivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
     direction = dim(boundary)
     h_inv = inverse_spacing(grid)[direction]
     return SbpOperators.boundary_operator(grid, scale(closure_stencil,h_inv), boundary)
 end
-NormalDerivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = NormalDerivative(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
-export NormalDerivative
+normal_derivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = normal_derivative(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
+export normal_derivative

--- a/src/SbpOperators/volumeops/derivatives/secondderivative.jl	Sat Feb 13 16:05:02 2021 +0100
+++ b/src/SbpOperators/volumeops/derivatives/secondderivative.jl	Mon Feb 15 11:13:12 2021 +0100
@@ -1,6 +1,6 @@
 """
-    SecondDerivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction)
-    SecondDerivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils)
+    second_derivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction)
+    second_derivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils)
 
 Creates the second-derivative operator `D2` as a `TensorMapping`
 
@@ -12,9 +12,9 @@
 one-dimensional operator with the `IdentityMapping`s in orthogonal coordinate dirrections.
 Also see the documentation of `SbpOperators.volume_operator(...)` for more details.
 """
-function SecondDerivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction) where Dim
+function second_derivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction) where Dim
     h_inv = inverse_spacing(grid)[direction]
     return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even, direction)
 end
-SecondDerivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils) = SecondDerivative(grid,inner_stencil,closure_stencils,1)
-export SecondDerivative
+second_derivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils) = second_derivative(grid,inner_stencil,closure_stencils,1)
+export second_derivative

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Mon Feb 15 11:13:12 2021 +0100
@@ -0,0 +1,29 @@
+"""
+    inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil)
+
+Creates the discrete inner product operator `H` as a `TensorMapping` on an equidistant
+grid, defined as `(u,v)  = u'Hv` for grid functions `u,v`.
+
+`inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil)` creates
+`H` on `grid` the using a set of stencils `closure_stencils` for the points in
+the closure regions and the stencil and `inner_stencil` in the interior. If
+`inner_stencil` is omitted a central interior stencil with weight 1 is used.
+
+On a 1-dimensional `grid`, `H` is a `VolumeOperator`. On a N-dimensional
+`grid`, `H` is the outer product of the 1-dimensional inner product operators in
+each coordinate direction. Also see the documentation of
+`SbpOperators.volume_operator(...)` for more details. On a 0-dimensional `grid`,
+`H` is a 0-dimensional `IdentityMapping`.
+"""
+function inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil = CenteredStencil(one(eltype(grid))))
+    h = spacing(grid)
+    H = SbpOperators.volume_operator(grid, scale(inner_stencil,h[1]), scale.(closure_stencils,h[1]), even, 1)
+    for i ∈ 2:dimension(grid)
+        Hᵢ = SbpOperators.volume_operator(grid, scale(inner_stencil,h[i]), scale.(closure_stencils,h[i]), even, i)
+        H = H∘Hᵢ
+    end
+    return H
+end
+export inner_product
+
+inner_product(grid::EquidistantGrid{0}, closure_stencils, inner_stencil) = IdentityMapping{eltype(grid)}()

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Mon Feb 15 11:13:12 2021 +0100
@@ -0,0 +1,43 @@
+"""
+    inverse_inner_product(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)
+    inverse_inner_product(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}})
+
+Creates the inverse inner product operator `H⁻¹` as a `TensorMapping` on an
+equidistant grid. `H⁻¹` is defined implicitly by `H⁻¹∘H = I`, where
+`H` is the corresponding inner product operator and `I` is the `IdentityMapping`.
+
+`inverse_inner_product(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)`
+constructs `H⁻¹` using a set of stencils `inv_closure_stencils` for the points
+in the closure regions and the stencil `inv_inner_stencil` in the interior. If
+`inv_closure_stencils` is omitted, a central interior stencil with weight 1 is used.
+
+`inverse_inner_product(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}})`
+constructs a diagonal inverse inner product operator where `closure_stencils` are the
+closure stencils of `H` (not `H⁻¹`!).
+
+On a 1-dimensional `grid`, `H⁻¹` is a `VolumeOperator`. On a N-dimensional
+`grid`, `H⁻¹` is the outer product of the 1-dimensional inverse inner product
+operators in each coordinate direction. Also see the documentation of
+`SbpOperators.volume_operator(...)` for more details. On a 0-dimensional `grid`,
+`H⁻¹` is a 0-dimensional `IdentityMapping`.
+"""
+function inverse_inner_product(grid::EquidistantGrid, inv_closure_stencils, inv_inner_stencil = CenteredStencil(one(eltype(grid))))
+    h⁻¹ = inverse_spacing(grid)
+    H⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[1]),scale.(inv_closure_stencils,h⁻¹[1]),even,1)
+    for i ∈ 2:dimension(grid)
+        Hᵢ⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[i]),scale.(inv_closure_stencils,h⁻¹[i]),even,i)
+        H⁻¹ = H⁻¹∘Hᵢ⁻¹
+    end
+    return H⁻¹
+end
+export inverse_inner_product
+
+inverse_inner_product(grid::EquidistantGrid{0}, inv_closure_stencils, inv_inner_stencil) = IdentityMapping{eltype(grid)}()
+
+function inverse_inner_product(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {M,T}
+     inv_closure_stencils = reciprocal_stencil.(closure_stencils)
+     inv_inner_stencil = CenteredStencil(one(T))
+     return inverse_inner_product(grid, inv_closure_stencils, inv_inner_stencil)
+end
+
+reciprocal_stencil(s::Stencil{T}) where T = Stencil(s.range,one(T)./s.weights)

--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Sat Feb 13 16:05:02 2021 +0100
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Mon Feb 15 11:13:12 2021 +0100
@@ -1,5 +1,5 @@
 """
-    Laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils)
+    laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils)
 
 Creates the Laplace operator operator `Δ` as a `TensorMapping`
 
@@ -7,14 +7,15 @@
 the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
 for the points in the closure regions.
 
-On a one-dimensional `grid`, `Δ` is a `SecondDerivative`. On a multi-dimensional `grid`, `Δ` is the sum of
-multi-dimensional `SecondDerivative`s where the sum is carried out lazily.
+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.
 """
-function Laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils) where Dim
-    Δ = SecondDerivative(grid, inner_stencil, closure_stencils, 1)
+function laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils) where Dim
+    Δ = second_derivative(grid, inner_stencil, closure_stencils, 1)
     for d = 2:Dim
-        Δ += SecondDerivative(grid, inner_stencil, closure_stencils, d)
+        Δ += second_derivative(grid, inner_stencil, closure_stencils, d)
     end
     return Δ
 end
-export Laplace
+export laplace

--- a/src/SbpOperators/volumeops/quadratures/inverse_quadrature.jl	Sat Feb 13 16:05:02 2021 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,41 +0,0 @@
-
-"""
-    InverseQuadrature(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)
-
-Creates the inverse `H⁻¹` of the quadrature operator as a `TensorMapping`
-
-The inverse quadrature approximates the integral operator on the grid using
-`inv_inner_stencil` in the interior and a set of stencils `inv_closure_stencils`
-for the points in the closure regions.
-
-On a one-dimensional `grid`, `H⁻¹` is a `VolumeOperator`. On a multi-dimensional
-`grid`, `H` is the outer product of the 1-dimensional inverse quadrature operators in
-each coordinate direction. Also see the documentation of
-`SbpOperators.volume_operator(...)` for more details.
-"""
-function InverseQuadrature(grid::EquidistantGrid{Dim}, inv_inner_stencil, inv_closure_stencils) where Dim
-    h⁻¹ = inverse_spacing(grid)
-    H⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[1]),scale.(inv_closure_stencils,h⁻¹[1]),even,1)
-    for i ∈ 2:Dim
-        Hᵢ⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[i]),scale.(inv_closure_stencils,h⁻¹[i]),even,i)
-        H⁻¹ = H⁻¹∘Hᵢ⁻¹
-    end
-    return H⁻¹
-end
-export InverseQuadrature
-
-"""
-    InverseDiagonalQuadrature(grid::EquidistantGrid, closure_stencils)
-
-Creates the inverse of the diagonal quadrature operator defined by the inner stencil
-1/h and a set of 1-element closure stencils in `closure_stencils`. Note that
-the closure stencils are those of the quadrature operator (and not the inverse).
-"""
-function InverseDiagonalQuadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {T,M}
-    inv_inner_stencil = Stencil(one(T), center=1)
-    inv_closure_stencils = reciprocal_stencil.(closure_stencils)
-    return InverseQuadrature(grid, inv_inner_stencil, inv_closure_stencils)
-end
-export InverseDiagonalQuadrature
-
-reciprocal_stencil(s::Stencil{T}) where T = Stencil(s.range,one(T)./s.weights)

--- a/src/SbpOperators/volumeops/quadratures/quadrature.jl	Sat Feb 13 16:05:02 2021 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,29 +0,0 @@
-"""
-    quadrature(grid::EquidistantGrid, closure_stencils, inner_stencil)
-    quadrature(grid::EquidistantGrid, closure_stencils)
-
-Creates the quadrature operator `H` as a `TensorMapping`
-
-`H` approximiates the integral operator on `grid` the using the stencil
-`inner_stencil` in the interior and a set of stencils `closure_stencils`
-for the points in the closure regions. If `inner_stencil` is omitted a central
-interior stencil with weight 1 is used.
-
-On a one-dimensional `grid`, `H` is a `VolumeOperator`. On a multi-dimensional
-`grid`, `H` is the outer product of the 1-dimensional quadrature operators in
-each coordinate direction. Also see the documentation of
-`SbpOperators.volume_operator(...)` for more details. On a 0-dimensional `grid`,
-`H` is a 0-dimensional `IdentityMapping`.
-"""
-function quadrature(grid::EquidistantGrid, closure_stencils, inner_stencil = CenteredStencil(one(eltype(grid))))
-    h = spacing(grid)
-    H = SbpOperators.volume_operator(grid, scale(inner_stencil,h[1]), scale.(closure_stencils,h[1]), even, 1)
-    for i ∈ 2:dimension(grid)
-        Hᵢ = SbpOperators.volume_operator(grid, scale(inner_stencil,h[i]), scale.(closure_stencils,h[i]), even, i)
-        H = H∘Hᵢ
-    end
-    return H
-end
-export quadrature
-
-quadrature(grid::EquidistantGrid{0}, closure_stencils, inner_stencil) = IdentityMapping{eltype(grid)}()

--- a/test/testSbpOperators.jl	Sat Feb 13 16:05:02 2021 +0100
+++ b/test/testSbpOperators.jl	Mon Feb 15 11:13:12 2021 +0100
@@ -241,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
@@ -272,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
@@ -283,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
@@ -309,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
@@ -320,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
@@ -339,16 +339,16 @@
     @testset "Constructors" begin
         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)
+            L = laplace(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)
+            L = laplace(g_3D, op.innerStencil, op.closureStencils)
             @test L isa TensorMapping{T,3,3} where T
-            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)
+            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
         end
     end
@@ -370,7 +370,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)
-            L = Laplace(g_3D,op.innerStencil,op.closureStencils)
+            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
@@ -381,7 +381,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)
-            L = Laplace(g_3D,op.innerStencil,op.closureStencils)
+            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
@@ -393,7 +393,7 @@
     end
 end
 
-@testset "Quadrature diagonal" begin
+@testset "Diagonal-stencil inner_product" begin
     Lx = π/2.
     Ly = Float64(π)
     Lz = 1.
@@ -401,23 +401,23 @@
     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 = quadrature(EquidistantGrid{Float64}(),op.quadratureClosure)
+            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)
+            H = inner_product(g_1D,op.quadratureClosure)
             inner_stencil = CenteredStencil(1.)
-            @test H == quadrature(g_1D,op.quadratureClosure,inner_stencil)
+            @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
@@ -426,12 +426,12 @@
     @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
@@ -448,7 +448,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
@@ -457,7 +457,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
@@ -471,13 +471,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
@@ -485,27 +485,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);
+            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
@@ -514,12 +519,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
@@ -531,15 +536,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
@@ -549,15 +554,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
@@ -578,7 +583,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
@@ -709,28 +714,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
@@ -738,8 +743,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)
@@ -750,10 +755,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)
@@ -766,25 +771,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
@@ -798,10 +803,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
@@ -811,10 +816,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