--- 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