Mercurial > repos > public > sbplib_julia

changeset 291:0f94dc29c4bf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in branch boundary_conditions
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 22 Jun 2020 21:43:05 +0200
parents fbabfd4e8f20 (current diff) 22d7bd87c20e (diff)
children 3747e5636eef 9c12d9eb38fd
files
diffstat 25 files changed, 1160 insertions(+), 496 deletions(-) [+]
line wrap: on
line diff
--- a/DiffOps/Manifest.toml	Wed Jun 26 15:07:47 2019 +0200
+++ b/DiffOps/Manifest.toml	Mon Jun 22 21:43:05 2020 +0200
@@ -17,6 +17,11 @@
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
+[[LazyTensors]]
+path = "../LazyTensors"
+uuid = "62fbed2c-918d-11e9-279b-eb3a325b37d3"
+version = "0.1.0"
+
 [[Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
--- a/DiffOps/Project.toml	Wed Jun 26 15:07:47 2019 +0200
+++ b/DiffOps/Project.toml	Mon Jun 22 21:43:05 2020 +0200
@@ -5,6 +5,7 @@
 
 [deps]
 Grids = "960fdf28-97ed-11e9-2a74-bd90bf2fab5a"
+LazyTensors = "62fbed2c-918d-11e9-279b-eb3a325b37d3"
 RegionIndices = "5d527584-97f1-11e9-084c-4540c7ecf219"
 SbpOperators = "204357d8-97fd-11e9-05e7-010897a14cd0"
 TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
--- a/DiffOps/src/DiffOps.jl	Wed Jun 26 15:07:47 2019 +0200
+++ b/DiffOps/src/DiffOps.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -3,6 +3,7 @@
 using RegionIndices
 using SbpOperators
 using Grids
+using LazyTensors
 
 """
     DiffOp
@@ -46,7 +47,7 @@
 
 # Maybe this should be split according to b3fbef345810 after all?! Seems like it makes performance more predictable
 function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
-    for I ∈ regionindices(D.grid.size, closureSize(D.op), (r1,r2))
+    for I ∈ regionindices(D.grid.size, closuresize(D.op), (r1,r2))
         @inbounds indextuple = (Index{r1}(I[1]), Index{r2}(I[2]))
         @inbounds u[I] = apply(D, v, indextuple)
     end
@@ -69,7 +70,7 @@
 
 using TiledIteration
 function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
-    ri = regionindices(D.grid.size, closureSize(D.op), (r1,r2))
+    ri = regionindices(D.grid.size, closuresize(D.op), (r1,r2))
     # TODO: Pass Tilesize to function
     for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2))
         for j ∈ tileaxs[2], i ∈ tileaxs[1]
@@ -97,7 +98,5 @@
 
 
 include("laplace.jl")
-export Laplace
-
 
 end # module
--- a/DiffOps/src/laplace.jl	Wed Jun 26 15:07:47 2019 +0200
+++ b/DiffOps/src/laplace.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -1,111 +1,205 @@
-struct NormalDerivative{N,M,K}
-	op::D2{Float64,N,M,K}
-	grid::EquidistantGrid
-	bId::CartesianBoundary
-end
-
-function apply_transpose(d::NormalDerivative, v::AbstractArray, I::Integer)
-	u = selectdim(v,3-dim(d.bId),I)
-	return apply_d(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId))
-end
-
-# Not correct abstraction level
-# TODO: Not type stable D:<
-function apply(d::NormalDerivative, v::AbstractArray, I::Tuple{Integer,Integer})
-	i = I[dim(d.bId)]
-	j = I[3-dim(d.bId)]
-	N_i = d.grid.size[dim(d.bId)]
-
-	r = getregion(i, closureSize(d.op), N_i)
-
-	if r != region(d.bId)
-		return 0
-	end
-
-	if r == Lower
-		# Note, closures are indexed by offset. Fix this D:<
-		return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[i-1]*v[j]
-	elseif r == Upper
-		return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[N_i-j]*v[j]
-	end
-end
+#TODO: move to sbpoperators.jl
+"""
+    Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
 
-struct BoundaryValue{N,M,K}
-	op::D2{Float64,N,M,K}
-	grid::EquidistantGrid
-	bId::CartesianBoundary
+Implements the Laplace operator `L` in Dim dimensions as a tensor operator
+The multi-dimensional tensor operator simply consists of a tuple of the 1D
+Laplace tensor operator as defined by ConstantLaplaceOperator.
+"""
+struct Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
+    tensorOps::NTuple(Dim,ConstantLaplaceOperator{T,N,M,K})
+    #TODO: Write a good constructor
 end
-
-function apply(e::BoundaryValue, v::AbstractArray, I::Tuple{Integer,Integer})
-	i = I[dim(e.bId)]
-	j = I[3-dim(e.bId)]
-	N_i = e.grid.size[dim(e.bId)]
-
-	r = getregion(i, closureSize(e.op), N_i)
-
-	if r != region(e.bId)
-		return 0
-	end
+export Laplace
 
-	if r == Lower
-		# Note, closures are indexed by offset. Fix this D:<
-		return e.op.eClosure[i-1]*v[j]
-	elseif r == Upper
-		return e.op.eClosure[N_i-j]*v[j]
-	end
-end
+LazyTensors.domain_size(H::Laplace{Dim}, range_size::NTuple{Dim,Integer}) = range_size
 
-function apply_transpose(e::BoundaryValue, v::AbstractArray, I::Integer)
-	u = selectdim(v,3-dim(e.bId),I)
-	return apply_e(e.op, u, region(e.bId))
-end
-
-struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim}
-    grid::EquidistantGrid{Dim,T}
-    a::T
-    op::D2{Float64,N,M,K}
-    # e::BoundaryValue
-    # d::NormalDerivative
-end
-
-function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim
+function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim}
     error("not implemented")
 end
 
 # u = L*v
-function apply(L::Laplace{1}, v::AbstractVector, i::Int)
-    uᵢ = L.a * SbpOperators.apply(L.op, L.grid.spacing[1], v, i)
+function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
+    return apply(L.tensorOps[1],v,I)
+end
+
+
+@inline function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
+    # 2nd x-derivative
+    @inbounds vx = view(v, :, Int(I[2]))
+    @inbounds uᵢ = apply(L.tensorOps[1], vx , (I[1],)) #Tuple conversion here is ugly. How to do it? Should we use indexing here?
+
+    # 2nd y-derivative
+    @inbounds vy = view(v, Int(I[1]), :)
+    @inbounds uᵢ += apply(L.tensorOps[2], vy , (I[2],)) #Tuple conversion here is ugly. How to do it?
+
     return uᵢ
 end
 
-@inline function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
-    # 2nd x-derivative
-    @inbounds vx = view(v, :, Int(I[2]))
-    @inbounds uᵢ = L.a*SbpOperators.apply(L.op, L.grid.inverse_spacing[1], vx , I[1])
-    # 2nd y-derivative
-    @inbounds vy = view(v, Int(I[1]), :)
-    @inbounds uᵢ += L.a*SbpOperators.apply(L.op, L.grid.inverse_spacing[2], vy, I[2])
-    # NOTE: the package qualifier 'SbpOperators' can problably be removed once all "applying" objects use LazyTensors
-    return uᵢ
+quadrature(L::Laplace) = Quadrature(L.op, L.grid)
+inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid)
+boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId)
+normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId)
+boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId)
+export quadrature
+
+# At the moment the grid property is used all over. It could possibly be removed if we implement all the 1D operators as TensorMappings
+"""
+    Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
+
+Implements the quadrature operator `H` of Dim dimension as a TensorMapping
+"""
+struct Quadrature{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
+    op::D2{T,N,M,K}
+    grid::EquidistantGrid{Dim,T}
+end
+export Quadrature
+
+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)
+    # Quadrature in x direction
+    @inbounds q = apply_quadrature(H.op, spacing(H.grid)[1], v[I] , I[1], N[1])
+    # Quadrature in y-direction
+    @inbounds q = apply_quadrature(H.op, spacing(H.grid)[2], q, I[2], N[2])
+    return q
+end
+
+LazyTensors.apply_transpose(H::Quadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T = LazyTensors.apply(H,v,I)
+
+
+"""
+    InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
+
+Implements the inverse quadrature operator `inv(H)` of Dim dimension as a TensorMapping
+"""
+struct InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
+    op::D2{T,N,M,K}
+    grid::EquidistantGrid{Dim,T}
+end
+export InverseQuadrature
+
+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)
+    # Inverse quadrature in x direction
+    @inbounds q_inv = apply_inverse_quadrature(H_inv.op, inverse_spacing(H_inv.grid)[1], v[I] , I[1], N[1])
+    # Inverse quadrature in y-direction
+    @inbounds q_inv = apply_inverse_quadrature(H_inv.op, inverse_spacing(H_inv.grid)[2], q_inv, I[2], N[2])
+    return q_inv
 end
 
-# Slow but maybe convenient?
-function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2})
-    I = Index{Unknown}.(Tuple(i))
-    apply(L, v, I)
+LazyTensors.apply_transpose(H_inv::InverseQuadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T = LazyTensors.apply(H_inv,v,I)
+
+"""
+    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{2}
+    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?
+function LazyTensors.range_size(e::BoundaryValue{T}, domain_size::NTuple{1,Integer}) where T
+    if dim(e.bId) == 1
+        return (UnknownDim, domain_size[1])
+    elseif dim(e.bId) == 2
+        return (domain_size[1], UnknownDim)
+    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
+    i = I[dim(e.bId)]
+    j = I[3-dim(e.bId)]
+    N_i = size(e.grid)[dim(e.bId)]
+    return apply_boundary_value(e.op, v[j], i, N_i, region(e.bId))
+end
+
+function LazyTensors.apply_transpose(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T
+    u = selectdim(v,3-dim(e.bId),Int(I[1]))
+    return apply_boundary_value_transpose(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{2}
+    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?
+function LazyTensors.range_size(e::NormalDerivative, domain_size::NTuple{1,Integer})
+    if dim(e.bId) == 1
+        return (UnknownDim, domain_size[1])
+    elseif dim(e.bId) == 2
+        return (domain_size[1], UnknownDim)
+    end
+end
+LazyTensors.domain_size(e::NormalDerivative, range_size::NTuple{2,Integer}) = (range_size[3-dim(e.bId)],)
+
+# TODO: Not type stable D:<
+# TODO: Make this independent of dimension
+function LazyTensors.apply(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T
+    i = I[dim(d.bId)]
+    j = I[3-dim(d.bId)]
+    N_i = size(d.grid)[dim(d.bId)]
+    h_inv = inverse_spacing(d.grid)[dim(d.bId)]
+    return apply_normal_derivative(d.op, h_inv, v[j], i, N_i, region(d.bId))
+end
+
+function LazyTensors.apply_transpose(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T
+    u = selectdim(v,3-dim(d.bId),Int(I[1]))
+    return apply_normal_derivative_transpose(d.op, inverse_spacing(d.grid)[dim(d.bId)], u, region(d.bId))
+end
+
+"""
+    BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1}
+
+Implements the boundary operator `q` as a TensorOperator
+"""
+struct BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1}
+    op::D2{T,N,M,K}
+    grid::EquidistantGrid{2}
+    bId::CartesianBoundary
+end
+export BoundaryQuadrature
+
+# TODO: Make this independent of dimension
+function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T
+    h = spacing(q.grid)[3-dim(q.bId)]
+    N = size(v)
+    return apply_quadrature(q.op, h, v[I[1]], I[1], N[1])
+end
+
+LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T = LazyTensors.apply(q,v,I)
+
+
 
 
 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}
-    e = BoundaryValue(L.op, L.grid, Bid())
-    d = NormalDerivative(L.op, L.grid, Bid())
-    Hᵧ = BoundaryQuadrature(L.op, L.grid, Bid())
-    # TODO: Implement BoundaryQuadrature method
-
-    return -L.Hi*e*Hᵧ*(d'*v - g)
-    # Need to handle d'*v - g so that it is an AbstractArray that TensorMappings can act on
+    e = boundary_value(L, Bid())
+    d = normal_derivative(L, Bid())
+    Hᵧ = boundary_quadrature(L, Bid())
+    H⁻¹ = inverse_quadrature(L)
+    return (-H⁻¹*e*Hᵧ*(d'*v - g))[I]
 end
 
 struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition
@@ -113,17 +207,16 @@
 end
 
 function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid}
-    e = BoundaryValue(L.op, L.grid, Bid())
-    d = NormalDerivative(L.op, L.grid, Bid())
-    Hᵧ = BoundaryQuadrature(L.op, L.grid, Bid())
-    # TODO: Implement BoundaryQuadrature method
-
-    return -L.Hi*(tau/h*e + d)*Hᵧ*(e'*v - g)
+    e = boundary_value(L, Bid())
+    d = normal_derivative(L, Bid())
+    Hᵧ = boundary_quadrature(L, Bid())
+    H⁻¹ = inverse_quadrature(L)
+    return (-H⁻¹*(tau/h*e + d)*Hᵧ*(e'*v - g))[I]
     # Need to handle scalar multiplication and addition of TensorMapping
 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/DiffOps/test/runtests.jl	Wed Jun 26 15:07:47 2019 +0200
+++ b/DiffOps/test/runtests.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -1,4 +1,269 @@
+using Test
 using DiffOps
-using Test
+using Grids
+using SbpOperators
+using RegionIndices
+using LazyTensors
+
+@testset "Laplace2D" begin
+    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    Lx = 3.5
+    Ly = 7.2
+    g = EquidistantGrid((42,41), (0.0, 0.0), (Lx,Ly))
+    L = Laplace(g, 1., op)
+    H = quadrature(L)
+
+    f0(x::Float64,y::Float64) = 2.
+    f1(x::Float64,y::Float64) = x+y
+    f2(x::Float64,y::Float64) = 1/2*x^2 + 1/2*y^2
+    f3(x::Float64,y::Float64) = 1/6*x^3 + 1/6*y^3
+    f4(x::Float64,y::Float64) = 1/24*x^4 + 1/24*y^4
+    f5(x::Float64,y::Float64) = sin(x) + cos(y)
+    f5ₓₓ(x::Float64,y::Float64) = -f5(x,y)
+
+    v0 = evalOn(g,f0)
+    v1 = evalOn(g,f1)
+    v2 = evalOn(g,f2)
+    v3 = evalOn(g,f3)
+    v4 = evalOn(g,f4)
+    v5 = evalOn(g,f5)
+    v5ₓₓ = evalOn(g,f5ₓₓ)
+
+    @test L isa TensorOperator{T,2} where T
+    @test L' isa TensorMapping{T,2,2} where T
+
+    # TODO: Should perhaps set tolerance level for isapporx instead?
+    #       Are these tolerance levels resonable or should tests be constructed
+    #       differently?
+    equalitytol = 0.5*1e-10
+    accuracytol = 0.5*1e-3
+    # 4th order interior stencil, 2nd order boundary stencil,
+    # implies that L*v should be exact for v - monomial up to order 3.
+    # Exact differentiation is measured point-wise. For other grid functions
+    # the error is measured in the H-norm.
+    @test all(abs.(collect(L*v0)) .<= equalitytol)
+    @test all(abs.(collect(L*v1)) .<= equalitytol)
+    @test all(collect(L*v2) .≈ v0) # Seems to be more accurate
+    @test all(abs.((collect(L*v3) - v1)) .<= equalitytol)
+    e4 = collect(L*v4) - v2
+    e5 = collect(L*v5) - v5ₓₓ
+    @test sum(collect(H*e4.^2)) <= accuracytol
+    @test sum(collect(H*e5.^2)) <= accuracytol
+end
+
+@testset "Quadrature" begin
+    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    Lx = 2.3
+    Ly = 5.2
+    g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
+    H = Quadrature(op,g)
+    v = ones(Float64, size(g))
+
+    @test H isa TensorOperator{T,2} where T
+    @test H' isa TensorMapping{T,2,2} where T
+    @test sum(collect(H*v)) ≈ (Lx*Ly)
+    @test collect(H*v) == collect(H'*v)
+end
+
+@testset "InverseQuadrature" begin
+    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    Lx = 7.3
+    Ly = 8.2
+    g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
+    H = Quadrature(op,g)
+    Hinv = InverseQuadrature(op,g)
+    v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
+
+    @test Hinv isa TensorOperator{T,2} where T
+    @test Hinv' isa TensorMapping{T,2,2} where T
+    @test collect(Hinv*H*v)  ≈ v
+    @test collect(Hinv*v) == collect(Hinv'*v)
+end
+
+@testset "BoundaryValue" begin
+    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
+
+    e_w = BoundaryValue(op, g, CartesianBoundary{1,Lower}())
+    e_e = BoundaryValue(op, g, CartesianBoundary{1,Upper}())
+    e_s = BoundaryValue(op, g, CartesianBoundary{2,Lower}())
+    e_n = BoundaryValue(op, g, CartesianBoundary{2,Upper}())
+
+    v = zeros(Float64, 4, 5)
+    v[:,5] = [1, 2, 3,4]
+    v[:,4] = [1, 2, 3,4]
+    v[:,3] = [4, 5, 6, 7]
+    v[:,2] = [7, 8, 9, 10]
+    v[:,1] = [10, 11, 12, 13]
+
+    @test e_w  isa TensorMapping{T,2,1} where T
+    @test e_w' isa TensorMapping{T,1,2} where T
+
+    @test domain_size(e_w, (3,2)) == (2,)
+    @test domain_size(e_e, (3,2)) == (2,)
+    @test domain_size(e_s, (3,2)) == (3,)
+    @test domain_size(e_n, (3,2)) == (3,)
+
+    @test size(e_w'*v) == (5,)
+    @test size(e_e'*v) == (5,)
+    @test size(e_s'*v) == (4,)
+    @test size(e_n'*v) == (4,)
+
+    @test collect(e_w'*v) == [10,7,4,1.0,1]
+    @test collect(e_e'*v) == [13,10,7,4,4.0]
+    @test collect(e_s'*v) == [10,11,12,13.0]
+    @test collect(e_n'*v) == [1,2,3,4.0]
+
+    g_x = [1,2,3,4.0]
+    g_y = [5,4,3,2,1.0]
+
+    G_w = zeros(Float64, (4,5))
+    G_w[1,:] = g_y
+
+    G_e = zeros(Float64, (4,5))
+    G_e[4,:] = g_y
+
+    G_s = zeros(Float64, (4,5))
+    G_s[:,1] = g_x
+
+    G_n = zeros(Float64, (4,5))
+    G_n[:,5] = g_x
+
+    @test size(e_w*g_y) == (UnknownDim,5)
+    @test size(e_e*g_y) == (UnknownDim,5)
+    @test size(e_s*g_x) == (4,UnknownDim)
+    @test size(e_n*g_x) == (4,UnknownDim)
 
-@test_broken false
+    # 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
+    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0))
+
+    d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}())
+    d_e = NormalDerivative(op, g, CartesianBoundary{1,Upper}())
+    d_s = NormalDerivative(op, g, CartesianBoundary{2,Lower}())
+    d_n = NormalDerivative(op, g, CartesianBoundary{2,Upper}())
+
+
+    v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
+    v∂x = evalOn(g, (x,y)-> 2*x + y)
+    v∂y = evalOn(g, (x,y)-> 2*(y-1) + x)
+
+    @test d_w  isa TensorMapping{T,2,1} where T
+    @test d_w' isa TensorMapping{T,1,2} where T
+
+    @test domain_size(d_w, (3,2)) == (2,)
+    @test domain_size(d_e, (3,2)) == (2,)
+    @test domain_size(d_s, (3,2)) == (3,)
+    @test domain_size(d_n, (3,2)) == (3,)
+
+    @test size(d_w'*v) == (6,)
+    @test size(d_e'*v) == (6,)
+    @test size(d_s'*v) == (5,)
+    @test size(d_n'*v) == (5,)
+
+    @test collect(d_w'*v) ≈ v∂x[1,:]
+    @test collect(d_e'*v) ≈ v∂x[5,:]
+    @test collect(d_s'*v) ≈ v∂y[:,1]
+    @test collect(d_n'*v) ≈ v∂y[:,6]
+
+
+    d_x_l = zeros(Float64, 5)
+    d_x_u = zeros(Float64, 5)
+    for i ∈ eachindex(d_x_l)
+        d_x_l[i] = op.dClosure[i-1]
+        d_x_u[i] = -op.dClosure[length(d_x_u)-i]
+    end
+
+    d_y_l = zeros(Float64, 6)
+    d_y_u = zeros(Float64, 6)
+    for i ∈ eachindex(d_y_l)
+        d_y_l[i] = op.dClosure[i-1]
+        d_y_u[i] = -op.dClosure[length(d_y_u)-i]
+    end
+
+    function prod_matrix(x,y)
+        G = zeros(Float64, length(x), length(y))
+        for I ∈ CartesianIndices(G)
+            G[I] = x[I[1]]*y[I[2]]
+        end
+
+        return G
+    end
+
+    g_x = [1,2,3,4.0,5]
+    g_y = [5,4,3,2,1.0,11]
+
+    G_w = prod_matrix(d_x_l, g_y)
+    G_e = prod_matrix(d_x_u, g_y)
+    G_s = prod_matrix(g_x, d_y_l)
+    G_n = prod_matrix(g_x, d_y_u)
+
+
+    @test size(d_w*g_y) == (UnknownDim,6)
+    @test size(d_e*g_y) == (UnknownDim,6)
+    @test size(d_s*g_x) == (5,UnknownDim)
+    @test size(d_n*g_x) == (5,UnknownDim)
+
+    # 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
+    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
+
+    @test collect(H_w'*v_w) == collect(H_w'*v_w)
+    @test collect(H_e'*v_e) == collect(H_e'*v_e)
+    @test collect(H_s'*v_s) == collect(H_s'*v_s)
+    @test collect(H_n'*v_n) == collect(H_n'*v_n)
+end
--- a/Grids/src/EquidistantGrid.jl	Wed Jun 26 15:07:47 2019 +0200
+++ b/Grids/src/EquidistantGrid.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -10,13 +10,13 @@
     size::NTuple{Dim, Int} # First coordinate direction stored first
     limit_lower::NTuple{Dim, T}
     limit_upper::NTuple{Dim, T}
-    inverse_spacing::NTuple{Dim, T} # The reciprocal of the grid spacing
+    inverse_spacing::NTuple{Dim, T} # Reciprocal of grid spacing
 
     # General constructor
     function EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T
         @assert all(size.>0)
         @assert all(limit_upper.-limit_lower .!= 0)
-        inverse_spacing = (size.-1)./abs.(limit_upper.-limit_lower)
+        inverse_spacing = (size.-1)./ abs.(limit_upper.-limit_lower)
         return new{Dim,T}(size, limit_lower, limit_upper, inverse_spacing)
     end
 end
@@ -25,6 +25,8 @@
     CartesianIndices(grid.size)
 end
 
+Base.size(g::EquidistantGrid) = g.size
+
 # Returns the number of dimensions of an EquidistantGrid.
 #
 # @Input: grid - an EquidistantGrid
@@ -33,11 +35,20 @@
     return length(grid.size)
 end
 
-# Returns the spacing of the grid
+# Returns the reciprocal of the spacing of the grid
 #
+function inverse_spacing(grid::EquidistantGrid)
+    return grid.inverse_spacing
+end
+export inverse_spacing
+
+# Returns the reciprocal of the spacing of the grid
+#
+# TODO: Evaluate if divisions affect performance
 function spacing(grid::EquidistantGrid)
     return 1.0./grid.inverse_spacing
 end
+export spacing
 
 # Computes the points of an EquidistantGrid as an array of tuples with
 # the same dimension as the grid.
--- a/Grids/src/Grids.jl	Wed Jun 26 15:07:47 2019 +0200
+++ b/Grids/src/Grids.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -9,6 +9,8 @@
 dim(::CartesianBoundary{Dim, R}) where {Dim, R} = Dim
 region(::CartesianBoundary{Dim, R}) where {Dim, R} = R
 
+export dim, region
+
 include("AbstractGrid.jl")
 include("EquidistantGrid.jl")
 
--- a/LazyTensors/Manifest.toml	Wed Jun 26 15:07:47 2019 +0200
+++ b/LazyTensors/Manifest.toml	Mon Jun 22 21:43:05 2020 +0200
@@ -1,2 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
+[[RegionIndices]]
+path = "../RegionIndices"
+uuid = "5d527584-97f1-11e9-084c-4540c7ecf219"
+version = "0.1.0"
--- a/LazyTensors/Project.toml	Wed Jun 26 15:07:47 2019 +0200
+++ b/LazyTensors/Project.toml	Mon Jun 22 21:43:05 2020 +0200
@@ -3,6 +3,9 @@
 authors = ["Jonatan Werpers <jonatan.werpers@it.uu.se>"]
 version = "0.1.0"
 
+[deps]
+RegionIndices = "5d527584-97f1-11e9-084c-4540c7ecf219"
+
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
--- a/LazyTensors/src/LazyTensors.jl	Wed Jun 26 15:07:47 2019 +0200
+++ b/LazyTensors/src/LazyTensors.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -1,6 +1,7 @@
 module LazyTensors
-
+using RegionIndices
 include("tensor_mapping.jl")
-include("lazy_operations.jl")
+include("lazy_array.jl")
+include("lazy_tensor_operations.jl")
 
 end # module
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/LazyTensors/src/lazy_array.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -0,0 +1,163 @@
+"""
+    LazyArray{T,D} <: AbstractArray{T,D}
+
+Array which is calcualted lazily when indexing.
+
+A subtype of `LazyArray` will use lazy version of `+`, `-`, `*`, `/`.
+"""
+abstract type LazyArray{T,D} <: AbstractArray{T,D} end
+export LazyArray
+
+"""
+    LazyElementwiseOperation{T,D,Op,T1,T2} <: LazyArray{T,D}
+Struct allowing for lazy evaluation of elementwise operations on AbstractArrays.
+
+A LazyElementwiseOperation contains two datatypes T1, and T2, together with an operation,
+where at least one of T1 and T2 is an AbstractArray, and one may be a Real.
+The operations are carried out when the LazyElementwiseOperation is indexed.
+"""
+struct LazyElementwiseOperation{T,D,Op,T1,T2} <: LazyArray{T,D}
+    a::T1
+    b::T2
+
+    @inline function LazyElementwiseOperation{T,D,Op}(a::T1,b::T2) where {T,D,Op,T1<:AbstractArray{T,D},T2<:AbstractArray{T,D}}
+        @boundscheck if size(a) != size(b)
+            throw(DimensionMismatch("dimensions must match"))
+        end
+        return new{T,D,Op,T1,T2}(a,b)
+    end
+
+    @inline function LazyElementwiseOperation{T,D,Op}(a::T1,b::T2) where {T,D,Op,T1<:AbstractArray{T,D},T2<:Real}
+        return new{T,D,Op,T1,T2}(a,b)
+    end
+
+    @inline function LazyElementwiseOperation{T,D,Op}(a::T1,b::T2) where {T,D,Op,T1<:Real,T2<:AbstractArray{T,D}}
+        return new{T,D,Op,T1,T2}(a,b)
+    end
+end
+# TODO: Move Op to be the first parameter? Compare to Binary operations
+
+Base.size(v::LazyElementwiseOperation) = size(v.a)
+
+# TODO: Make sure boundschecking is done properly and that the lenght of the vectors are equal
+# NOTE: Boundschecking in getindex functions now assumes that the size of the
+# vectors in the LazyElementwiseOperation are the same size. If we remove the
+# size assertion in the constructor we might have to handle
+# boundschecking differently.
+Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:+,T1,T2}, I::Vararg{Int,D}) where {T,D,T1<:AbstractArray{T,D},T2<:AbstractArray{T,D}}
+    @boundscheck if !checkbounds(Bool,leo.a,I...)
+        throw(BoundsError([leo],I...))
+    end
+    return leo.a[I...] + leo.b[I...]
+end
+
+Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:-,T1,T2}, I::Vararg{Int,D}) where {T,D,T1<:AbstractArray{T,D},T2<:AbstractArray{T,D}}
+    @boundscheck if !checkbounds(Bool,leo.a,I...)
+        throw(BoundsError([leo],I...))
+    end
+    return leo.a[I...] - leo.b[I...]
+end
+
+Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:*,T1,T2}, I::Vararg{Int,D}) where {T,D,T1<:AbstractArray{T,D},T2<:AbstractArray{T,D}}
+    @boundscheck if !checkbounds(Bool,leo.a,I...)
+        throw(BoundsError([leo],I...))
+    end
+    return leo.a[I...] * leo.b[I...]
+end
+
+Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:/,T1,T2}, I::Vararg{Int,D}) where {T,D,T1<:AbstractArray{T,D},T2<:AbstractArray{T,D}}
+    @boundscheck if !checkbounds(Bool,leo.a,I...)
+        throw(BoundsError([leo],I...))
+    end
+    return leo.a[I...] / leo.b[I...]
+end
+
+Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:+,T1,T2}, I::Vararg{Int,D}) where {T,D,T1<:AbstractArray{T,D},T2<:Real}
+    @boundscheck if !checkbounds(Bool,leo.a,I...)
+        throw(BoundsError([leo],I...))
+    end
+    return leo.a[I...] + leo.b
+end
+
+Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:-,T1,T2}, I::Vararg{Int,D}) where {T,D,T1<:AbstractArray{T,D},T2<:Real}
+    @boundscheck if !checkbounds(Bool,leo.a,I...)
+        throw(BoundsError([leo],I...))
+    end
+    return leo.a[I...] - leo.b
+end
+
+Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:*,T1,T2}, I::Vararg{Int,D}) where {T,D,T1<:AbstractArray{T,D},T2<:Real}
+    @boundscheck if !checkbounds(Bool,leo.a,I...)
+        throw(BoundsError([leo],I...))
+    end
+    return leo.a[I...] * leo.b
+end
+
+Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:/,T1,T2}, I::Vararg{Int,D}) where {T,D,T1<:AbstractArray{T,D},T2<:Real}
+    @boundscheck if !checkbounds(Bool,leo.a,I...)
+        throw(BoundsError([leo],I...))
+    end
+    return leo.a[I...] / leo.b
+end
+
+Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:+,T1,T2}, I::Vararg{Int,D}) where {T,D,T1<:Real,T2<:AbstractArray{T,D}}
+    @boundscheck if !checkbounds(Bool,leo.b,I...)
+        throw(BoundsError([leo],I...))
+    end
+    return leo.a + leo.b[I...]
+end
+
+Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:-,T1,T2}, I::Vararg{Int,D}) where {T,D,T1<:Real,T2<:AbstractArray{T,D}}
+    @boundscheck if !checkbounds(Bool,leo.b,I...)
+        throw(BoundsError([leo],I...))
+    end
+    return leo.a - leo.b[I...]
+end
+
+Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:*,T1,T2}, I::Vararg{Int,D}) where {T,D,T1<:Real,T2<:AbstractArray{T,D}}
+    @boundscheck if !checkbounds(Bool,leo.b,I...)
+        throw(BoundsError([leo],I...))
+    end
+    return leo.a * leo.b[I...]
+end
+
+Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:/,T1,T2}, I::Vararg{Int,D}) where {T,D,T1<:Real,T2<:AbstractArray{T,D}}
+    @boundscheck if !checkbounds(Bool,leo.b,I...)
+        throw(BoundsError([leo],I...))
+    end
+    return leo.a / leo.b[I...]
+end
+
+# Define lazy operations for AbstractArrays. Operations constructs a LazyElementwiseOperation which
+# can later be indexed into. Lazy operations are denoted by the usual operator followed by a tilde
+Base.@propagate_inbounds +̃(a::AbstractArray{T,D}, b::AbstractArray{T,D}) where {T,D} = LazyElementwiseOperation{T,D,:+}(a,b)
+Base.@propagate_inbounds -̃(a::AbstractArray{T,D}, b::AbstractArray{T,D}) where {T,D} = LazyElementwiseOperation{T,D,:-}(a,b)
+Base.@propagate_inbounds *̃(a::AbstractArray{T,D}, b::AbstractArray{T,D}) where {T,D} = LazyElementwiseOperation{T,D,:*}(a,b)
+Base.@propagate_inbounds /̃(a::AbstractArray{T,D}, b::AbstractArray{T,D}) where {T,D} = LazyElementwiseOperation{T,D,:/}(a,b)
+
+Base.@propagate_inbounds +̃(a::AbstractArray{T,D}, b::Real) where {T,D} = LazyElementwiseOperation{T,D,:+}(a,b)
+Base.@propagate_inbounds -̃(a::AbstractArray{T,D}, b::Real) where {T,D} = LazyElementwiseOperation{T,D,:-}(a,b)
+Base.@propagate_inbounds *̃(a::AbstractArray{T,D}, b::Real) where {T,D} = LazyElementwiseOperation{T,D,:*}(a,b)
+Base.@propagate_inbounds /̃(a::AbstractArray{T,D}, b::Real) where {T,D} = LazyElementwiseOperation{T,D,:/}(a,b)
+
+Base.@propagate_inbounds +̃(a::Real, b::AbstractArray{T,D}) where {T,D} = LazyElementwiseOperation{T,D,:+}(a,b)
+Base.@propagate_inbounds -̃(a::Real, b::AbstractArray{T,D}) where {T,D} = LazyElementwiseOperation{T,D,:-}(a,b)
+Base.@propagate_inbounds *̃(a::Real, b::AbstractArray{T,D}) where {T,D} = LazyElementwiseOperation{T,D,:*}(a,b)
+Base.@propagate_inbounds /̃(a::Real, b::AbstractArray{T,D}) where {T,D} = LazyElementwiseOperation{T,D,:/}(a,b)
+
+
+
+# NOTE: Är det knas att vi har till exempel * istället för .* ??
+# Oklart om det ens går att lösa..
+Base.@propagate_inbounds Base.:+(a::LazyArray{T,D}, b::LazyArray{T,D}) where {T,D} = a +̃ b
+Base.@propagate_inbounds Base.:+(a::LazyArray{T,D}, b::AbstractArray{T,D}) where {T,D} = a +̃ b
+Base.@propagate_inbounds Base.:+(a::AbstractArray{T,D}, b::LazyArray{T,D}) where {T,D} = a +̃ b
+
+Base.@propagate_inbounds Base.:-(a::LazyArray{T,D}, b::LazyArray{T,D}) where {T,D} = a -̃ b
+Base.@propagate_inbounds Base.:-(a::LazyArray{T,D}, b::AbstractArray{T,D}) where {T,D} = a -̃ b
+Base.@propagate_inbounds Base.:-(a::AbstractArray{T,D}, b::LazyArray{T,D}) where {T,D} = a -̃ b
+
+# Element wise operation for `*` and `\` are not overloaded due to conflicts with the behavior
+# of regular `*` and `/` for AbstractArrays. Use tilde versions instead.
+
+export +̃, -̃, *̃, /̃
--- a/LazyTensors/src/lazy_operations.jl	Wed Jun 26 15:07:47 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,193 +0,0 @@
-"""
-    LazyArray{T,D} <: AbstractArray{T,D}
-
-Array which is calcualted lazily when indexing.
-
-A subtype of `LazyArray` will use lazy version of `+`, `-`, `*`, `/`.
-"""
-abstract type LazyArray{T,D} <: AbstractArray{T,D} end
-export LazyArray
-
-
-
-"""
-    LazyTensorMappingApplication{T,R,D} <: LazyArray{T,R}
-
-Struct for lazy application of a TensorMapping. Created using `*`.
-
-Allows the result of a `TensorMapping` applied to a vector to be treated as an `AbstractArray`.
-With a mapping `m` and a vector `v` the LazyTensorMappingApplication object can be created by `m*v`.
-The actual result will be calcualted when indexing into `m*v`.
-"""
-struct LazyTensorMappingApplication{T,R,D} <: LazyArray{T,R}
-    t::TensorMapping{T,R,D}
-    o::AbstractArray{T,D}
-end
-export LazyTensorMappingApplication
-
-Base.:*(tm::TensorMapping{T,R,D}, o::AbstractArray{T,D}) where {T,R,D} = LazyTensorMappingApplication(tm,o)
-
-Base.getindex(ta::LazyTensorMappingApplication{T,R,D}, I::Vararg) where {T,R,D} = apply(ta.t, ta.o, I...)
-Base.size(ta::LazyTensorMappingApplication{T,R,D}) where {T,R,D} = range_size(ta.t,size(ta.o))
-# TODO: What else is needed to implement the AbstractArray interface?
-
-# # We need the associativity to be a→b→c = a→(b→c), which is the case for '→'
-Base.:*(args::Union{TensorMapping{T}, AbstractArray{T}}...) where T = foldr(*,args)
-# # Should we overload some other infix binary operator?
-# →(tm::TensorMapping{T,R,D}, o::AbstractArray{T,D}) where {T,R,D} = LazyTensorMappingApplication(tm,o)
-# TODO: We need to be really careful about good error messages.
-# For example what happens if you try to multiply LazyTensorMappingApplication with a TensorMapping(wrong order)?
-
-
-
-"""
-    LazyElementwiseOperation{T,D,Op, T1<:AbstractArray{T,D}, T2 <: AbstractArray{T,D}} <: AbstractArray{T,D}
-
-Struct allowing for lazy evaluation of elementwise operations on AbstractArrays.
-
-A LazyElementwiseOperation contains two AbstractArrays of equal size,
-together with an operation. The operations are carried out when the
-LazyElementwiseOperation is indexed.
-"""
-struct LazyElementwiseOperation{T,D,Op, T1<:AbstractArray{T,D}, T2 <: AbstractArray{T,D}} <: LazyArray{T,D}
-    a::T1
-    b::T2
-
-    @inline function LazyElementwiseOperation{T,D,Op}(a::T1,b::T2) where {T,D,Op, T1<:AbstractArray{T,D}, T2<:AbstractArray{T,D}}
-        @boundscheck if size(a) != size(b)
-            throw(DimensionMismatch("dimensions must match"))
-        end
-        return new{T,D,Op,T1,T2}(a,b)
-    end
-end
-# TODO: Move Op to be the first parameter? Compare to Binary operations
-
-Base.size(v::LazyElementwiseOperation) = size(v.a)
-
-# TODO: Make sure boundschecking is done properly and that the lenght of the vectors are equal
-# NOTE: Boundschecking in getindex functions now assumes that the size of the
-# vectors in the LazyElementwiseOperation are the same size. If we remove the
-# size assertion in the constructor we might have to handle
-# boundschecking differently.
-Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:+}, I...) where {T,D}
-    @boundscheck if !checkbounds(Bool,leo.a,I...)
-        throw(BoundsError([leo],[I...]))
-    end
-    return leo.a[I...] + leo.b[I...]
-end
-Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:-}, I...) where {T,D}
-    @boundscheck if !checkbounds(Bool,leo.a,I...)
-        throw(BoundsError([leo],[I...]))
-    end
-    return leo.a[I...] - leo.b[I...]
-end
-Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:*}, I...) where {T,D}
-    @boundscheck if !checkbounds(Bool,leo.a,I...)
-        throw(BoundsError([leo],[I...]))
-    end
-    return leo.a[I...] * leo.b[I...]
-end
-Base.@propagate_inbounds @inline function Base.getindex(leo::LazyElementwiseOperation{T,D,:/}, I...) where {T,D}
-    @boundscheck if !checkbounds(Bool,leo.a,I...)
-        throw(BoundsError([leo],[I...]))
-    end
-    return leo.a[I...] / leo.b[I...]
-end
-
-# Define lazy operations for AbstractArrays. Operations constructs a LazyElementwiseOperation which
-# can later be indexed into. Lazy operations are denoted by the usual operator followed by a tilde
-Base.@propagate_inbounds +̃(a::AbstractArray{T,D}, b::AbstractArray{T,D}) where {T,D} = LazyElementwiseOperation{T,D,:+}(a,b)
-Base.@propagate_inbounds -̃(a::AbstractArray{T,D}, b::AbstractArray{T,D}) where {T,D} = LazyElementwiseOperation{T,D,:-}(a,b)
-Base.@propagate_inbounds *̃(a::AbstractArray{T,D}, b::AbstractArray{T,D}) where {T,D} = LazyElementwiseOperation{T,D,:*}(a,b)
-Base.@propagate_inbounds /̃(a::AbstractArray{T,D}, b::AbstractArray{T,D}) where {T,D} = LazyElementwiseOperation{T,D,:/}(a,b)
-
-# NOTE: Är det knas att vi har till exempel * istället för .* ??
-# Oklart om det ens går att lösa..
-Base.@propagate_inbounds Base.:+(a::LazyArray{T,D}, b::LazyArray{T,D}) where {T,D} = a +̃ b
-Base.@propagate_inbounds Base.:+(a::LazyArray{T,D}, b::AbstractArray{T,D}) where {T,D} = a +̃ b
-Base.@propagate_inbounds Base.:+(a::AbstractArray{T,D}, b::LazyArray{T,D}) where {T,D} = a +̃ b
-
-Base.@propagate_inbounds Base.:-(a::LazyArray{T,D}, b::LazyArray{T,D}) where {T,D} = a -̃ b
-Base.@propagate_inbounds Base.:-(a::LazyArray{T,D}, b::AbstractArray{T,D}) where {T,D} = a -̃ b
-Base.@propagate_inbounds Base.:-(a::AbstractArray{T,D}, b::LazyArray{T,D}) where {T,D} = a -̃ b
-
-# Element wise operation for `*` and `\` are not overloaded due to conflicts with the behavior
-# of regular `*` and `/` for AbstractArrays. Use tilde versions instead.
-
-export +̃, -̃, *̃, /̃
-
-
-
-"""
-    LazyTensorMappingTranspose{T,R,D} <: TensorMapping{T,D,R}
-
-Struct for lazy transpose of a TensorMapping.
-
-If a mapping implements the the `apply_transpose` method this allows working with
-the transpose of mapping `m` by using `m'`. `m'` will work as a regular TensorMapping lazily calling
-the appropriate methods of `m`.
-"""
-struct LazyTensorMappingTranspose{T,R,D} <: TensorMapping{T,D,R}
-    tm::TensorMapping{T,R,D}
-end
-export LazyTensorMappingTranspose
-
-# # TBD: Should this be implemented on a type by type basis or through a trait to provide earlier errors?
-Base.adjoint(t::TensorMapping) = LazyTensorMappingTranspose(t)
-Base.adjoint(t::LazyTensorMappingTranspose) = t.tm
-
-apply(tm::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,R}, I::Vararg) where {T,R,D} = apply_transpose(tm.tm, v, I...)
-apply_transpose(tm::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {T,R,D} = apply(tm.tm, v, I...)
-
-range_size(tmt::LazyTensorMappingTranspose{T,R,D}, d_size::NTuple{R,Integer}) where {T,R,D} = domain_size(tmt.tm, domain_size)
-domain_size(tmt::LazyTensorMappingTranspose{T,R,D}, r_size::NTuple{D,Integer}) where {T,R,D} = range_size(tmt.tm, range_size)
-
-
-
-
-struct LazyTensorMappingBinaryOperation{Op,T,R,D,T1<:TensorMapping{T,R,D},T2<:TensorMapping{T,R,D}} <: TensorMapping{T,D,R}
-    A::T1
-    B::T2
-
-    @inline function LazyTensorMappingBinaryOperation{Op,T,R,D}(A::T1,B::T2) where {Op,T,R,D, T1<:TensorMapping{T,R,D},T2<:TensorMapping{T,R,D}}
-        return new{Op,T,R,D,T1,T2}(A,B)
-    end
-end
-
-apply(mb::LazyTensorMappingBinaryOperation{:+,T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {T,R,D} = apply(mb.A, v, I...) + apply(mb.B,v,I...)
-apply(mb::LazyTensorMappingBinaryOperation{:-,T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {T,R,D} = apply(mb.A, v, I...) - apply(mb.B,v,I...)
-
-range_size(mp::LazyTensorMappingBinaryOperation{Op,T,R,D}, domain_size::NTuple{D,Integer}) where {Op,T,R,D} = range_size(mp.A, domain_size)
-domain_size(mp::LazyTensorMappingBinaryOperation{Op,T,R,D}, range_size::NTuple{R,Integer}) where {Op,T,R,D} = domain_size(mp.A, range_size)
-
-Base.:+(A::TensorMapping{T,R,D}, B::TensorMapping{T,R,D}) where {T,R,D} = LazyTensorMappingBinaryOperation{:+,T,R,D}(A,B)
-Base.:-(A::TensorMapping{T,R,D}, B::TensorMapping{T,R,D}) where {T,R,D} = LazyTensorMappingBinaryOperation{:-,T,R,D}(A,B)
-
-
-# TODO: Write tests and documentation for LazyTensorMappingComposition
-# struct LazyTensorMappingComposition{T,R,K,D} <: TensorMapping{T,R,D}
-#     t1::TensorMapping{T,R,K}
-#     t2::TensorMapping{T,K,D}
-# end
-
-# Base.:∘(s::TensorMapping{T,R,K}, t::TensorMapping{T,K,D}) where {T,R,K,D} = LazyTensorMappingComposition(s,t)
-
-# function range_size(tm::LazyTensorMappingComposition{T,R,K,D}, domain_size::NTuple{D,Integer}) where {T,R,K,D}
-#     range_size(tm.t1, domain_size(tm.t2, domain_size))
-# end
-
-# function domain_size(tm::LazyTensorMappingComposition{T,R,K,D}, range_size::NTuple{R,Integer}) where {T,R,K,D}
-#     domain_size(tm.t1, domain_size(tm.t2, range_size))
-# end
-
-# function apply(c::LazyTensorMappingComposition{T,R,K,D}, v::AbstractArray{T,D}, I::Vararg) where {T,R,K,D}
-#     apply(c.t1, LazyTensorMappingApplication(c.t2,v), I...)
-# end
-
-# function apply_transpose(c::LazyTensorMappingComposition{T,R,K,D}, v::AbstractArray{T,D}, I::Vararg) where {T,R,K,D}
-#     apply_transpose(c.t2, LazyTensorMappingApplication(c.t1',v), I...)
-# end
-
-# # Have i gone too crazy with the type parameters? Maybe they aren't all needed?
-
-# export →
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/LazyTensors/src/lazy_tensor_operations.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -0,0 +1,101 @@
+"""
+    LazyTensorMappingApplication{T,R,D} <: LazyArray{T,R}
+
+Struct for lazy application of a TensorMapping. Created using `*`.
+
+Allows the result of a `TensorMapping` applied to a vector to be treated as an `AbstractArray`.
+With a mapping `m` and a vector `v` the LazyTensorMappingApplication object can be created by `m*v`.
+The actual result will be calcualted when indexing into `m*v`.
+"""
+struct LazyTensorMappingApplication{T,R,D} <: LazyArray{T,R}
+    t::TensorMapping{T,R,D}
+    o::AbstractArray{T,D}
+end
+export LazyTensorMappingApplication
+
+Base.:*(tm::TensorMapping{T,R,D}, o::AbstractArray{T,D}) where {T,R,D} = LazyTensorMappingApplication(tm,o)
+Base.getindex(ta::LazyTensorMappingApplication{T,R,D}, I::Vararg{Index,R}) where {T,R,D} = apply(ta.t, ta.o, I)
+Base.getindex(ta::LazyTensorMappingApplication{T,R,D}, I::Vararg{Int,R}) where {T,R,D} = apply(ta.t, ta.o, Index{Unknown}.(I))
+Base.size(ta::LazyTensorMappingApplication{T,R,D}) where {T,R,D} = range_size(ta.t,size(ta.o))
+# TODO: What else is needed to implement the AbstractArray interface?
+
+# # We need the associativity to be a→b→c = a→(b→c), which is the case for '→'
+Base.:*(a::TensorMapping{T,R,D}, b::TensorMapping{T,D,K}, args::Union{TensorMapping{T}, AbstractArray{T}}...) where {T,R,D,K} = foldr(*,(a,b,args...))
+# # Should we overload some other infix binary opesrator?
+# →(tm::TensorMapping{T,R,D}, o::AbstractArray{T,D}) where {T,R,D} = LazyTensorMappingApplication(tm,o)
+# TODO: We need to be really careful about good error messages.
+# For example what happens if you try to multiply LazyTensorMappingApplication with a TensorMapping(wrong order)?
+
+"""
+    LazyTensorMappingTranspose{T,R,D} <: TensorMapping{T,D,R}
+
+Struct for lazy transpose of a TensorMapping.
+
+If a mapping implements the the `apply_transpose` method this allows working with
+the transpose of mapping `m` by using `m'`. `m'` will work as a regular TensorMapping lazily calling
+the appropriate methods of `m`.
+"""
+struct LazyTensorMappingTranspose{T,R,D} <: TensorMapping{T,D,R}
+    tm::TensorMapping{T,R,D}
+end
+export LazyTensorMappingTranspose
+
+# # TBD: Should this be implemented on a type by type basis or through a trait to provide earlier errors?
+Base.adjoint(tm::TensorMapping) = LazyTensorMappingTranspose(tm)
+Base.adjoint(tmt::LazyTensorMappingTranspose) = tmt.tm
+
+apply(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,R}, I::NTuple{D,Index}) where {T,R,D} = apply_transpose(tmt.tm, v, I)
+apply_transpose(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,D}, I::NTuple{R,Index}) where {T,R,D} = apply(tmt.tm, v, I)
+
+range_size(tmt::LazyTensorMappingTranspose{T,R,D}, d_size::NTuple{R,Integer}) where {T,R,D} = domain_size(tmt.tm, d_size)
+domain_size(tmt::LazyTensorMappingTranspose{T,R,D}, r_size::NTuple{D,Integer}) where {T,R,D} = range_size(tmt.tm, r_size)
+
+
+
+
+struct LazyTensorMappingBinaryOperation{Op,T,R,D,T1<:TensorMapping{T,R,D},T2<:TensorMapping{T,R,D}} <: TensorMapping{T,D,R}
+    tm1::T1
+    tm2::T2
+
+    @inline function LazyTensorMappingBinaryOperation{Op,T,R,D}(tm1::T1,tm2::T2) where {Op,T,R,D, T1<:TensorMapping{T,R,D},T2<:TensorMapping{T,R,D}}
+        return new{Op,T,R,D,T1,T2}(tm1,tm2)
+    end
+end
+
+apply(tmBinOp::LazyTensorMappingBinaryOperation{:+,T,R,D}, v::AbstractArray{T,D}, I::NTuple{R,Index}) where {T,R,D} = apply(tmBinOp.tm1, v, I) + apply(tmBinOp.tm2, v, I)
+apply(tmBinOp::LazyTensorMappingBinaryOperation{:-,T,R,D}, v::AbstractArray{T,D}, I::NTuple{R,Index}) where {T,R,D} = apply(tmBinOp.tm1, v, I) - apply(tmBinOp.tm2, v, I)
+
+range_size(tmBinOp::LazyTensorMappingBinaryOperation{Op,T,R,D}, domain_size::NTuple{D,Integer}) where {Op,T,R,D} = range_size(tmBinOp.tm1, domain_size)
+domain_size(tmBinOp::LazyTensorMappingBinaryOperation{Op,T,R,D}, range_size::NTuple{R,Integer}) where {Op,T,R,D} = domain_size(tmBinOp.tm2, range_size)
+
+Base.:+(tm1::TensorMapping{T,R,D}, tm2::TensorMapping{T,R,D}) where {T,R,D} = LazyTensorMappingBinaryOperation{:+,T,R,D}(tm1,tm2)
+Base.:-(tm1::TensorMapping{T,R,D}, tm2::TensorMapping{T,R,D}) where {T,R,D} = LazyTensorMappingBinaryOperation{:-,T,R,D}(tm1,tm2)
+
+
+# TODO: Write tests and documentation for LazyTensorMappingComposition
+# struct LazyTensorMappingComposition{T,R,K,D} <: TensorMapping{T,R,D}
+#     t1::TensorMapping{T,R,K}
+#     t2::TensorMapping{T,K,D}
+# end
+
+# Base.:∘(s::TensorMapping{T,R,K}, t::TensorMapping{T,K,D}) where {T,R,K,D} = LazyTensorMappingComposition(s,t)
+
+# function range_size(tm::LazyTensorMappingComposition{T,R,K,D}, domain_size::NTuple{D,Integer}) where {T,R,K,D}
+#     range_size(tm.t1, domain_size(tm.t2, domain_size))
+# end
+
+# function domain_size(tm::LazyTensorMappingComposition{T,R,K,D}, range_size::NTuple{R,Integer}) where {T,R,K,D}
+#     domain_size(tm.t1, domain_size(tm.t2, range_size))
+# end
+
+# function apply(c::LazyTensorMappingComposition{T,R,K,D}, v::AbstractArray{T,D}, I::NTuple{R,Int}) where {T,R,K,D}
+#     apply(c.t1, LazyTensorMappingApplication(c.t2,v), I...)
+# end
+
+# function apply_transpose(c::LazyTensorMappingComposition{T,R,K,D}, v::AbstractArray{T,D}, I::NTuple{D,Int}) where {T,R,K,D}
+#     apply_transpose(c.t2, LazyTensorMappingApplication(c.t1',v), I...)
+# end
+
+# # Have i gone too crazy with the type parameters? Maybe they aren't all needed?
+
+# export →
--- a/LazyTensors/src/tensor_mapping.jl	Wed Jun 26 15:07:47 2019 +0200
+++ b/LazyTensors/src/tensor_mapping.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -62,7 +62,12 @@
 """
 function domain_size end
 
-export range_size, domain_size
+"""
+    Dummy type for representing dimensions of tensormappings when domain_size is unknown
+"""
+struct UnknownDim end
+export range_size, domain_size, TensorMappingDim, UnknownDim
+
 # TODO: Think about boundschecking!
 
 
--- a/LazyTensors/test/runtests.jl	Wed Jun 26 15:07:47 2019 +0200
+++ b/LazyTensors/test/runtests.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -1,12 +1,13 @@
 using Test
 using LazyTensors
+using RegionIndices
 
 @testset "Generic Mapping methods" begin
     struct DummyMapping{T,R,D} <: TensorMapping{T,R,D} end
-    LazyTensors.apply(m::DummyMapping{T,R,D}, v, i) where {T,R,D} = :apply
+    LazyTensors.apply(m::DummyMapping{T,R,D}, v, i::NTuple{R,Index{<:Region}}) where {T,R,D} = :apply
     @test range_dim(DummyMapping{Int,2,3}()) == 2
     @test domain_dim(DummyMapping{Int,2,3}()) == 3
-    @test apply(DummyMapping{Int,2,3}(), zeros(Int, (0,0,0)),0) == :apply
+    @test apply(DummyMapping{Int,2,3}(), zeros(Int, (0,0,0)),(Index{Unknown}(0),Index{Unknown}(0))) == :apply
 end
 
 @testset "Generic Operator methods" begin
@@ -18,17 +19,19 @@
 @testset "Mapping transpose" begin
     struct DummyMapping{T,R,D} <: TensorMapping{T,R,D} end
 
-    LazyTensors.apply(m::DummyMapping{T,R,D}, v, i) where {T,R,D} = :apply
-    LazyTensors.apply_transpose(m::DummyMapping{T,R,D}, v, i) where {T,R,D} = :apply_transpose
+    LazyTensors.apply(m::DummyMapping{T,R,D}, v, I::NTuple{R,Index{<:Region}}) where {T,R,D} = :apply
+    LazyTensors.apply_transpose(m::DummyMapping{T,R,D}, v, I::NTuple{D,Index{<:Region}}) where {T,R,D} = :apply_transpose
 
-    LazyTensors.range_size(m::DummyMapping{T,R,D}, domain_size) where {T,R,D} = :range_size
-    LazyTensors.domain_size(m::DummyMapping{T,R,D}, range_size) where {T,R,D} = :domain_size
+    LazyTensors.range_size(m::DummyMapping{T,R,D}, domain_size::NTuple{D,Integer}) where {T,R,D} = :range_size
+    LazyTensors.domain_size(m::DummyMapping{T,R,D}, range_size::NTuple{R,Integer}) where {T,R,D} = :domain_size
 
     m = DummyMapping{Float64,2,3}()
+    I = Index{Unknown}(0)
+    @test m' isa TensorMapping{Float64, 3,2}
     @test m'' == m
-    @test apply(m',zeros(Float64,(0,0)),0) == :apply_transpose
-    @test apply(m'',zeros(Float64,(0,0,0)),0) == :apply
-    @test apply_transpose(m', zeros(Float64,(0,0,0)),0) == :apply
+    @test apply(m',zeros(Float64,(0,0)), (I,I,I)) == :apply_transpose
+    @test apply(m'',zeros(Float64,(0,0,0)),(I,I)) == :apply
+    @test apply_transpose(m', zeros(Float64,(0,0,0)),(I,I)) == :apply
 
     @test range_size(m', (0,0)) == :domain_size
     @test domain_size(m', (0,0,0)) == :range_size
@@ -37,24 +40,53 @@
 @testset "TensorApplication" begin
     struct DummyMapping{T,R,D} <: TensorMapping{T,R,D} end
 
-    LazyTensors.apply(m::DummyMapping{T,R,D}, v, i) where {T,R,D} = (:apply,v,i)
-    LazyTensors.apply_transpose(m::DummyMapping{T,R,D}, v, i) where {T,R,D} = :apply_transpose
-
-    LazyTensors.range_size(m::DummyMapping{T,R,D}, domain_size) where {T,R,D} = 2 .* domain_size
-    LazyTensors.domain_size(m::DummyMapping{T,R,D}, range_size) where {T,R,D} = range_size.÷2
+    LazyTensors.apply(m::DummyMapping{T,R,D}, v, i::NTuple{R,Index{<:Region}}) where {T,R,D} = (:apply,v,i)
+    LazyTensors.range_size(m::DummyMapping{T,R,D}, domain_size::NTuple{D,Integer}) where {T,R,D} = 2 .* domain_size
+    LazyTensors.domain_size(m::DummyMapping{T,R,D}, range_size::NTuple{R,Integer}) where {T,R,D} = range_size.÷2
 
 
     m = DummyMapping{Int, 1, 1}()
     v = [0,1,2]
     @test m*v isa AbstractVector{Int}
     @test size(m*v) == 2 .*size(v)
-    @test (m*v)[0] == (:apply,v,0)
+    @test (m*v)[Index{Upper}(0)] == (:apply,v,(Index{Upper}(0),))
+    @test (m*v)[0] == (:apply,v,(Index{Unknown}(0),))
     @test m*m*v isa AbstractVector{Int}
-    @test (m*m*v)[1] == (:apply,m*v,1)
-    @test (m*m*v)[3] == (:apply,m*v,3)
-    @test (m*m*v)[6] == (:apply,m*v,6)
+    @test (m*m*v)[Index{Upper}(1)] == (:apply,m*v,(Index{Upper}(1),))
+    @test (m*m*v)[1] == (:apply,m*v,(Index{Unknown}(1),))
+    @test (m*m*v)[Index{Interior}(3)] == (:apply,m*v,(Index{Interior}(3),))
+    @test (m*m*v)[3] == (:apply,m*v,(Index{Unknown}(3),))
+    @test (m*m*v)[Index{Lower}(6)] == (:apply,m*v,(Index{Lower}(6),))
+    @test (m*m*v)[6] == (:apply,m*v,(Index{Unknown}(6),))
     @test_broken BoundsError == (m*m*v)[0]
     @test_broken BoundsError == (m*m*v)[7]
+
+    m = DummyMapping{Int, 2, 1}()
+    @test_throws MethodError m*ones(Int,2,2)
+    @test_throws MethodError m*m*v
+
+    m = DummyMapping{Float64, 2, 2}()
+    v = ones(3,3)
+    I = (Index{Lower}(1),Index{Interior}(2));
+    @test size(m*v) == 2 .*size(v)
+    @test (m*v)[I] == (:apply,v,I)
+
+    struct ScalingOperator{T,D} <: TensorOperator{T,D}
+        λ::T
+    end
+
+    LazyTensors.apply(m::ScalingOperator{T,D}, v, I::NTuple{D, Index}) where {T,D} = m.λ*v[I]
+
+    m = ScalingOperator{Int,1}(2)
+    v = [1,2,3]
+    @test m*v isa AbstractVector
+    @test m*v == [2,4,6]
+
+    m = ScalingOperator{Int,2}(2)
+    v = [[1 2];[3 4]]
+    @test m*v == [[2 4];[6 8]]
+    I = (Index{Upper}(2),Index{Lower}(1))
+    @test (m*v)[I] == 6
 end
 
 @testset "TensorMapping binary operations" begin
@@ -62,7 +94,7 @@
         λ::T
     end
 
-    LazyTensors.apply(m::ScalarMapping{T,R,D}, v, i) where {T,R,D} = m.λ*v[i]
+    LazyTensors.apply(m::ScalarMapping{T,R,D}, v, I::Tuple{Index{<:Region}}) where {T,R,D} = m.λ*v[I...]
     LazyTensors.range_size(m::ScalarMapping, domain_size) = domain_size
     LazyTensors.domain_size(m::ScalarMapping, range_sizes) = range_sizes
 
@@ -70,7 +102,6 @@
     B = ScalarMapping{Float64,1,1}(3.0)
 
     v = [1.1,1.2,1.3]
-
     for i ∈ eachindex(v)
         @test ((A+B)*v)[i] == 2*v[i] + 3*v[i]
     end
@@ -88,24 +119,45 @@
         data::T1
     end
     Base.size(v::DummyArray) = size(v.data)
-    Base.getindex(v::DummyArray, I...) = v.data[I...]
+    Base.getindex(v::DummyArray{T,D}, I::Vararg{Int,D}) where {T,D} = v.data[I...]
 
     # Test lazy operations
     v1 = [1, 2.3, 4]
     v2 = [1., 2, 3]
-    r_add = v1 .+ v2
-    r_sub = v1 .- v2
-    r_times = v1 .* v2
-    r_div = v1 ./ v2
+    s = 3.4
+    r_add_v = v1 .+ v2
+    r_sub_v = v1 .- v2
+    r_times_v = v1 .* v2
+    r_div_v = v1 ./ v2
+    r_add_s = v1 .+ s
+    r_sub_s = v1 .- s
+    r_times_s = v1 .* s
+    r_div_s = v1 ./ s
     @test isa(v1 +̃ v2, LazyArray)
     @test isa(v1 -̃ v2, LazyArray)
     @test isa(v1 *̃ v2, LazyArray)
     @test isa(v1 /̃ v2, LazyArray)
+    @test isa(v1 +̃ s, LazyArray)
+    @test isa(v1 -̃ s, LazyArray)
+    @test isa(v1 *̃ s, LazyArray)
+    @test isa(v1 /̃ s, LazyArray)
+    @test isa(s +̃ v1, LazyArray)
+    @test isa(s -̃ v1, LazyArray)
+    @test isa(s *̃ v1, LazyArray)
+    @test isa(s /̃ v1, LazyArray)
     for i ∈ eachindex(v1)
-        @test (v1 +̃ v2)[i] == r_add[i]
-        @test (v1 -̃ v2)[i] == r_sub[i]
-        @test (v1 *̃ v2)[i] == r_times[i]
-        @test (v1 /̃ v2)[i] == r_div[i]
+        @test (v1 +̃ v2)[i] == r_add_v[i]
+        @test (v1 -̃ v2)[i] == r_sub_v[i]
+        @test (v1 *̃ v2)[i] == r_times_v[i]
+        @test (v1 /̃ v2)[i] == r_div_v[i]
+        @test (v1 +̃ s)[i] == r_add_s[i]
+        @test (v1 -̃ s)[i] == r_sub_s[i]
+        @test (v1 *̃ s)[i] == r_times_s[i]
+        @test (v1 /̃ s)[i] == r_div_s[i]
+        @test (s +̃ v1)[i] == r_add_s[i]
+        @test (s -̃ v1)[i] == -r_sub_s[i]
+        @test (s *̃ v1)[i] == r_times_s[i]
+        @test (s /̃ v1)[i] == 1/r_div_s[i]
     end
     @test_throws BoundsError (v1 +̃  v2)[4]
     v2 = [1., 2, 3, 4]
@@ -120,8 +172,8 @@
     @test isa(v1 - v2, LazyArray)
     @test isa(v2 - v1, LazyArray)
     for i ∈ eachindex(v2)
-        @test (v1 + v2)[i] == (v2 + v1)[i] == r_add[i]
-        @test (v1 - v2)[i] == -(v2 - v1)[i] == r_sub[i]
+        @test (v1 + v2)[i] == (v2 + v1)[i] == r_add_v[i]
+        @test (v1 - v2)[i] == -(v2 - v1)[i] == r_sub_v[i]
     end
     @test_throws BoundsError (v1 + v2)[4]
     v2 = [1., 2, 3, 4]
--- a/Manifest.toml	Wed Jun 26 15:07:47 2019 +0200
+++ b/Manifest.toml	Mon Jun 22 21:43:05 2020 +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"
@@ -24,6 +24,7 @@
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
 [[LazyTensors]]
+deps = ["RegionIndices"]
 path = "LazyTensors"
 uuid = "62fbed2c-918d-11e9-279b-eb3a325b37d3"
 version = "0.1.0"
@@ -36,16 +37,16 @@
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
 [[OffsetArrays]]
-git-tree-sha1 = "1af2f79c7eaac3e019a0de41ef63335ff26a0a57"
+git-tree-sha1 = "87d0a91efe29352d5caaa271ae3927083c096e33"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "0.11.1"
+version = "0.11.4"
 
 [[Random]]
 deps = ["Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[RegionIndices]]
-path = "RegionIndices"
+path = "/Users/vidar/dev/hg/sbpteam/sbplib_julia/RegionIndices"
 uuid = "5d527584-97f1-11e9-084c-4540c7ecf219"
 version = "0.1.0"
 
--- a/RegionIndices/src/RegionIndices.jl	Wed Jun 26 15:07:47 2019 +0200
+++ b/RegionIndices/src/RegionIndices.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -30,6 +30,8 @@
 Base.convert(::Type{CartesianIndex}, I::NTuple{N,Index} where N) = CartesianIndex(convert.(Int, I))
 
 Base.Int(I::Index) = I.i
+Base.to_index(I::Index) = Int(I) #How to get this to work for all cases??
+Base.getindex(A::AbstractArray{T,N}, I::NTuple{N,Index}) where {T,N} = A[I...] #Is this ok??
 
 function Index(i::Integer, boundary_width::Integer, dim_size::Integer)
     return Index{getregion(i,boundary_width,dim_size)}(i)
@@ -65,6 +67,8 @@
     end
 end
 
+export getregion
+
 function getrange(gridsize::Integer, closuresize::Integer, region::DataType)
     if region == Lower
         r = 1:closuresize
--- a/SbpOperators/src/SbpOperators.jl	Wed Jun 26 15:07:47 2019 +0200
+++ b/SbpOperators/src/SbpOperators.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -3,161 +3,8 @@
 using RegionIndices
 
 include("stencil.jl")
-
-export D2, closureSize, apply, readOperator, apply_e, apply_d
-
-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(op::D2, v::AbstractVector, ::Type{Lower})
-    apply(op.eClosure,v,1)
-end
-
-function apply_e(op::D2, v::AbstractVector, ::Type{Upper})
-    apply(flip(op.eClosure),v,length(v))
-end
-
-
-function apply_d(op::D2, h_inv::Real, v::AbstractVector, ::Type{Lower})
-    -h_inv*apply(op.dClosure,v,1)
-end
-
-function apply_d(op::D2, h_inv::Real, v::AbstractVector, ::Type{Upper})
-    -h_inv*apply(flip(op.dClosure),v,length(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
+include("constantstenciloperator.jl")
+include("d2.jl")
+include("readoperator.jl")
 
 end # module
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/constantlaplace.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -0,0 +1,54 @@
+#TODO: Naming?! What is this? It is a 1D tensor operator but what is then the
+# potentially multi-D laplace tensor mapping then?
+# Ideally I would like the below to be the laplace operator in 1D, while the
+# multi-D operator is a a tuple of the 1D-operator. Possible via recursive
+# definitions? Or just bad design?
+"""
+    ConstantLaplaceOperator{T<:Real,N,M,K} <: TensorOperator{T,1}
+Implements the Laplace tensor operator `L` with constant grid spacing and coefficients
+in 1D dimension
+"""
+struct ConstantLaplaceOperator{T<:Real,N,M,K} <: TensorOperator{T,1}
+    h_inv::T # The grid spacing could be included in the stencil already. Preferable?
+    a::T # TODO: Better name?
+    innerStencil::Stencil{T,N}
+    closureStencils::NTuple{M,Stencil{T,K}}
+    parity::Parity
+    #TODO: Write a nice constructor
+end
+
+@enum Parity begin
+    odd = -1
+    even = 1
+end
+
+LazyTensors.domain_size(L::ConstantLaplaceOperator, range_size::NTuple{1,Integer}) = range_size
+
+function LazyTensors.apply(L::ConstantLaplaceOperator{T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
+    return apply(L, v, I[1])
+end
+
+# Apply for different regions Lower/Interior/Upper or Unknown region
+@inline function LazyTensors.apply(L::ConstantLaplaceOperator, v::AbstractVector, i::Index{Lower})
+    return @inbounds L.a*L.h_inv*L.h_inv*apply_stencil(L.closureStencils[Int(i)], v, Int(i))
+end
+
+@inline function LazyTensors.apply(L::ConstantLaplaceOperator, v::AbstractVector, i::Index{Interior})
+    return @inbounds L.a*L.h_inv*L.h_inv*apply_stencil(L.innerStencil, v, Int(i))
+end
+
+@inline function LazyTensors.apply(L::ConstantLaplaceOperator, v::AbstractVector, i::Index{Upper})
+    N = length(v) # TODO: Use domain_size here instead?
+    return @inbounds L.a*L.h_inv*L.h_inv*Int(L.parity)*apply_stencil_backwards(L.closureStencils[N-Int(i)+1], v, Int(i))
+end
+
+@inline function LazyTensors.apply(L::ConstantLaplaceOperator, v::AbstractVector, index::Index{Unknown})
+    N = length(v)  # TODO: Use domain_size here instead?
+    r = getregion(Int(index), closuresize(L), N)
+    i = Index(Int(index), r)
+    return apply(L, v, i)
+end
+
+function closuresize(L::ConstantLaplaceOperator{T<:Real,N,M,K}) where T,N,M,K
+    return M
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/constantstenciloperator.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -0,0 +1,109 @@
+abstract type ConstantStencilOperator end
+
+# Apply for different regions Lower/Interior/Upper or Unknown region
+@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Lower})
+    return @inbounds h_inv*h_inv*apply_stencil(op.closureStencils[Int(i)], v, Int(i))
+end
+
+@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Interior})
+    return @inbounds h_inv*h_inv*apply_stencil(op.innerStencil, v, Int(i))
+end
+
+@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Upper})
+    N = length(v)
+    return @inbounds h_inv*h_inv*Int(op.parity)*apply_stencil_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))
+end
+
+@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, index::Index{Unknown})
+    N = length(v)
+    r = getregion(Int(index), closuresize(op), N)
+    i = Index(Int(index), r)
+    return apply_2nd_derivative(op, h_inv, v, i)
+end
+export apply_2nd_derivative
+
+apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Lower}, N::Integer) where T = v*h*op.quadratureClosure[Int(i)]
+apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Upper}, N::Integer) where T = v*h*op.quadratureClosure[N-Int(i)+1]
+apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Interior}, N::Integer) where T = v*h
+
+function apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, index::Index{Unknown}, N::Integer) where T
+    r = getregion(Int(index), closuresize(op), N)
+    i = Index(Int(index), r)
+    return apply_quadrature(op, h, v, i, N)
+end
+export apply_quadrature
+
+# TODO: Evaluate if divisions affect performance
+apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Index{Lower}, N::Integer) where T = h_inv*v/op.quadratureClosure[Int(i)]
+apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Index{Upper}, N::Integer) where T = h_inv*v/op.quadratureClosure[N-Int(i)+1]
+apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Index{Interior}, N::Integer) where T = v*h_inv
+
+function apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, index::Index{Unknown}, N::Integer) where T
+    r = getregion(Int(index), closuresize(op), N)
+    i = Index(Int(index), r)
+    return apply_inverse_quadrature(op, h_inv, v, i, N)
+end
+
+export apply_inverse_quadrature
+
+function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Lower})
+    @boundscheck if length(v) < closuresize(op)
+        throw(BoundsError())
+    end
+    apply_stencil(op.eClosure,v,1)
+end
+
+function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Upper})
+    @boundscheck if length(v) < closuresize(op)
+        throw(BoundsError())
+    end
+    apply_stencil_backwards(op.eClosure,v,length(v))
+end
+export apply_boundary_value_transpose
+
+function apply_boundary_value(op::ConstantStencilOperator, v::Number, i::Index, N::Integer, ::Type{Lower})
+    @boundscheck if !(0<length(Int(i)) <= N)
+        throw(BoundsError())
+    end
+    op.eClosure[Int(i)-1]*v
+end
+
+function apply_boundary_value(op::ConstantStencilOperator, v::Number, i::Index,  N::Integer, ::Type{Upper})
+    @boundscheck if !(0<length(Int(i)) <= N)
+        throw(BoundsError())
+    end
+    op.eClosure[N-Int(i)]*v
+end
+export apply_boundary_value
+
+function apply_normal_derivative_transpose(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, ::Type{Lower})
+    @boundscheck if length(v) < closuresize(op)
+        throw(BoundsError())
+    end
+    h_inv*apply_stencil(op.dClosure,v,1)
+end
+
+function apply_normal_derivative_transpose(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, ::Type{Upper})
+    @boundscheck if length(v) < closuresize(op)
+        throw(BoundsError())
+    end
+    -h_inv*apply_stencil_backwards(op.dClosure,v,length(v))
+end
+
+export apply_normal_derivative_transpose
+
+function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, i::Index, N::Integer, ::Type{Lower})
+    @boundscheck if !(0<length(Int(i)) <= N)
+        throw(BoundsError())
+    end
+    h_inv*op.dClosure[Int(i)-1]*v
+end
+
+function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, i::Index, N::Integer, ::Type{Upper})
+    @boundscheck if !(0<length(Int(i)) <= N)
+        throw(BoundsError())
+    end
+    -h_inv*op.dClosure[N-Int(i)]*v
+end
+
+export apply_normal_derivative
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/d2.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -0,0 +1,19 @@
+export D2, closuresize, readOperator
+
+@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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/readoperator.jl	Mon Jun 22 21:43:05 2020 +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
--- a/SbpOperators/src/stencil.jl	Wed Jun 26 15:07:47 2019 +0200
+++ b/SbpOperators/src/stencil.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -21,7 +21,7 @@
     return s.weights[1 + i - s.range[1]]
 end
 
-Base.@propagate_inbounds @inline function apply(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
+Base.@propagate_inbounds @inline function apply_stencil(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
     w = s.weights[1]*v[i + s.range[1]]
     @simd for k ∈ 2:N
         w += s.weights[k]*v[i + s.range[1] + k-1]
@@ -29,7 +29,7 @@
     return w
 end
 
-Base.@propagate_inbounds @inline function apply_backwards(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
+Base.@propagate_inbounds @inline function apply_stencil_backwards(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
     w = s.weights[N]*v[i - s.range[2]]
     @simd for k ∈ N-1:-1:1
         w += s.weights[k]*v[i - s.range[1] - k + 1]
--- a/SbpOperators/test/runtests.jl	Wed Jun 26 15:07:47 2019 +0200
+++ b/SbpOperators/test/runtests.jl	Mon Jun 22 21:43:05 2020 +0200
@@ -1,4 +1,23 @@
 using SbpOperators
 using Test
 
-@test_broken false
+@testset "apply_quadrature" begin
+    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    h = 0.5
+
+    @test apply_quadrature(op, h, 1.0, 10, 100) == h
+
+    N = 10
+    qc = op.quadratureClosure
+    q = h.*(qc..., ones(N-2*closuresize(op))..., reverse(qc)...)
+    @assert length(q) == N
+
+    for i ∈ 1:N
+        @test apply_quadrature(op, h, 1.0, i, N) == q[i]
+    end
+
+    v = [2.,3.,2.,4.,5.,4.,3.,4.,5.,4.5]
+    for i ∈ 1:N
+        @test apply_quadrature(op, h, v[i], i, N) == q[i]*v[i]
+    end
+end
--- a/TODO.txt	Wed Jun 26 15:07:47 2019 +0200
+++ b/TODO.txt	Mon Jun 22 21:43:05 2020 +0200
@@ -8,7 +8,26 @@
 
 Profilera
 
-Konvertera till paket
 Skriv tester
 
 Specificera operatorer i TOML eller något liknande?
+
+Borde det finns motsvarande apply_stencil för apply_quadrature,
+apply_boundary_value och apply_normal_derivative?
+
+Borde man alltid skicka in N som parameter i apply_2nd_derivative, t.ex som i
+apply_quadrature?
+
+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)