--- a/DiffOps/src/laplace.jl	Fri Jun 28 14:11:57 2019 +0200
+++ b/DiffOps/src/laplace.jl	Fri Jun 28 14:17:13 2019 +0200
@@ -33,7 +33,7 @@
 
 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) = throw(MethodError) # TODO: Implement this
+boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId)
 
 
 """
@@ -100,6 +100,28 @@
     return apply_d_T(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId))
 end
 
+"""
+    BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1}
+
+Implements the boundary operator `q` as a TensorOperator
+"""
+struct BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1}
+    op::D2{T,N,M,K}
+    grid::EquidistantGrid{2}
+    bId::CartesianBoundary
+end
+export BoundaryQuadrature
+
+# TODO: Make this independent of dimension
+function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Int}) where T
+    h = spacing(q.grid)[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)
+
+
 
 
 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end

--- a/DiffOps/test/runtests.jl	Fri Jun 28 14:11:57 2019 +0200
+++ b/DiffOps/test/runtests.jl	Fri Jun 28 14:17:13 2019 +0200
@@ -140,3 +140,48 @@
     @test collect(d_s*g_x) ≈ G_s
     @test collect(d_n*g_x) ≈ G_n
 end
+
+@testset "BoundaryQuadrature" begin
+    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0))
+
+    H_w = BoundaryQuadrature(op, g, CartesianBoundary{1,Lower}())
+    H_e = BoundaryQuadrature(op, g, CartesianBoundary{1,Upper}())
+    H_s = BoundaryQuadrature(op, g, CartesianBoundary{2,Lower}())
+    H_n = BoundaryQuadrature(op, g, CartesianBoundary{2,Upper}())
+
+    v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
+
+    function get_quadrature(N)
+        qc = op.quadratureClosure
+        q = (qc..., ones(N-2*closuresize(op))..., reverse(qc)...)
+        @assert length(q) == N
+        return q
+    end
+
+    v_w = v[1,:]
+    v_e = v[10,:]
+    v_s = v[:,1]
+    v_n = v[:,11]
+
+    q_x = spacing(g)[1].*get_quadrature(10)
+    q_y = spacing(g)[2].*get_quadrature(11)
+
+    @test H_w isa TensorOperator{T,1} where T
+
+    @test domain_size(H_w, (3,)) == (3,)
+    @test domain_size(H_n, (3,)) == (3,)
+
+    @test range_size(H_w, (3,)) == (3,)
+    @test range_size(H_n, (3,)) == (3,)
+
+    @test size(H_w*v_w) == (11,)
+    @test size(H_e*v_e) == (11,)
+    @test size(H_s*v_s) == (10,)
+    @test size(H_n*v_n) == (10,)
+
+    @test collect(H_w*v_w) ≈ q_y.*v_w
+    @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
+end