Mercurial > repos > public > sbplib_julia

changeset 231:fbabfd4e8f20

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in boundary_conditions
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 26 Jun 2019 15:07:47 +0200
parents ce56727e4232 (current diff) 3d94f60f6fae (diff)
children 0f94dc29c4bf
files AbstractGrid.jl EquidistantGrid.jl d2_2nd.txt d2_4th.txt diffOp.jl grid.jl h_2nd.txt h_4th.txt index.jl sbpD2.jl stencil.jl
diffstat 44 files changed, 1302 insertions(+), 516 deletions(-) [+]
line wrap: on
line diff
--- a/AbstractGrid.jl	Fri Feb 22 15:22:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,15 +0,0 @@
-abstract type AbstractGrid end
-
-function dimension(grid::AbstractGrid)
-    error("Not implemented for abstact type AbstractGrid")
-end
-
-function points(grid::AbstractGrid)
-    error("Not implemented for abstact type AbstractGrid")
-end
-
-# Evaluate function f on the grid g
-function evalOn(g::AbstractGrid, f::Function)
-    F(x) = f(x...)
-    return F.(points(g))
-end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DiffOps/Manifest.toml	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,61 @@
+# This file is machine-generated - editing it directly is not advised
+
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[Grids]]
+deps = ["RegionIndices"]
+path = "../Grids"
+uuid = "960fdf28-97ed-11e9-2a74-bd90bf2fab5a"
+version = "0.1.0"
+
+[[InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[OffsetArrays]]
+git-tree-sha1 = "1af2f79c7eaac3e019a0de41ef63335ff26a0a57"
+uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+version = "0.11.1"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[RegionIndices]]
+path = "../RegionIndices"
+uuid = "5d527584-97f1-11e9-084c-4540c7ecf219"
+version = "0.1.0"
+
+[[SbpOperators]]
+deps = ["RegionIndices"]
+path = "../SbpOperators"
+uuid = "204357d8-97fd-11e9-05e7-010897a14cd0"
+version = "0.1.0"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[Test]]
+deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[TiledIteration]]
+deps = ["OffsetArrays", "Test"]
+git-tree-sha1 = "58f6f07d3b54a363ec283a8f5fc9fb4ecebde656"
+uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
+version = "0.2.3"
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DiffOps/Project.toml	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,16 @@
+name = "DiffOps"
+uuid = "39474f48-97ec-11e9-01fc-6ddcbe5918df"
+authors = ["Jonatan Werpers <jonatan.werpers@it.uu.se>"]
+version = "0.1.0"
+
+[deps]
+Grids = "960fdf28-97ed-11e9-2a74-bd90bf2fab5a"
+RegionIndices = "5d527584-97f1-11e9-084c-4540c7ecf219"
+SbpOperators = "204357d8-97fd-11e9-05e7-010897a14cd0"
+TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DiffOps/src/DiffOps.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,103 @@
+module DiffOps
+
+using RegionIndices
+using SbpOperators
+using Grids
+
+"""
+    DiffOp
+
+Supertype of differential operator discretisations.
+The action of the DiffOp is defined in the method
+    apply(D::DiffOp, v::AbstractVector, I...)
+"""
+abstract type DiffOp end
+
+function apply end
+
+function matrixRepresentation(D::DiffOp)
+    error("not implemented")
+end
+
+abstract type DiffOpCartesian{Dim} <: DiffOp end
+
+# DiffOp must have a grid of dimension Dim!!!
+function apply!(D::DiffOpCartesian{Dim}, u::AbstractArray{T,Dim}, v::AbstractArray{T,Dim}) where {T,Dim}
+    for I ∈ eachindex(D.grid)
+        u[I] = apply(D, v, I)
+    end
+
+    return nothing
+end
+export apply!
+
+function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
+    apply_region!(D, u, v, Lower, Lower)
+    apply_region!(D, u, v, Lower, Interior)
+    apply_region!(D, u, v, Lower, Upper)
+    apply_region!(D, u, v, Interior, Lower)
+    apply_region!(D, u, v, Interior, Interior)
+    apply_region!(D, u, v, Interior, Upper)
+    apply_region!(D, u, v, Upper, Lower)
+    apply_region!(D, u, v, Upper, Interior)
+    apply_region!(D, u, v, Upper, Upper)
+    return nothing
+end
+
+# 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))
+        @inbounds indextuple = (Index{r1}(I[1]), Index{r2}(I[2]))
+        @inbounds u[I] = apply(D, v, indextuple)
+    end
+    return nothing
+end
+export apply_region!
+
+function apply_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
+    apply_region_tiled!(D, u, v, Lower, Lower)
+    apply_region_tiled!(D, u, v, Lower, Interior)
+    apply_region_tiled!(D, u, v, Lower, Upper)
+    apply_region_tiled!(D, u, v, Interior, Lower)
+    apply_region_tiled!(D, u, v, Interior, Interior)
+    apply_region_tiled!(D, u, v, Interior, Upper)
+    apply_region_tiled!(D, u, v, Upper, Lower)
+    apply_region_tiled!(D, u, v, Upper, Interior)
+    apply_region_tiled!(D, u, v, Upper, Upper)
+    return nothing
+end
+
+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))
+    # TODO: Pass Tilesize to function
+    for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2))
+        for j ∈ tileaxs[2], i ∈ tileaxs[1]
+            I = ri[i,j]
+            u[I] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2])))
+        end
+    end
+    return nothing
+end
+export apply_region_tiled!
+
+function apply(D::DiffOp, v::AbstractVector)::AbstractVector
+    u = zeros(eltype(v), size(v))
+    apply!(D,v,u)
+    return u
+end
+
+export apply
+
+"""
+A BoundaryCondition should implement the method
+    sat(::DiffOp, v::AbstractArray, data::AbstractArray, ...)
+"""
+abstract type BoundaryCondition end
+
+
+include("laplace.jl")
+export Laplace
+
+
+end # module
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DiffOps/src/laplace.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,131 @@
+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
+
+struct BoundaryValue{N,M,K}
+	op::D2{Float64,N,M,K}
+	grid::EquidistantGrid
+	bId::CartesianBoundary
+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
+
+	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
+
+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
+    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)
+    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ᵢ
+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)
+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}
+    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
+end
+
+struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition
+    tau::Float64
+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)
+    # 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) +
+# 		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) +
+# 		sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau),  v, s.g_n, i)
+# end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DiffOps/test/runtests.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,4 @@
+using DiffOps
+using Test
+
+@test_broken false
--- a/EquidistantGrid.jl	Fri Feb 22 15:22:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,56 +0,0 @@
-# EquidistantGrid is a grid with equidistant grid spacing per coordinat
-# direction. The domain is defined through the two points P1 = x̄₁, P2 = x̄₂
-# by the exterior product of the vectors obtained by projecting (x̄₂-x̄₁) onto
-# the coordinate directions. E.g for a 2D grid with x̄₁=(-1,0) and x̄₂=(1,2)
-# the domain is defined as (-1,1)x(0,2).
-
-struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid
-    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
-
-    # 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)
-        return new{Dim,T}(size, limit_lower, limit_upper, inverse_spacing)
-    end
-end
-
-function Base.eachindex(grid::EquidistantGrid)
-    CartesianIndices(grid.size)
-end
-
-# Returns the number of dimensions of an EquidistantGrid.
-#
-# @Input: grid - an EquidistantGrid
-# @Return: dimension - The dimension of the grid
-function dimension(grid::EquidistantGrid)
-    return length(grid.size)
-end
-
-# Returns the spacing of the grid
-#
-function spacing(grid::EquidistantGrid)
-    return 1.0./grid.inverse_spacing
-end
-
-# Computes the points of an EquidistantGrid as an array of tuples with
-# the same dimension as the grid.
-#
-# @Input: grid - an EquidistantGrid
-# @Return: points - the points of the grid.
-function points(grid::EquidistantGrid)
-    # TODO: Make this return an abstract array?
-    indices = Tuple.(CartesianIndices(grid.size))
-    h = spacing(grid)
-    return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices)
-end
-
-function pointsalongdim(grid::EquidistantGrid, dim::Integer)
-    @assert dim<=dimension(grid)
-    @assert dim>0
-    points = collect(range(grid.limit_lower[dim],stop=grid.limit_upper[dim],length=grid.size[dim]))
-end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Grids/Manifest.toml	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +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"
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Grids/Project.toml	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,13 @@
+name = "Grids"
+uuid = "960fdf28-97ed-11e9-2a74-bd90bf2fab5a"
+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"
+
+[targets]
+test = ["Test"]
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Grids/src/AbstractGrid.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,24 @@
+"""
+     AbstractGrid
+
+Should implement
+    dimension(grid::AbstractGrid)
+    points(grid::AbstractGrid)
+
+"""
+abstract type AbstractGrid end
+
+function dimension end
+function points end
+export dimension, points
+
+"""
+    evalOn(g::AbstractGrid, f::Function)
+
+Evaluate function f on the grid g
+"""
+function evalOn(g::AbstractGrid, f::Function)
+    F(x) = f(x...)
+    return F.(points(g))
+end
+export evalOn
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Grids/src/EquidistantGrid.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,58 @@
+# EquidistantGrid is a grid with equidistant grid spacing per coordinat
+# direction. The domain is defined through the two points P1 = x̄₁, P2 = x̄₂
+# by the exterior product of the vectors obtained by projecting (x̄₂-x̄₁) onto
+# the coordinate directions. E.g for a 2D grid with x̄₁=(-1,0) and x̄₂=(1,2)
+# the domain is defined as (-1,1)x(0,2).
+
+export EquidistantGrid
+
+struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid
+    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
+
+    # 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)
+        return new{Dim,T}(size, limit_lower, limit_upper, inverse_spacing)
+    end
+end
+
+function Base.eachindex(grid::EquidistantGrid)
+    CartesianIndices(grid.size)
+end
+
+# Returns the number of dimensions of an EquidistantGrid.
+#
+# @Input: grid - an EquidistantGrid
+# @Return: dimension - The dimension of the grid
+function dimension(grid::EquidistantGrid)
+    return length(grid.size)
+end
+
+# Returns the spacing of the grid
+#
+function spacing(grid::EquidistantGrid)
+    return 1.0./grid.inverse_spacing
+end
+
+# Computes the points of an EquidistantGrid as an array of tuples with
+# the same dimension as the grid.
+#
+# @Input: grid - an EquidistantGrid
+# @Return: points - the points of the grid.
+function points(grid::EquidistantGrid)
+    # TODO: Make this return an abstract array?
+    indices = Tuple.(CartesianIndices(grid.size))
+    h = spacing(grid)
+    return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices)
+end
+
+function pointsalongdim(grid::EquidistantGrid, dim::Integer)
+    @assert dim<=dimension(grid)
+    @assert dim>0
+    points = collect(range(grid.limit_lower[dim],stop=grid.limit_upper[dim],length=grid.size[dim]))
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Grids/src/Grids.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,15 @@
+module Grids
+
+using RegionIndices
+
+export BoundaryIdentifier, CartesianBoundary
+
+abstract type BoundaryIdentifier end
+struct CartesianBoundary{Dim, R<:Region} <: BoundaryIdentifier end
+dim(::CartesianBoundary{Dim, R}) where {Dim, R} = Dim
+region(::CartesianBoundary{Dim, R}) where {Dim, R} = R
+
+include("AbstractGrid.jl")
+include("EquidistantGrid.jl")
+
+end # module
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Grids/test/runtests.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,4 @@
+using Grids
+using Test
+
+@test_broken false
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/LazyTensors/Manifest.toml	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,2 @@
+# This file is machine-generated - editing it directly is not advised
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/LazyTensors/Project.toml	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,10 @@
+name = "LazyTensors"
+uuid = "62fbed2c-918d-11e9-279b-eb3a325b37d3"
+authors = ["Jonatan Werpers <jonatan.werpers@it.uu.se>"]
+version = "0.1.0"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/LazyTensors/src/LazyTensors.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,6 @@
+module LazyTensors
+
+include("tensor_mapping.jl")
+include("lazy_operations.jl")
+
+end # module
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/LazyTensors/src/lazy_operations.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,193 @@
+"""
+    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/tensor_mapping.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,78 @@
+"""
+    TensorMapping{T,R,D}
+
+Describes a mapping of a D dimension tensor to an R dimension tensor.
+The action of the mapping is implemented through the method
+
+    apply(t::TensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {R,D,T}
+
+The size of range tensor should be dependent on the size of the domain tensor
+and the type should implement the methods
+
+    range_size(::TensorMapping{T,R,D}, domain_size::NTuple{D,Integer}) where {T,R,D}
+    domain_size(::TensorMapping{T,R,D}, range_size::NTuple{R,Integer}) where {T,R,D}
+
+to allow querying for one or the other.
+
+Optionally the action of the transpose may be defined through
+    apply_transpose(t::TensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {R,D,T}
+"""
+abstract type TensorMapping{T,R,D} end
+export TensorMapping
+
+"""
+    apply(t::TensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {R,D,T}
+
+Return the result of the mapping for a given index.
+"""
+function apply end
+export apply
+
+"""
+    apply_transpose(t::TensorMapping{T,R,D}, v::AbstractArray{T,R}, I::Vararg) where {R,D,T}
+
+Return the result of the transposed mapping for a given index.
+"""
+function apply_transpose end
+export apply_transpose
+
+"""
+Return the dimension of the range space of a given mapping
+"""
+range_dim(::TensorMapping{T,R,D}) where {T,R,D} = R
+
+"""
+Return the dimension of the domain space of a given mapping
+"""
+domain_dim(::TensorMapping{T,R,D}) where {T,R,D} = D
+
+export range_dim, domain_dim
+
+"""
+    range_size(M::TensorMapping, domain_size)
+
+Return the resulting range size for the mapping applied to a given domain_size
+"""
+function range_size end
+
+"""
+    domain_size(M::TensorMapping, range_size)
+
+Return the resulting domain size for the mapping applied to a given range_size
+"""
+function domain_size end
+
+export range_size, domain_size
+# TODO: Think about boundschecking!
+
+
+"""
+    TensorOperator{T,D}
+
+A `TensorMapping{T,D,D}` where the range and domain tensor have the same number of
+dimensions and the same size.
+"""
+abstract type TensorOperator{T,D} <: TensorMapping{T,D,D} end
+export TensorOperator
+domain_size(::TensorOperator{T,D}, range_size::NTuple{D,Integer}) where {T,D} = range_size
+range_size(::TensorOperator{T,D}, domain_size::NTuple{D,Integer}) where {T,D} = domain_size
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/LazyTensors/test/runtests.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,130 @@
+using Test
+using LazyTensors
+
+@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
+    @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
+end
+
+@testset "Generic Operator methods" begin
+    struct DummyOperator{T,D} <: TensorOperator{T,D} end
+    @test range_size(DummyOperator{Int,2}(), (3,5)) == (3,5)
+    @test domain_size(DummyOperator{Float64, 3}(), (3,3,1)) == (3,3,1)
+end
+
+@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.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
+
+    m = DummyMapping{Float64,2,3}()
+    @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 range_size(m', (0,0)) == :domain_size
+    @test domain_size(m', (0,0,0)) == :range_size
+end
+
+@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
+
+
+    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*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_broken BoundsError == (m*m*v)[0]
+    @test_broken BoundsError == (m*m*v)[7]
+end
+
+@testset "TensorMapping binary operations" begin
+    struct ScalarMapping{T,R,D} <: TensorMapping{T,R,D}
+        λ::T
+    end
+
+    LazyTensors.apply(m::ScalarMapping{T,R,D}, v, i) 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
+
+    A = ScalarMapping{Float64,1,1}(2.0)
+    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
+
+    for i ∈ eachindex(v)
+        @test ((A-B)*v)[i] == 2*v[i] - 3*v[i]
+    end
+
+    @test range_size(A+B, (3,)) == range_size(A, (3,)) == range_size(B,(3,))
+    @test domain_size(A+B, (3,)) == domain_size(A, (3,)) == domain_size(B,(3,))
+end
+
+@testset "LazyArray" begin
+    struct DummyArray{T,D, T1<:AbstractArray{T,D}} <: LazyArray{T,D}
+        data::T1
+    end
+    Base.size(v::DummyArray) = size(v.data)
+    Base.getindex(v::DummyArray, I...) = 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
+    @test isa(v1 +̃ v2, LazyArray)
+    @test isa(v1 -̃ v2, LazyArray)
+    @test isa(v1 *̃ v2, LazyArray)
+    @test isa(v1 /̃ v2, 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]
+    end
+    @test_throws BoundsError (v1 +̃  v2)[4]
+    v2 = [1., 2, 3, 4]
+    # Test that size of arrays is asserted when not specified inbounds
+    @test_throws DimensionMismatch v1 +̃ v2
+
+    # Test operations on LazyArray
+    v1 = DummyArray([1, 2.3, 4])
+    v2 = [1., 2, 3]
+    @test isa(v1 + v2, LazyArray)
+    @test isa(v2 + v1, LazyArray)
+    @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]
+    end
+    @test_throws BoundsError (v1 + v2)[4]
+    v2 = [1., 2, 3, 4]
+    # Test that size of arrays is asserted when not specified inbounds
+    @test_throws DimensionMismatch v1 + v2
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Manifest.toml	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,72 @@
+# This file is machine-generated - editing it directly is not advised
+
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[DiffOps]]
+deps = ["Grids", "RegionIndices", "SbpOperators", "TiledIteration"]
+path = "DiffOps"
+uuid = "39474f48-97ec-11e9-01fc-6ddcbe5918df"
+version = "0.1.0"
+
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[Grids]]
+deps = ["RegionIndices"]
+path = "Grids"
+uuid = "960fdf28-97ed-11e9-2a74-bd90bf2fab5a"
+version = "0.1.0"
+
+[[InteractiveUtils]]
+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"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[OffsetArrays]]
+git-tree-sha1 = "1af2f79c7eaac3e019a0de41ef63335ff26a0a57"
+uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+version = "0.11.1"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[RegionIndices]]
+path = "RegionIndices"
+uuid = "5d527584-97f1-11e9-084c-4540c7ecf219"
+version = "0.1.0"
+
+[[SbpOperators]]
+deps = ["RegionIndices"]
+path = "SbpOperators"
+uuid = "204357d8-97fd-11e9-05e7-010897a14cd0"
+version = "0.1.0"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[Test]]
+deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[TiledIteration]]
+deps = ["OffsetArrays", "Test"]
+git-tree-sha1 = "58f6f07d3b54a363ec283a8f5fc9fb4ecebde656"
+uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
+version = "0.2.3"
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Project.toml	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,6 @@
+[deps]
+DiffOps = "39474f48-97ec-11e9-01fc-6ddcbe5918df"
+Grids = "960fdf28-97ed-11e9-2a74-bd90bf2fab5a"
+LazyTensors = "62fbed2c-918d-11e9-279b-eb3a325b37d3"
+RegionIndices = "5d527584-97f1-11e9-084c-4540c7ecf219"
+SbpOperators = "204357d8-97fd-11e9-05e7-010897a14cd0"
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/RegionIndices/Project.toml	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,10 @@
+name = "RegionIndices"
+uuid = "5d527584-97f1-11e9-084c-4540c7ecf219"
+authors = ["Jonatan Werpers <jonatan.werpers@it.uu.se>"]
+version = "0.1.0"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/RegionIndices/src/RegionIndices.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,81 @@
+module RegionIndices
+
+abstract type Region end
+struct Interior <: Region end
+struct Lower    <: Region end
+struct Upper    <: Region end
+struct Unknown  <: Region end
+
+export Region, Interior, Lower, Upper, Unknown
+
+struct Index{R<:Region, T<:Integer}
+    i::T
+
+    Index{R,T}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i)
+    Index{R}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i)
+    Index(i::T, ::Type{R}) where {R<:Region,T<:Integer} = Index{R,T}(i)
+    Index(t::Tuple{T, DataType}) where {R<:Region,T<:Integer} = Index{t[2],T}(t[1]) # TBD: This is not very specific in what types are allowed in t[2]. Can this be fixed?
+end
+
+export Index
+
+# Index(R::Type{<:Region}) = Index{R}
+
+## Vill kunna skriva
+## IndexTupleType(Int, (Lower, Interior))
+Index(R::Type{<:Region}, T::Type{<:Integer}) = Index{R,T}
+IndexTupleType(T::Type{<:Integer},R::NTuple{N, DataType} where N) = Tuple{Index.(R, T)...}
+
+Base.convert(::Type{T}, i::Index{R,T} where R) where T = i.i
+Base.convert(::Type{CartesianIndex}, I::NTuple{N,Index} where N) = CartesianIndex(convert.(Int, I))
+
+Base.Int(I::Index) = I.i
+
+function Index(i::Integer, boundary_width::Integer, dim_size::Integer)
+    return Index{getregion(i,boundary_width,dim_size)}(i)
+end
+
+IndexTuple(t::Vararg{Tuple{T, DataType}}) where T<:Integer = Index.(t)
+export IndexTuple
+
+# TODO: Use the values of the region structs, e.g. Lower(), for the region parameter instead of the types.
+# For example the following works:
+#   (Lower(),Upper()) isa NTuple{2, Region} -> true
+#   typeof((Lower(),Upper()))               -> Tuple{Lower,Upper}
+function regionindices(gridsize::NTuple{Dim,Integer}, closuresize::Integer, region::NTuple{Dim,DataType}) where Dim
+    return regionindices(gridsize, ntuple(x->closuresize,Dim), region)
+end
+
+function regionindices(gridsize::NTuple{Dim,Integer}, closuresize::NTuple{Dim,Integer}, region::NTuple{Dim,DataType}) where Dim
+    regions = map(getrange,gridsize,closuresize,region)
+    return CartesianIndices(regions)
+end
+
+export regionindices
+
+function getregion(i::Integer, boundary_width::Integer, dim_size::Integer)
+	if 0 < i <= boundary_width
+        return Lower
+    elseif boundary_width < i <= dim_size-boundary_width
+        return Interior
+    elseif dim_size-boundary_width < i <= dim_size
+        return Upper
+    else
+        error("Bounds error") # TODO: Make this more standard
+    end
+end
+
+function getrange(gridsize::Integer, closuresize::Integer, region::DataType)
+    if region == Lower
+        r = 1:closuresize
+    elseif region == Interior
+        r = (closuresize+1):(gridsize - closuresize)
+    elseif region == Upper
+        r = (gridsize - closuresize + 1):gridsize
+    else
+        error("Unspecified region")
+    end
+    return r
+end
+
+end # module
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/RegionIndices/test/runtests.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,4 @@
+using RegionIndices
+using Test
+
+@test_broken false
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/Manifest.toml	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +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"
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/Project.toml	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,13 @@
+name = "SbpOperators"
+uuid = "204357d8-97fd-11e9-05e7-010897a14cd0"
+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"
+
+[targets]
+test = ["Test"]
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/operators/d2_2nd.txt	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,13 @@
+# D2 order 2
+
+boundary_stencils
+1 -2 1
+
+inner_stencil
+1 -2 1
+
+e
+1
+
+d1
+-3/2 2 -1/2
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/operators/d2_4th.txt	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,16 @@
+# D2 order 4
+
+boundary_stencils
+     2    -5       4      -1     0     0
+     1    -2       1       0     0     0
+ -4/43 59/43 -110/43   59/43 -4/43     0
+ -1/49     0   59/49 -118/49 64/49 -4/49
+
+inner_stencil
+-1/12     4/3  -5/2   4/3 -1/12
+
+e
+1
+
+d1
+-11/6 3 -3/2 1/3
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/operators/h_2nd.txt	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,7 @@
+# H order 2
+
+closure
+1/2
+
+inner_stencil
+1
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/operators/h_4th.txt	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,7 @@
+# H order 4
+
+closure
+17/48 59/48 43/48 49/48
+
+inner_stencil
+1
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/SbpOperators.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,163 @@
+module SbpOperators
+
+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
+
+end # module
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/stencil.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,38 @@
+struct Stencil{T<:Real,N}
+    range::Tuple{Int,Int}
+    weights::NTuple{N,T}
+
+    function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T <: Real, N}
+        @assert range[2]-range[1]+1 == N
+        new{T,N}(range,weights)
+    end
+end
+
+function flip(s::Stencil)
+    range = (-s.range[2], -s.range[1])
+    return Stencil(range, reverse(s.weights))
+end
+
+# Provides index into the Stencil based on offset for the root element
+@inline function Base.getindex(s::Stencil, i::Int)
+    @boundscheck if i < s.range[1] || s.range[2] < i
+        return eltype(s.weights)(0)
+    end
+    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}
+    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]
+    end
+    return w
+end
+
+Base.@propagate_inbounds @inline function apply_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]
+    end
+    return w
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/test/runtests.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -0,0 +1,4 @@
+using SbpOperators
+using Test
+
+@test_broken false
--- a/TODO.txt	Fri Feb 22 15:22:34 2019 +0100
+++ b/TODO.txt	Wed Jun 26 15:07:47 2019 +0200
@@ -10,3 +10,5 @@
 
 Konvertera till paket
 Skriv tester
+
+Specificera operatorer i TOML eller något liknande?
--- a/d2_2nd.txt	Fri Feb 22 15:22:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,13 +0,0 @@
-# D2 order 2
-
-boundary_stencils
-1 -2 1
-
-inner_stencil
-1 -2 1
-
-e
-1
-
-d1
--3/2 2 -1/2
--- a/d2_4th.txt	Fri Feb 22 15:22:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,16 +0,0 @@
-# D2 order 4
-
-boundary_stencils
-     2    -5       4      -1     0     0
-     1    -2       1       0     0     0
- -4/43 59/43 -110/43   59/43 -4/43     0
- -1/49     0   59/49 -118/49 64/49 -4/49
-
-inner_stencil
--1/12     4/3  -5/2   4/3 -1/12
-
-e
-1
-
-d1
--11/6 3 -3/2 1/3
\ No newline at end of file
--- a/diffOp.jl	Fri Feb 22 15:22:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,131 +0,0 @@
-abstract type DiffOp end
-
-# TBD: The "error("not implemented")" thing seems to be hiding good error information. How to fix that? Different way of saying that these should be implemented?
-function apply(D::DiffOp, v::AbstractVector, i::Int)
-    error("not implemented")
-end
-
-function innerProduct(D::DiffOp, u::AbstractVector, v::AbstractVector)::Real
-    error("not implemented")
-end
-
-function matrixRepresentation(D::DiffOp)
-    error("not implemented")
-end
-
-function boundaryCondition(D::DiffOp,b::Grid.BoundaryId,type)::(Closure, Penalty)
-    error("not implemented")
-end
-
-function interface(Du::DiffOp, Dv::DiffOp, b::Grid.BoundaryId; type)
-    error("not implemented")
-end
-
-abstract type Closure end
-
-function apply(c::Closure, v::AbstractVector, i::Int)
-    error("not implemented")
-end
-
-abstract type Penalty end
-
-function apply(c::Penalty, g, i::Int)
-    error("not implemented")
-end
-
-abstract type DiffOpCartesian{Dim} <: DiffOp end
-
-# DiffOp must have a grid of dimension Dim!!!
-function apply!(D::DiffOpCartesian{Dim}, u::AbstractArray{T,Dim}, v::AbstractArray{T,Dim}) where {T,Dim}
-    for I ∈ eachindex(D.grid)
-        u[I] = apply(D, v, I)
-    end
-
-    return nothing
-end
-
-function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
-    apply_region!(D, u, v, Lower, Lower)
-    apply_region!(D, u, v, Lower, Interior)
-    apply_region!(D, u, v, Lower, Upper)
-    apply_region!(D, u, v, Interior, Lower)
-    apply_region!(D, u, v, Interior, Interior)
-    apply_region!(D, u, v, Interior, Upper)
-    apply_region!(D, u, v, Upper, Lower)
-    apply_region!(D, u, v, Upper, Interior)
-    apply_region!(D, u, v, Upper, Upper)
-    return nothing
-end
-
-# 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))
-        @inbounds indextuple = (Index{r1}(I[1]), Index{r2}(I[2]))
-        @inbounds u[I] = apply(D, v, indextuple)
-    end
-    return nothing
-end
-
-function apply_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
-    apply_region_tiled!(D, u, v, Lower, Lower)
-    apply_region_tiled!(D, u, v, Lower, Interior)
-    apply_region_tiled!(D, u, v, Lower, Upper)
-    apply_region_tiled!(D, u, v, Interior, Lower)
-    apply_region_tiled!(D, u, v, Interior, Interior)
-    apply_region_tiled!(D, u, v, Interior, Upper)
-    apply_region_tiled!(D, u, v, Upper, Lower)
-    apply_region_tiled!(D, u, v, Upper, Interior)
-    apply_region_tiled!(D, u, v, Upper, Upper)
-    return nothing
-end
-
-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))
-    # TODO: Pass Tilesize to function
-    for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2))
-        for j ∈ tileaxs[2], i ∈ tileaxs[1]
-            I = ri[i,j]
-            u[I] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2])))
-        end
-    end
-    return nothing
-end
-
-function apply(D::DiffOp, v::AbstractVector)::AbstractVector
-    u = zeros(eltype(v), size(v))
-    apply!(D,v,u)
-    return u
-end
-
-struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim}
-    grid::Grid.EquidistantGrid{Dim,T}
-    a::T
-    op::D2{Float64,N,M,K}
-end
-
-function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim
-    error("not implemented")
-end
-
-# u = L*v
-function apply(L::Laplace{1}, v::AbstractVector, i::Int)
-    uᵢ = L.a * apply(L.op, L.grid.spacing[1], v, i)
-    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*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*apply(L.op, L.grid.inverse_spacing[2], vy, I[2])
-    return uᵢ
-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)
-end
--- a/grid.jl	Fri Feb 22 15:22:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,9 +0,0 @@
-module Grid
-
-# TODO: Where is this used?
-abstract type BoundaryId end
-
-include("AbstractGrid.jl")
-include("EquidistantGrid.jl")
-
-end
--- a/h_2nd.txt	Fri Feb 22 15:22:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,7 +0,0 @@
-# H order 2
-
-closure
-1/2
-
-inner_stencil
-1
--- a/h_4th.txt	Fri Feb 22 15:22:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,7 +0,0 @@
-# H order 4
-
-closure
-17/48 59/48 43/48 49/48
-
-inner_stencil
-1
\ No newline at end of file
--- a/index.jl	Fri Feb 22 15:22:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,66 +0,0 @@
-abstract type Region end
-struct Interior <: Region end
-struct Lower    <: Region end
-struct Upper    <: Region end
-struct Unknown  <: Region end
-
-struct Index{R<:Region, T<:Integer}
-    i::T
-
-    Index{R,T}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i)
-    Index{R}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i)
-    Index(i::T, ::Type{R}) where {R<:Region,T<:Integer} = Index{R,T}(i)
-    Index(t::Tuple{T, DataType}) where {R<:Region,T<:Integer} = Index{t[2],T}(t[1]) # TBD: This is not very specific in what types are allowed in t[2]. Can this be fixed?
-end
-
-# Index(R::Type{<:Region}) = Index{R}
-
-## Vill kunna skriva
-## IndexTupleType(Int, (Lower, Interior))
-Index(R::Type{<:Region}, T::Type{<:Integer}) = Index{R,T}
-IndexTupleType(T::Type{<:Integer},R::NTuple{N, DataType} where N) = Tuple{Index.(R, T)...}
-
-Base.convert(::Type{T}, i::Index{R,T} where R) where T = i.i
-Base.convert(::Type{CartesianIndex}, I::NTuple{N,Index} where N) = CartesianIndex(convert.(Int, I))
-
-Base.Int(I::Index) = I.i
-
-function Index(i::Integer, boundary_width::Integer, dim_size::Integer)
-    if 0 < i <= boundary_width
-        return Index{Lower}(i)
-    elseif boundary_width < i <= dim_size-boundary_width
-        return Index{Interior}(i)
-    elseif dim_size-boundary_width < i <= dim_size
-        return Index{Upper}(i)
-    else
-        error("Bounds error") # TODO: Make this more standard
-    end
-end
-
-IndexTuple(t::Vararg{Tuple{T, DataType}}) where T<:Integer = Index.(t)
-
-# TODO: Use the values of the region structs, e.g. Lower(), for the region parameter instead of the types.
-# For example the following works:
-#   (Lower(),Upper()) isa NTuple{2, Region} -> true
-#   typeof((Lower(),Upper()))               -> Tuple{Lower,Upper}
-function regionindices(gridsize::NTuple{Dim,Integer}, closuresize::Integer, region::NTuple{Dim,DataType}) where Dim
-    return regionindices(gridsize, ntuple(x->closuresize,Dim), region)
-end
-
-function regionindices(gridsize::NTuple{Dim,Integer}, closuresize::NTuple{Dim,Integer}, region::NTuple{Dim,DataType}) where Dim
-    regions = map(getrange,gridsize,closuresize,region)
-    return CartesianIndices(regions)
-end
-
-function getrange(gridsize::Integer, closuresize::Integer, region::DataType)
-    if region == Lower
-        r = 1:closuresize
-    elseif region == Interior
-        r = (closuresize+1):(gridsize - closuresize)
-    elseif region == Upper
-        r = (gridsize - closuresize + 1):gridsize
-    else
-        error("Unspecified region")
-    end
-    return r
-end
--- a/sbp.jl	Fri Feb 22 15:22:34 2019 +0100
+++ b/sbp.jl	Wed Jun 26 15:07:47 2019 +0200
@@ -1,8 +1,9 @@
 module sbp
-include("index.jl")
-include("grid.jl")
-include("stencil.jl")
-include("sbpD2.jl")
-include("diffOp.jl")
+
+using Grids
+using RegionIndices
+using SbpOperators
+using DiffOps
+
 include("TimeStepper.jl")
 end  # module
--- a/sbpD2.jl	Fri Feb 22 15:22:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,153 +0,0 @@
-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::Real, v::AbstractVector, ::Type{Lower})
-    -apply(op.dClosure,v,1)/h
-end
-
-function apply_d(op::D2, h::Real, v::AbstractVector, ::Type{Upper})
-    -apply(flip(op.dClosure),v,length(v))/h
-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
--- a/stencil.jl	Fri Feb 22 15:22:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,38 +0,0 @@
-struct Stencil{T<:Real,N}
-    range::Tuple{Int,Int}
-    weights::NTuple{N,T}
-
-    function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T <: Real, N}
-        @assert range[2]-range[1]+1 == N
-        new{T,N}(range,weights)
-    end
-end
-
-function flip(s::Stencil)
-    range = (-s.range[2], -s.range[1])
-    return Stencil(range, reverse(s.weights))
-end
-
-# Provides index into the Stencil based on offset for the root element
-@inline function Base.getindex(s::Stencil, i::Int)
-    @boundscheck if i < s.range[1] || s.range[2] < i
-        return eltype(s.weights)(0)
-    end
-    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}
-    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]
-    end
-    return w
-end
-
-Base.@propagate_inbounds @inline function apply_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]
-    end
-    return w
-end