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