--- a/test/testSbpOperators.jl	Tue Dec 01 16:41:59 2020 +0100
+++ b/test/testSbpOperators.jl	Tue Dec 01 16:48:34 2020 +0100
@@ -179,21 +179,23 @@
     g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
 
     @testset "Constructors" begin
-        # 1D
-        e_l = BoundaryRestriction{Float64,Lower,4}(op.eClosure,size(g_1D)[1])
-        @test e_l == BoundaryRestriction(g_1D,op.eClosure,Lower())
-        @test e_l == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
-        @test e_l isa TensorMapping{T,0,1} where T
+        @testset "1D" begin
+            e_l = BoundaryRestriction{Float64,Lower,4}(op.eClosure,size(g_1D)[1])
+            @test e_l == BoundaryRestriction(g_1D,op.eClosure,Lower())
+            @test e_l == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
+            @test e_l isa TensorMapping{T,0,1} where T
 
-        e_r = BoundaryRestriction{Float64,Upper,4}(op.eClosure,size(g_1D)[1])
-        @test e_r == BoundaryRestriction(g_1D,op.eClosure,Upper())
-        @test e_r == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}())
-        @test e_r isa TensorMapping{T,0,1} where T
+            e_r = BoundaryRestriction{Float64,Upper,4}(op.eClosure,size(g_1D)[1])
+            @test e_r == BoundaryRestriction(g_1D,op.eClosure,Upper())
+            @test e_r == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}())
+            @test e_r isa TensorMapping{T,0,1} where T
+        end
 
-        # 2D
-        e_w = boundary_restriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
-        @test e_w isa InflatedTensorMapping
-        @test e_w isa TensorMapping{T,1,2} where T
+        @testset "2D" begin
+            e_w = boundary_restriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
+            @test e_w isa InflatedTensorMapping
+            @test e_w isa TensorMapping{T,1,2} where T
+        end
     end
 
     e_l = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
@@ -205,66 +207,72 @@
     e_n = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}())
 
     @testset "Sizes" begin
-        # 1D
-        @test domain_size(e_l) == (11,)
-        @test domain_size(e_r) == (11,)
+        @testset "1D" begin
+            @test domain_size(e_l) == (11,)
+            @test domain_size(e_r) == (11,)
 
-        @test range_size(e_l) == ()
-        @test range_size(e_r) == ()
+            @test range_size(e_l) == ()
+            @test range_size(e_r) == ()
+        end
 
-        # 2D
-        @test domain_size(e_w) == (11,15)
-        @test domain_size(e_e) == (11,15)
-        @test domain_size(e_s) == (11,15)
-        @test domain_size(e_n) == (11,15)
+        @testset "2D" begin
+            @test domain_size(e_w) == (11,15)
+            @test domain_size(e_e) == (11,15)
+            @test domain_size(e_s) == (11,15)
+            @test domain_size(e_n) == (11,15)
 
-        @test range_size(e_w) == (15,)
-        @test range_size(e_e) == (15,)
-        @test range_size(e_s) == (11,)
-        @test range_size(e_n) == (11,)
+            @test range_size(e_w) == (15,)
+            @test range_size(e_e) == (15,)
+            @test range_size(e_s) == (11,)
+            @test range_size(e_n) == (11,)
+        end
     end
 
 
     @testset "Application" begin
-        # 1D
-        v = evalOn(g_1D,x->1+x^2)
-        u = fill(3.124)
-        @test (e_l*v)[] == v[1]
-        @test (e_r*v)[] == v[end]
-        @test (e_r*v)[1] == v[end]
-        @test e_l'*u == [u[]; zeros(10)]
-        @test e_r'*u == [zeros(10); u[]]
+        @testset "1D" begin
+            v = evalOn(g_1D,x->1+x^2)
+            u = fill(3.124)
+            @test (e_l*v)[] == v[1]
+            @test (e_r*v)[] == v[end]
+            @test (e_r*v)[1] == v[end]
+            @test e_l'*u == [u[]; zeros(10)]
+            @test e_r'*u == [zeros(10); u[]]
+        end
 
-        # 2D
-        v = rand(11, 15)
+        @testset "2D" begin
+            v = rand(11, 15)
+            u = fill(3.124)
 
-        @test e_w*v == v[1,:]
-        @test e_e*v == v[end,:]
-        @test e_s*v == v[:,1]
-        @test e_n*v == v[:,end]
+            @test e_w*v == v[1,:]
+            @test e_e*v == v[end,:]
+            @test e_s*v == v[:,1]
+            @test e_n*v == v[:,end]
 
 
-       g_x = rand(11)
-       g_y = rand(15)
+           g_x = rand(11)
+           g_y = rand(15)
 
-       G_w = zeros(Float64, (11,15))
-       G_w[1,:] = g_y
+           G_w = zeros(Float64, (11,15))
+           G_w[1,:] = g_y
 
-       G_e = zeros(Float64, (11,15))
-       G_e[end,:] = g_y
+           G_e = zeros(Float64, (11,15))
+           G_e[end,:] = g_y
 
-       G_s = zeros(Float64, (11,15))
-       G_s[:,1] = g_x
+           G_s = zeros(Float64, (11,15))
+           G_s[:,1] = g_x
+
+           G_n = zeros(Float64, (11,15))
+           G_n[:,end] = g_x
 
-       G_n = zeros(Float64, (11,15))
-       G_n[:,end] = g_x
-
-       @test e_w'*g_y == G_w
-       @test e_e'*g_y == G_e
-       @test e_s'*g_x == G_s
-       @test e_n'*g_x == G_n
+           @test e_w'*g_y == G_w
+           @test e_e'*g_y == G_e
+           @test e_s'*g_x == G_s
+           @test e_n'*g_x == G_n
+       end
 
        @testset "Regions" begin
+           u = fill(3.124)
            @test (e_l'*u)[Index(1,Lower)] == 3.124
            @test (e_l'*u)[Index(2,Lower)] == 0
            @test (e_l'*u)[Index(6,Interior)] == 0