--- a/DiffOps/src/laplace.jl	Thu Jun 27 12:49:58 2019 +0200
+++ b/DiffOps/src/laplace.jl	Thu Jun 27 13:36:04 2019 +0200
@@ -1,69 +1,3 @@
-"""
-    NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
-
-Implements the boundary operator `d` as a TensorMapping
-"""
-struct NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
-	op::D2{T,N,M,K}
-	grid::EquidistantGrid
-	bId::CartesianBoundary
-end
-export NormalDerivative
-
-# 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{T}, domain_size::NTuple{1,Integer}) where T = size(e.grid)
-LazyTensors.domain_size(e::NormalDerivative{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
-
-# Not correct abstraction level
-# TODO: Not type stable D:<
-function LazyTensors.apply(d::NormalDerivative, v::AbstractArray, I::NTuple{2,Int})
-	i = I[dim(d.bId)]
-	j = I[3-dim(d.bId)]
-	N_i = size(d.grid)[dim(d.bId)]
-
-    h_inv = d.grid.inverse_spacing[dim(d.bId)]
-    return apply_d(d.op, h_inv, v[j], N_i, i, region(d.bId))
-end
-
-function LazyTensors.apply_transpose(d::NormalDerivative, v::AbstractArray, I::NTuple{1,Int})
-    u = selectdim(v,3-dim(d.bId),I[1])
-    return apply_d_T(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId))
-end
-
-
-"""
-    BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
-
-Implements the boundary operator `e` as a TensorMapping
-"""
-struct BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
-	op::D2{T,N,M,K}
-	grid::EquidistantGrid
-	bId::CartesianBoundary
-end
-export BoundaryValue
-
-# 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)
-LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
-
-function LazyTensors.apply(e::BoundaryValue, v::AbstractArray, I::NTuple{2,Int})
-	i = I[dim(e.bId)]
-	j = I[3-dim(e.bId)]
-	N_i = size(e.grid)[dim(e.bId)]
-
-    return apply_e(e.op, v[j], N_i, i, region(e.bId))
-end
-
-function LazyTensors.apply_transpose(e::BoundaryValue, v::AbstractArray, I::NTuple{1,Int})
-	u = selectdim(v,3-dim(e.bId),I[1])
-	return apply_e_T(e.op, u, region(e.bId))
-end
-
-
-
 struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim}
     grid::EquidistantGrid{Dim,T}
     a::T
@@ -100,6 +34,72 @@
 end
 
 
+
+"""
+    BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
+
+    Implements the boundary operator `e` as a TensorMapping
+"""
+struct BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
+    op::D2{T,N,M,K}
+    grid::EquidistantGrid
+    bId::CartesianBoundary
+end
+export BoundaryValue
+
+# 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)
+LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
+
+function LazyTensors.apply(e::BoundaryValue, v::AbstractArray, I::NTuple{2,Int})
+    i = I[dim(e.bId)]
+    j = I[3-dim(e.bId)]
+    N_i = size(e.grid)[dim(e.bId)]
+    return apply_e(e.op, v[j], N_i, i, region(e.bId))
+end
+
+function LazyTensors.apply_transpose(e::BoundaryValue, v::AbstractArray, I::NTuple{1,Int})
+    u = selectdim(v,3-dim(e.bId),I[1])
+    return apply_e_T(e.op, u, region(e.bId))
+end
+
+
+
+"""
+    NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
+
+    Implements the boundary operator `d` as a TensorMapping
+"""
+struct NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
+    op::D2{T,N,M,K}
+    grid::EquidistantGrid
+    bId::CartesianBoundary
+end
+export NormalDerivative
+
+# 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{T}, domain_size::NTuple{1,Integer}) where T = size(e.grid)
+LazyTensors.domain_size(e::NormalDerivative{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
+
+# Not correct abstraction level
+# TODO: Not type stable D:<
+function LazyTensors.apply(d::NormalDerivative, v::AbstractArray, I::NTuple{2,Int})
+    i = I[dim(d.bId)]
+    j = I[3-dim(d.bId)]
+    N_i = size(d.grid)[dim(d.bId)]
+    h_inv = d.grid.inverse_spacing[dim(d.bId)]
+    return apply_d(d.op, h_inv, v[j], N_i, i, region(d.bId))
+end
+
+function LazyTensors.apply_transpose(d::NormalDerivative, v::AbstractArray, I::NTuple{1,Int})
+    u = selectdim(v,3-dim(d.bId),I[1])
+    return apply_d_T(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId))
+end
+
+
+
 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end
 
 function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid}
@@ -127,7 +127,7 @@
 end
 
 # function apply(s::MyWaveEq{D},  v::AbstractArray{T,D}, i::CartesianIndex{D}) where D
-# 	return apply(s.L, v, i) +
+    #   return apply(s.L, v, i) +
 # 		sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau),  v, s.g_w, i) +
 # 		sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau),  v, s.g_e, i) +
 # 		sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau),  v, s.g_s, i) +

--- a/Manifest.toml	Thu Jun 27 12:49:58 2019 +0200
+++ b/Manifest.toml	Thu Jun 27 13:36:04 2019 +0200
@@ -4,7 +4,7 @@
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
 [[DiffOps]]
-deps = ["Grids", "RegionIndices", "SbpOperators", "TiledIteration"]
+deps = ["Grids", "LazyTensors", "RegionIndices", "SbpOperators", "TiledIteration"]
 path = "DiffOps"
 uuid = "39474f48-97ec-11e9-01fc-6ddcbe5918df"
 version = "0.1.0"