--- a/DiffOps/src/laplace.jl	Tue Nov 26 08:19:22 2019 -0800
+++ b/DiffOps/src/laplace.jl	Tue Nov 26 08:28:26 2019 -0800
@@ -31,10 +31,65 @@
     apply(L, v, I)
 end
 
+quadrature(L::Laplace) = Quadrature(L.op, L.grid)
+inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid)
 boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId)
 normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId)
 boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId)
 
+"""
+    Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
+
+Implements the quadrature operator `H` of Dim dimension as a TensorMapping
+"""
+struct Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
+    op::D2{T,N,M,K}
+    grid::EquidistantGrid{Dim,T}
+end
+export Quadrature
+
+LazyTensors.range_size(H::Quadrature{2}, domain_size::NTuple{2,Integer}) where T = size(H.grid)
+LazyTensors.domain_size(H::Quadrature{2}, range_size::NTuple{2,Integer}) where T = size(H.grid)
+
+# TODO: Dispatch on Tuple{Index{R1},Index{R2}}?
+@inline function LazyTensors.apply(H::Quadrature{2}, v::AbstractArray{T,2} where T, I::NTuple{2,Integer})
+    I = CartesianIndex(I);
+    N = size(H.grid)
+    # Quadrature in x direction
+    @inbounds q = apply_quadrature(H.op, H.grid.spacing[1], v[I] , I[1], N[1])
+    # Quadrature in y-direction
+    @inbounds q = apply_quadrature(H.op, H.grid.spacing[2], q, I[2], N[2])
+    return q
+end
+
+LazyTensors.apply_transpose(H::Quadrature{2}, v::AbstractArray{T,2} where T, I::NTuple{2,Integer}) = LazyTensors.apply(H,v,I)
+
+"""
+    InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
+
+Implements the inverse quadrature operator `inv(H)` of Dim dimension as a TensorMapping
+"""
+struct InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
+    op::D2{T,N,M,K}
+    grid::EquidistantGrid{Dim,T}
+end
+export InverseQuadrature
+
+LazyTensors.range_size(H_inv::InverseQuadrature{2}, domain_size::NTuple{2,Integer}) where T = size(H_inv.grid)
+LazyTensors.domain_size(H_inv::InverseQuadrature{2}, range_size::NTuple{2,Integer}) where T = size(H_inv.grid)
+
+# TODO: Dispatch on Tuple{Index{R1},Index{R2}}?
+@inline function LazyTensors.apply(H_inv::InverseQuadrature{2}, v::AbstractArray{T,2} where T, I::NTuple{2,Integer})
+    I = CartesianIndex(I);
+    N = size(H_inv.grid)
+    # Inverse quadrature in x direction
+    @inbounds q_inv = apply_inverse_quadrature(H_inv.op, H_inv.grid.inverse_spacing[1], v[I] , I[1], N[1])
+    # Inverse quadrature in y-direction
+    @inbounds q_inv = apply_inverse_quadrature(H_inv.op, H_inv.grid.inverse_spacing[2], q_inv, I[2], N[2])
+    return q_inv
+end
+
+LazyTensors.apply_transpose(H_inv::InverseQuadrature{2}, v::AbstractArray{T,2} where T, I::NTuple{2,Integer}) = LazyTensors.apply(H_inv,v,I)
 
 """
     BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
@@ -66,8 +121,6 @@
     return apply_e_T(e.op, u, region(e.bId))
 end
 
-
-
 """
     NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
 
@@ -113,13 +166,14 @@
 export BoundaryQuadrature
 
 # TODO: Make this independent of dimension
+# TODO: Dispatch directly on Index{R}?
 function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Int}) where T
-    h = spacing(q.grid)[3-dim(q.bId)]
+    h = q.grid.spacing[3-dim(q.bId)]
     N = size(v)
     return apply_quadrature(q.op, h, v[I[1]], I[1], N[1])
 end
 
-LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Int}) where T = apply(q,v,I)
+LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Int}) where T = LazyTensors.apply(q,v,I)
 
 
 

--- a/DiffOps/test/runtests.jl	Tue Nov 26 08:19:22 2019 -0800
+++ b/DiffOps/test/runtests.jl	Tue Nov 26 08:28:26 2019 -0800
@@ -5,6 +5,35 @@
 using RegionIndices
 using LazyTensors
 
+@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 TensorMapping{T,2,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 "InverseQuadrature" begin
+    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    Lx = 7.3
+    Ly = 8.2
+    g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
+    H = Quadrature(op,g)
+    Hinv = InverseQuadrature(op,g)
+    v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
+
+    @test Hinv isa TensorMapping{T,2,2} where T
+    @test Hinv' isa TensorMapping{T,2,2} where T
+    @test collect(Hinv*H*v)  ≈ v
+    @test collect(Hinv*v) == collect(Hinv'*v)
+end
+
 @testset "BoundaryValue" begin
     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
     g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
@@ -184,4 +213,9 @@
     @test collect(H_e*v_e) ≈ q_y.*v_e
     @test collect(H_s*v_s) ≈ q_x.*v_s
     @test collect(H_n*v_n) ≈ q_x.*v_n
+
+    @test collect(H_w'*v_w) == collect(H_w'*v_w)
+    @test collect(H_e'*v_e) == collect(H_e'*v_e)
+    @test collect(H_s'*v_s) == collect(H_s'*v_s)
+    @test collect(H_n'*v_n) == collect(H_n'*v_n)
 end