--- a/SbpOperators/src/InverseQuadrature.jl	Thu Sep 24 21:42:54 2020 +0200
+++ b/SbpOperators/src/InverseQuadrature.jl	Thu Sep 24 22:31:48 2020 +0200
@@ -34,9 +34,9 @@
 end
 
 """
-    InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
+    DiagonalNorm{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
 
-Implements the quadrature operator `Hi` of Dim dimension as a TensorMapping
+Implements the inverse diagnoal norm operator `Hi` of Dim dimension as a TensorMapping
 """
 export InverseDiagonalNorm, closuresize
 struct InverseDiagonalNorm{T<:Real,N,M} <: TensorOperator{T,1}

--- a/SbpOperators/src/SbpOperators.jl	Thu Sep 24 21:42:54 2020 +0200
+++ b/SbpOperators/src/SbpOperators.jl	Thu Sep 24 22:31:48 2020 +0200
@@ -9,4 +9,6 @@
 include("readoperator.jl")
 include("laplace/secondderivative.jl")
 include("laplace/laplace.jl")
+include("quadrature/diagonal_inner_product.jl")
+include("quadrature/quadrature.jl")
 end # module

--- a/SbpOperators/src/laplace/laplace.jl	Thu Sep 24 21:42:54 2020 +0200
+++ b/SbpOperators/src/laplace/laplace.jl	Thu Sep 24 22:31:48 2020 +0200
@@ -12,7 +12,7 @@
     #TODO: Write a good constructor
 end
 
-LazyTensors.domain_size(H::Laplace{Dim}, range_size::NTuple{Dim,Integer}) where {Dim} = range_size
+LazyTensors.domain_size(L::Laplace{Dim}, range_size::NTuple{Dim,Integer}) where {Dim} = range_size
 
 function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim}
     error("not implemented")
@@ -24,8 +24,7 @@
     return u
 end
 
-# TODO: Fix dispatch on tuples!
-@inline function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::Index, J::Index) where T
+function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::Index, J::Index) where T
     # 2nd x-derivative
     @inbounds vx = view(v, :, Int(J))
     @inbounds uᵢ = LazyTensors.apply(L.D2[1], vx , I)
@@ -37,9 +36,7 @@
     return uᵢ
 end
 
-@inline function LazyTensors.apply_transpose(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim}
-    return LazyTensors.apply(L, v, I)
-end
+LazyTensors.apply_transpose(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim} = LazyTensors.apply(L, v, I...)
 
 # quadrature(L::Laplace) = Quadrature(L.op, L.grid)
 # inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid)

--- a/SbpOperators/src/laplace/secondderivative.jl	Thu Sep 24 21:42:54 2020 +0200
+++ b/SbpOperators/src/laplace/secondderivative.jl	Thu Sep 24 22:31:48 2020 +0200
@@ -20,32 +20,26 @@
 #      I thought I::Vararg{Index,R} fell back to just Index for R = 1
 
 # Apply for different regions Lower/Interior/Upper or Unknown region
-@inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Lower}) where T
+function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Lower}) where T
     return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.closureStencils[Int(I)], v, Int(I))
 end
 
-@inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Interior}) where T
+function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Interior}) where T
     return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.innerStencil, v, Int(I))
 end
 
-@inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Upper}) where T
+function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Upper}) where T
     N = length(v) # TODO: Use domain_size here instead? N = domain_size(D2,size(v))
     return @inbounds D2.h_inv*D2.h_inv*Int(D2.parity)*apply_stencil_backwards(D2.closureStencils[N-Int(I)+1], v, Int(I))
 end
 
-@inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, index::Index{Unknown}) where T
+function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, index::Index{Unknown}) where T
     N = length(v)  # TODO: Use domain_size here instead?
     r = getregion(Int(index), closuresize(D2), N)
     I = Index(Int(index), r)
     return LazyTensors.apply(D2, v, I)
 end
 
-
-@inline function LazyTensors.apply_transpose(D2::SecondDerivative, v::AbstractVector, I::Index)
-    return LazyTensors.apply(D2, v, I)
-end
+LazyTensors.apply_transpose(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index) where {T} = LazyTensors.apply(D2, v, I)
 
-
-function closuresize(D2::SecondDerivative{T,N,M,K}) where {T<:Real,N,M,K}
-    return M
-end
+closuresize(D2::SecondDerivative{T,N,M,K}) where {T<:Real,N,M,K} = M

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/quadrature/diagonal_inner_product.jl	Thu Sep 24 22:31:48 2020 +0200
@@ -0,0 +1,41 @@
+export DiagonalInnerProduct, closuresize
+"""
+    DiagonalInnerProduct{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
+
+Implements the diagnoal norm operator `H` of Dim dimension as a TensorMapping
+"""
+struct DiagonalInnerProduct{T,M} <: TensorOperator{T,1}
+    h::T # The grid spacing could be included in the stencil already. Preferable?
+    closure::NTuple{M,T}
+    #TODO: Write a nice constructor
+end
+
+LazyTensors.domain_size(H::DiagonalInnerProduct, range_size::NTuple{1,Integer}) = range_size
+
+function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index) where T
+    return @inbounds apply(H, v, I)
+end
+
+function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index{Lower}) where T
+    return @inbounds H.h*H.closure[Int(I)]*v[Int(I)]
+end
+
+function LazyTensors.apply(H::DiagonalInnerProduct{T},v::AbstractVector{T}, I::Index{Upper}) where T
+    N = length(v);
+    return @inbounds H.h*H.closure[N-Int(I)+1]v[Int(I)]
+end
+
+function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index{Interior}) where T
+    return @inbounds H.h*v[Int(I)]
+end
+
+function LazyTensors.apply(H::DiagonalInnerProduct{T},  v::AbstractVector{T}, index::Index{Unknown}) where T
+    N = length(v);
+    r = getregion(Int(index), closuresize(H), N)
+    i = Index(Int(index), r)
+    return LazyTensors.apply(H, v, i)
+end
+
+LazyTensors.apply_transpose(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index) where T = LazyTensors.apply(H,v,I)
+
+closuresize(H::DiagonalInnerProduct{T,M}) where {T,M} = M

--- a/SbpOperators/src/quadrature/quadrature.jl	Thu Sep 24 21:42:54 2020 +0200
+++ b/SbpOperators/src/quadrature/quadrature.jl	Thu Sep 24 22:31:48 2020 +0200
@@ -1,30 +1,27 @@
-# At the moment the grid property is used all over. It could possibly be removed if we implement all the 1D operators as TensorMappings
+export Quadrature
 """
     Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
 
 Implements the quadrature operator `Q` of Dim dimension as a TensorMapping
-The multi-dimensional tensor operator consists of a tuple of 1D DiagonalNorm H
+The multi-dimensional tensor operator consists of a tuple of 1D DiagonalInnerProduct H
 tensor operators.
 """
-export Quadrature
-struct Quadrature{Dim,T<:Real,N,M} <: TensorOperator{T,Dim}
-    H::NTuple{Dim,DiagonalNorm{T,N,M}}
+struct Quadrature{Dim,T<:Real,M} <: TensorOperator{T,Dim}
+    H::NTuple{Dim,DiagonalInnerProduct{T,M}}
 end
 
-LazyTensors.domain_size(Q::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
+LazyTensors.domain_size(Q::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where {Dim} = range_size
 
 function LazyTensors.apply(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim}
     error("not implemented")
 end
 
-LazyTensors.apply_transpose(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {Dim,T} = LazyTensors.apply(Q,v,I)
-
-@inline function LazyTensors.apply(Q::Quadrature{1,T}, v::AbstractVector{T}, I::Index) where T
+function LazyTensors.apply(Q::Quadrature{1,T}, v::AbstractVector{T}, I::Index) where T
     @inbounds q = apply(Q.H[1], v , I)
     return q
 end
 
-@inline function LazyTensors.apply(Q::Quadrature{2,T}, v::AbstractArray{T,2}, I::Index, J::Index) where T
+function LazyTensors.apply(Q::Quadrature{2,T}, v::AbstractArray{T,2}, I::Index, J::Index) where T
     # Quadrature in x direction
     @inbounds vx = view(v, :, Int(J))
     @inbounds qx = apply(Q.H[1], vx , I)
@@ -34,43 +31,4 @@
     return qx*qy
 end
 
-"""
-    DiagonalNorm{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
-
-Implements the diagnoal norm operator `H` of Dim dimension as a TensorMapping
-"""
-export DiagonalNorm, closuresize, LazyTensors.apply
-struct DiagonalNorm{T<:Real,N,M} <: TensorOperator{T,1}
-    h::T # The grid spacing could be included in the stencil already. Preferable?
-    closure::NTuple{M,T}
-    #TODO: Write a nice constructor
-end
-
-@inline function LazyTensors.apply(H::DiagonalNorm{T}, v::AbstractVector{T}, I::Index) where T
-    return @inbounds apply(H, v, I)
-end
-
-LazyTensors.apply_transpose(H::Quadrature{Dim,T}, v::AbstractArray{T,2}, I::Index) where T = LazyTensors.apply(H,v,I)
-
-@inline LazyTensors.apply(H::DiagonalNorm, v::AbstractVector{T}, I::Index{Lower}) where T
-    return @inbounds H.h*H.closure[Int(I)]*v[Int(I)]
-end
-@inline LazyTensors.apply(H::DiagonalNorm,v::AbstractVector{T}, I::Index{Upper}) where T
-    N = length(v);
-    return @inbounds H.h*H.closure[N-Int(I)+1]v[Int(I)]
-end
-
-@inline LazyTensors.apply(H::DiagonalNorm, v::AbstractVector{T}, I::Index{Interior}) where T
-    return @inbounds H.h*v[Int(I)]
-end
-
-function LazyTensors.apply(H::DiagonalNorm,  v::AbstractVector{T}, index::Index{Unknown}) where T
-    N = length(v);
-    r = getregion(Int(index), closuresize(H), N)
-    i = Index(Int(index), r)
-    return LazyTensors.apply(H, v, i)
-end
-
-function closuresize(H::DiagonalNorm{T<:Real,N,M}) where {T,N,M}
-    return M
-end
+LazyTensors.apply_transpose(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {Dim,T} = LazyTensors.apply(Q,v,I...)

--- a/SbpOperators/test/runtests.jl	Thu Sep 24 21:42:54 2020 +0200
+++ b/SbpOperators/test/runtests.jl	Thu Sep 24 22:31:48 2020 +0200
@@ -119,20 +119,39 @@
     @test sqrt(prod(h)*sum(collect(e4.^2))) <= accuracytol
     @test sqrt(prod(h)*sum(collect(e5.^2))) <= accuracytol
 end
-#
-# @testset "Quadrature" begin
-#     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
-#     Lx = 2.3
-#     Ly = 5.2
-#     g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
-#     H = Quadrature(op,g)
-#     v = ones(Float64, size(g))
-#
-#     @test H isa TensorOperator{T,2} where T
-#     @test H' isa TensorMapping{T,2,2} where T
-#     @test sum(collect(H*v)) ≈ (Lx*Ly)
-#     @test collect(H*v) == collect(H'*v)
-# end
+
+@testset "DiagonalInnerProduct" begin
+    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    L = 2.3
+    g = EquidistantGrid((77,), (0.0,), (L,))
+    h = spacing(g)
+    H = DiagonalInnerProduct(h[1],op.quadratureClosure)
+    v = ones(Float64, size(g))
+
+    @test H isa TensorOperator{T,1} where T
+    @test H' isa TensorMapping{T,1,1} where T
+    @test sum(collect(H*v)) ≈ L
+    @test collect(H*v) == collect(H'*v)
+end
+
+@testset "Quadrature" begin
+    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    Lx = 2.3
+    Ly = 5.2
+    g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
+
+    h = spacing(g)
+    Hx = DiagonalInnerProduct(h[1],op.quadratureClosure);
+    Hy = DiagonalInnerProduct(h[2],op.quadratureClosure);
+    Q = Quadrature((Hx,Hy))
+
+    v = ones(Float64, size(g))
+
+    @test Q isa TensorOperator{T,2} where T
+    @test Q' isa TensorMapping{T,2,2} where T
+    @test sum(collect(Q*v)) ≈ (Lx*Ly)
+    @test collect(Q*v) == collect(Q'*v)
+end
 #
 # @testset "InverseQuadrature" begin
 #     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")