--- a/DiffOps/src/laplace.jl	Thu Jun 27 14:33:39 2019 +0200
+++ b/DiffOps/src/laplace.jl	Thu Jun 27 14:37:21 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
@@ -102,6 +36,71 @@
 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)],)
+
+# 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 +126,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 14:33:39 2019 +0200
+++ b/Manifest.toml	Thu Jun 27 14:37:21 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"

--- a/SbpOperators/src/SbpOperators.jl	Thu Jun 27 14:33:39 2019 +0200
+++ b/SbpOperators/src/SbpOperators.jl	Thu Jun 27 14:37:21 2019 +0200
@@ -3,205 +3,8 @@
 using RegionIndices
 
 include("stencil.jl")
-
-export D2, closureSize, apply, readOperator, apply_e, apply_d, apply_e_T, apply_d_T
-
-abstract type ConstantStencilOperator end
-
-# Apply for different regions Lower/Interior/Upper or Unknown region
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Lower})
-    return @inbounds h*h*apply(op.closureStencils[Int(i)], v, Int(i))
-end
-
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Interior})
-    return @inbounds h*h*apply(op.innerStencil, v, Int(i))
-end
-
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Upper})
-    N = length(v)
-    return @inbounds h*h*Int(op.parity)*apply_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))
-end
-
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, index::Index{Unknown})
-    cSize = closureSize(op)
-    N = length(v)
-
-    i = Int(index)
-
-    if 0 < i <= cSize
-        return apply(op, h, v, Index{Lower}(i))
-    elseif cSize < i <= N-cSize
-        return apply(op, h, v, Index{Interior}(i))
-    elseif N-cSize < i <= N
-        return apply(op, h, v, Index{Upper}(i))
-    else
-        error("Bounds error") # TODO: Make this more standard
-    end
-end
-
-
-# Wrapper functions for using regular indecies without specifying regions
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int)
-    return apply(op, h, v, Index{Unknown}(i))
-end
-
-@enum Parity begin
-    odd = -1
-    even = 1
-end
-
-struct D2{T,N,M,K} <: ConstantStencilOperator
-    quadratureClosure::NTuple{M,T}
-    innerStencil::Stencil{T,N}
-    closureStencils::NTuple{M,Stencil{T,K}}
-    eClosure::Stencil{T,M}
-    dClosure::Stencil{T,M}
-    parity::Parity
-end
-
-function closureSize(D::D2)::Int
-    return length(D.quadratureClosure)
-end
-
-function readOperator(D2fn, Hfn)
-    d = readSectionedFile(D2fn)
-    h = readSectionedFile(Hfn)
-
-    # Create inner stencil
-    innerStencilWeights = stringToTuple(Float64, d["inner_stencil"][1])
-    width = length(innerStencilWeights)
-    r = (-div(width,2), div(width,2))
-
-    innerStencil = Stencil(r, innerStencilWeights)
-
-    # Create boundary stencils
-    boundarySize = length(d["boundary_stencils"])
-    closureStencils = Vector{typeof(innerStencil)}() # TBD: is the the right way to get the correct type?
-
-    for i ∈ 1:boundarySize
-        stencilWeights = stringToTuple(Float64, d["boundary_stencils"][i])
-        width = length(stencilWeights)
-        r = (1-i,width-i)
-        closureStencils = (closureStencils..., Stencil(r, stencilWeights))
-    end
-
-    quadratureClosure = pad_tuple(stringToTuple(Float64, h["closure"][1]), boundarySize)
-    eClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["e"][1]), boundarySize))
-    dClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["d1"][1]), boundarySize))
-
-    d2 = D2(
-        quadratureClosure,
-        innerStencil,
-        closureStencils,
-        eClosure,
-        dClosure,
-        even
-    )
-
-    return d2
-end
-
-
-function apply_e_T(op::D2, v::AbstractVector, ::Type{Lower})
-    @boundscheck if length(v) < closureSize(op)
-        throw(BoundsError())
-    end
-    apply(op.eClosure,v,1)
-end
-
-function apply_e_T(op::D2, v::AbstractVector, ::Type{Upper})
-    @boundscheck if length(v) < closureSize(op)
-        throw(BoundsError())
-    end
-    apply(flip(op.eClosure),v,length(v))
-end
-
-
-function apply_e(op::D2, v::Number, N::Integer, i::Integer, ::Type{Lower})
-    @boundscheck if !(0<length(i) <= N)
-        throw(BoundsError())
-    end
-    op.eClosure[i-1]*v
-end
-
-function apply_e(op::D2, v::Number, N::Integer, i::Integer, ::Type{Upper})
-    @boundscheck if !(0<length(i) <= N)
-        throw(BoundsError())
-    end
-    op.eClosure[N-i]*v
-end
-
-function apply_d_T(op::D2, h_inv::Real, v::AbstractVector, ::Type{Lower})
-    @boundscheck if length(v) < closureSize(op)
-        throw(BoundsError())
-    end
-    h_inv*apply(op.dClosure,v,1)
-end
-
-function apply_d_T(op::D2, h_inv::Real, v::AbstractVector, ::Type{Upper})
-    @boundscheck if length(v) < closureSize(op)
-        throw(BoundsError())
-    end
-    -h_inv*apply(flip(op.dClosure),v,length(v))
-end
-
-function apply_d(op::D2, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Lower})
-    @boundscheck if !(0<length(i) <= N)
-        throw(BoundsError())
-    end
-    h_inv*op.dClosure[i-1]*v
-end
-
-function apply_d(op::D2, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Upper})
-    @boundscheck if !(0<length(i) <= N)
-        throw(BoundsError())
-    end
-    -h_inv*op.dClosure[N-i]*v
-end
-
-
-function readSectionedFile(filename)::Dict{String, Vector{String}}
-    f = open(filename)
-    sections = Dict{String, Vector{String}}()
-    currentKey = ""
-
-    for ln ∈ eachline(f)
-        if ln == "" || ln[1] == '#' # Skip comments and empty lines
-            continue
-        end
-
-        if isletter(ln[1]) # Found start of new section
-            if ~haskey(sections, ln)
-                sections[ln] =  Vector{String}()
-            end
-            currentKey = ln
-            continue
-        end
-
-        push!(sections[currentKey], ln)
-    end
-
-    return sections
-end
-
-function stringToTuple(T::DataType, s::String)
-    return Tuple(stringToVector(T,s))
-end
-
-function stringToVector(T::DataType, s::String)
-    return T.(eval.(Meta.parse.(split(s))))
-end
-
-
-function pad_tuple(t::NTuple{N, T}, n::Integer) where {N,T}
-    if N >= n
-        return t
-    else
-        return pad_tuple((t..., zero(T)), n)
-    end
-end
-
-sbp_operators_path() = (@__DIR__) * "/../operators/"
-export sbp_operators_path
+include("constantstenciloperator.jl")
+include("d2.jl")
+include("readoperator.jl")
 
 end # module

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/constantstenciloperator.jl	Thu Jun 27 14:37:21 2019 +0200
@@ -0,0 +1,39 @@
+export apply
+
+abstract type ConstantStencilOperator end
+
+# Apply for different regions Lower/Interior/Upper or Unknown region
+@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Lower})
+    return @inbounds h*h*apply(op.closureStencils[Int(i)], v, Int(i))
+end
+
+@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Interior})
+    return @inbounds h*h*apply(op.innerStencil, v, Int(i))
+end
+
+@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Upper})
+    N = length(v)
+    return @inbounds h*h*Int(op.parity)*apply_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))
+end
+
+@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, index::Index{Unknown})
+    cSize = closureSize(op)
+    N = length(v)
+
+    i = Int(index)
+
+    if 0 < i <= cSize
+        return apply(op, h, v, Index{Lower}(i))
+    elseif cSize < i <= N-cSize
+        return apply(op, h, v, Index{Interior}(i))
+    elseif N-cSize < i <= N
+        return apply(op, h, v, Index{Upper}(i))
+    else
+        error("Bounds error") # TODO: Make this more standard
+    end
+end
+
+# Wrapper functions for using regular indecies without specifying regions
+@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int)
+    return apply(op, h, v, Index{Unknown}(i))
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/d2.jl	Thu Jun 27 14:37:21 2019 +0200
@@ -0,0 +1,76 @@
+export D2, closureSize, readOperator, apply_e, apply_d, apply_e_T, apply_d_T
+
+@enum Parity begin
+    odd = -1
+    even = 1
+end
+
+struct D2{T,N,M,K} <: ConstantStencilOperator
+    quadratureClosure::NTuple{M,T}
+    innerStencil::Stencil{T,N}
+    closureStencils::NTuple{M,Stencil{T,K}}
+    eClosure::Stencil{T,M}
+    dClosure::Stencil{T,M}
+    parity::Parity
+end
+
+function closureSize(D::D2)::Int
+    return length(D.quadratureClosure)
+end
+
+function apply_e_T(op::D2, v::AbstractVector, ::Type{Lower})
+    @boundscheck if length(v) < closureSize(op)
+        throw(BoundsError())
+    end
+    apply(op.eClosure,v,1)
+end
+
+function apply_e_T(op::D2, v::AbstractVector, ::Type{Upper})
+    @boundscheck if length(v) < closureSize(op)
+        throw(BoundsError())
+    end
+    apply(flip(op.eClosure),v,length(v))
+end
+
+
+function apply_e(op::D2, v::Number, N::Integer, i::Integer, ::Type{Lower})
+    @boundscheck if !(0<length(i) <= N)
+        throw(BoundsError())
+    end
+    op.eClosure[i-1]*v
+end
+
+function apply_e(op::D2, v::Number, N::Integer, i::Integer, ::Type{Upper})
+    @boundscheck if !(0<length(i) <= N)
+        throw(BoundsError())
+    end
+    op.eClosure[N-i]*v
+end
+
+function apply_d_T(op::D2, h_inv::Real, v::AbstractVector, ::Type{Lower})
+    @boundscheck if length(v) < closureSize(op)
+        throw(BoundsError())
+    end
+    h_inv*apply(op.dClosure,v,1)
+end
+
+function apply_d_T(op::D2, h_inv::Real, v::AbstractVector, ::Type{Upper})
+    @boundscheck if length(v) < closureSize(op)
+        throw(BoundsError())
+    end
+    -h_inv*apply(flip(op.dClosure),v,length(v))
+end
+
+function apply_d(op::D2, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Lower})
+    @boundscheck if !(0<length(i) <= N)
+        throw(BoundsError())
+    end
+    h_inv*op.dClosure[i-1]*v
+end
+
+function apply_d(op::D2, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Upper})
+    @boundscheck if !(0<length(i) <= N)
+        throw(BoundsError())
+    end
+    -h_inv*op.dClosure[N-i]*v
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/readoperator.jl	Thu Jun 27 14:37:21 2019 +0200
@@ -0,0 +1,80 @@
+function readOperator(D2fn, Hfn)
+    d = readSectionedFile(D2fn)
+    h = readSectionedFile(Hfn)
+
+    # Create inner stencil
+    innerStencilWeights = stringToTuple(Float64, d["inner_stencil"][1])
+    width = length(innerStencilWeights)
+    r = (-div(width,2), div(width,2))
+
+    innerStencil = Stencil(r, innerStencilWeights)
+
+    # Create boundary stencils
+    boundarySize = length(d["boundary_stencils"])
+    closureStencils = Vector{typeof(innerStencil)}() # TBD: is the the right way to get the correct type?
+
+    for i ∈ 1:boundarySize
+        stencilWeights = stringToTuple(Float64, d["boundary_stencils"][i])
+        width = length(stencilWeights)
+        r = (1-i,width-i)
+        closureStencils = (closureStencils..., Stencil(r, stencilWeights))
+    end
+
+    quadratureClosure = pad_tuple(stringToTuple(Float64, h["closure"][1]), boundarySize)
+    eClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["e"][1]), boundarySize))
+    dClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["d1"][1]), boundarySize))
+
+    d2 = D2(
+        quadratureClosure,
+        innerStencil,
+        closureStencils,
+        eClosure,
+        dClosure,
+        even
+    )
+
+    return d2
+end
+
+function readSectionedFile(filename)::Dict{String, Vector{String}}
+    f = open(filename)
+    sections = Dict{String, Vector{String}}()
+    currentKey = ""
+
+    for ln ∈ eachline(f)
+        if ln == "" || ln[1] == '#' # Skip comments and empty lines
+            continue
+        end
+
+        if isletter(ln[1]) # Found start of new section
+            if ~haskey(sections, ln)
+                sections[ln] =  Vector{String}()
+            end
+            currentKey = ln
+            continue
+        end
+
+        push!(sections[currentKey], ln)
+    end
+
+    return sections
+end
+
+function stringToTuple(T::DataType, s::String)
+    return Tuple(stringToVector(T,s))
+end
+
+function stringToVector(T::DataType, s::String)
+    return T.(eval.(Meta.parse.(split(s))))
+end
+
+function pad_tuple(t::NTuple{N, T}, n::Integer) where {N,T}
+    if N >= n
+        return t
+    else
+        return pad_tuple((t..., zero(T)), n)
+    end
+end
+
+sbp_operators_path() = (@__DIR__) * "/../operators/"
+export sbp_operators_path