--- a/DiffOps/src/laplace.jl	Thu Jan 09 13:38:06 2020 +0100
+++ b/DiffOps/src/laplace.jl	Thu Jun 18 22:07:10 2020 +0200
@@ -1,11 +1,13 @@
 struct Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
-    grid::EquidistantGrid{Dim,T}
+    grid::EquidistantGrid{Dim,T}  # TODO: Should this be here? Should probably be possible to applay a Laplace object to any size grid function
     a::T # TODO: Better name?
     op::D2{T,N,M,K}
 end
 export Laplace
 
-LazyTensors.domain_size(H::Laplace{Dim}, range_size::NTuple{Dim,Integer}) where Dim = size(L.grid)
+# At the moment the grid property is used all over. It could possibly be removed if we implement all the 1D operators as TensorMappings
+
+LazyTensors.domain_size(H::Laplace{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
 
 function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim}
     error("not implemented")
@@ -46,7 +48,7 @@
 end
 export Quadrature
 
-LazyTensors.domain_size(H::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = size(H.grid)
+LazyTensors.domain_size(H::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
 
 @inline function LazyTensors.apply(H::Quadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
     N = size(H.grid)
@@ -71,7 +73,7 @@
 end
 export InverseQuadrature
 
-LazyTensors.domain_size(H_inv::InverseQuadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = size(H_inv.grid)
+LazyTensors.domain_size(H_inv::InverseQuadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
 
 @inline function LazyTensors.apply(H_inv::InverseQuadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
     N = size(H_inv.grid)
@@ -98,8 +100,15 @@
 
 # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
 # Can we give special treatment to TensorMappings that go to a higher dim?
-LazyTensors.range_size(e::BoundaryValue{T}, domain_size::NTuple{1,Integer}) where T = size(e.grid)
+function LazyTensors.range_size(e::BoundaryValue{T}, domain_size::NTuple{1,Integer}) where T
+    if dim(e.bId) == 1
+        return (missing, domain_size[1])
+    elseif dim(e.bId) == 2
+        return (domain_size[1], missing)
+    end
+end
 LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
+# TODO: Make a nicer solution for 3-dim(e.bId)
 
 # TODO: Make this independent of dimension
 function LazyTensors.apply(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T
@@ -128,7 +137,13 @@
 
 # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
 # Can we give special treatment to TensorMappings that go to a higher dim?
-LazyTensors.range_size(e::NormalDerivative, domain_size::NTuple{1,Integer}) = size(e.grid)
+function LazyTensors.range_size(e::NormalDerivative, domain_size::NTuple{1,Integer})
+    if dim(e.bId) == 1
+        return (missing, domain_size[1])
+    elseif dim(e.bId) == 2
+        return (domain_size[1], missing)
+    end
+end
 LazyTensors.domain_size(e::NormalDerivative, range_size::NTuple{2,Integer}) = (range_size[3-dim(e.bId)],)
 
 # TODO: Not type stable D:<

--- a/DiffOps/test/runtests.jl	Thu Jan 09 13:38:06 2020 +0100
+++ b/DiffOps/test/runtests.jl	Thu Jun 18 22:07:10 2020 +0200
@@ -129,15 +129,16 @@
     G_n = zeros(Float64, (4,5))
     G_n[:,5] = g_x
 
-    @test size(e_w*g_y) == (4,5)
-    @test size(e_e*g_y) == (4,5)
-    @test size(e_s*g_x) == (4,5)
-    @test size(e_n*g_x) == (4,5)
+    @test size(e_w*g_y) === (missing,5)
+    @test size(e_e*g_y) === (missing,5)
+    @test size(e_s*g_x) === (4,missing)
+    @test size(e_n*g_x) === (4,missing)
 
-    @test collect(e_w*g_y) == G_w
-    @test collect(e_e*g_y) == G_e
-    @test collect(e_s*g_x) == G_s
-    @test collect(e_n*g_x) == G_n
+    # These tests should be moved to where they are possible (i.e we know what the grid should be)
+    @test_broken collect(e_w*g_y) == G_w
+    @test_broken collect(e_e*g_y) == G_e
+    @test_broken collect(e_s*g_x) == G_s
+    @test_broken collect(e_n*g_x) == G_n
 end
 
 @testset "NormalDerivative" begin
@@ -205,15 +206,16 @@
     G_n = prod_matrix(g_x, d_y_u)
 
 
-    @test size(d_w*g_y) == (5,6)
-    @test size(d_e*g_y) == (5,6)
-    @test size(d_s*g_x) == (5,6)
-    @test size(d_n*g_x) == (5,6)
+    @test size(d_w*g_y) === (missing,6)
+    @test size(d_e*g_y) === (missing,6)
+    @test size(d_s*g_x) === (5,missing)
+    @test size(d_n*g_x) === (5,missing)
 
-    @test collect(d_w*g_y) ≈ G_w
-    @test collect(d_e*g_y) ≈ G_e
-    @test collect(d_s*g_x) ≈ G_s
-    @test collect(d_n*g_x) ≈ G_n
+    # These tests should be moved to where they are possible (i.e we know what the grid should be)
+    @test_broken collect(d_w*g_y) ≈ G_w
+    @test_broken collect(d_e*g_y) ≈ G_e
+    @test_broken collect(d_s*g_x) ≈ G_s
+    @test_broken collect(d_n*g_x) ≈ G_n
 end
 
 @testset "BoundaryQuadrature" begin

--- a/TODO.txt	Thu Jan 09 13:38:06 2020 +0100
+++ b/TODO.txt	Thu Jun 18 22:07:10 2020 +0200
@@ -21,3 +21,13 @@
 Just nu agerar apply_normal_derivative, apply_boundary_value på inte på v som
 en vektor, utan randvärdet plockas ut utanför. Känns inte konsistent med övrig
 design
+
+## TODO 2020-06-18
+ * Remove grid as a property of the Laplace tensor operator
+ * Add 1D operators (D1, D2, e, d ... ) as TensorOperators
+ * Move Laplace tensor operator to different package
+ * Add new Laplace opertor to DiffOps, probably named WaveEqOp(?!!?)
+ * Decide: Should there be some kind of collection struct for SBP operators (as TensorOperators), providing easy access to all parts (D2, e, d ,
+ H.. H_gamma etc.)
+ * Is "missing" a good value for unknown dimension sizes (of e*g for example)
+ * Formalize how range_size() and domain_size() are supposed to work in TensorMappings where dim(domain) != dim(range) (add tests or document)