diff -r f01b70b81e95 -r 6b6d921e8f05 AbstractGrid.jl
--- a/AbstractGrid.jl	Thu Feb 07 16:34:55 2019 +0100
+++ b/AbstractGrid.jl	Thu Feb 21 14:19:25 2019 +0100
@@ -1,10 +1,6 @@
 abstract type AbstractGrid end
 
-function numberOfDimensions(grid::AbstractGrid)
-    error("Not implemented for abstact type AbstractGrid")
-end
-
-function numberOfPoints(grid::AbstractGrid)
+function dimension(grid::AbstractGrid)
     error("Not implemented for abstact type AbstractGrid")
 end
 

diff -r f01b70b81e95 -r 6b6d921e8f05 EquidistantGrid.jl
--- a/EquidistantGrid.jl	Thu Feb 07 16:34:55 2019 +0100
+++ b/EquidistantGrid.jl	Thu Feb 21 14:19:25 2019 +0100
@@ -5,85 +5,52 @@
 # the domain is defined as (-1,1)x(0,2).
 
 struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid
-    numberOfPointsPerDim::NTuple{Dim, Int} # First coordinate direction stored first, then
-
+    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(nPointsPerDim::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T
-        @assert all(nPointsPerDim.>0)
+    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)
-        return new{Dim,T}(nPointsPerDim, limit_lower, limit_upper)
+        inverse_spacing = (size.-1)./abs.(limit_upper.-limit_lower)
+        return new{Dim,T}(size, limit_lower, limit_upper, inverse_spacing)
     end
+end
 
-    # # 1D constructor which can be called as EquidistantGrid(m, (xl,xr))
-    # function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Real})
-    #     return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],)))
-    # end
-
+function Base.eachindex(grid::EquidistantGrid)
+    CartesianIndices(grid.size)
 end
 
 # Returns the number of dimensions of an EquidistantGrid.
 #
 # @Input: grid - an EquidistantGrid
-# @Return: numberOfPoints - The number of dimensions
-function numberOfDimensions(grid::EquidistantGrid)
-    return length(grid.numberOfPointsPerDim)
-end
-
-# Computes the total number of points of an EquidistantGrid.
-#
-# @Input: grid - an EquidistantGrid
-# @Return: numberOfPoints - The total number of points
-function numberOfPoints(grid::EquidistantGrid)
-    return prod(grid.numberOfPointsPerDim)
+# @Return: dimension - The dimension of the grid
+function dimension(grid::EquidistantGrid)
+    return length(grid.size)
 end
 
-# Computes the grid spacing of an EquidistantGrid, i.e the unsigned distance
-# between two points for each coordinate direction.
+# Returns the spacing of the grid
 #
-# @Input: grid - an EquidistantGrid
-# @Return: h̄ - Grid spacing for each coordinate direction stored in a tuple.
-function spacings(grid::EquidistantGrid)
-    return abs.(grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1)
+function spacing(grid::EquidistantGrid)
+    return 1.0./grid.inverse_spacing
 end
 
-function Base.eachindex(grid::EquidistantGrid)
-    CartesianIndices(grid.numberOfPointsPerDim)
-end
-
-# Computes the points of an EquidistantGrid as a vector of tuples. The vector is ordered
-# such that points in the first coordinate direction varies first, then the second
-# and lastely the third (if applicable)
+# 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?
-    physical_domain_size = (grid.limit_upper .- grid.limit_lower)
-    indices = Tuple.(CartesianIndices(grid.numberOfPointsPerDim))
-    return broadcast(I -> grid.limit_lower .+ physical_domain_size.*(I.-1), indices)
+    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<=numberOfDimensions(grid)
+    @assert dim<=dimension(grid)
     @assert dim>0
-    points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.numberOfPointsPerDim[dim])
+    points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.size[dim])
 end
-
-using PyPlot, PyCall
-
-function plotgridfunction(grid::EquidistantGrid, gridfunction)
-    if numberOfDimensions(grid) == 1
-        plot(pointsalongdim(grid,1), gridfunction, linewidth=2.0)
-    elseif numberOfDimensions(grid) == 2
-        mx = grid.numberOfPointsPerDim[1]
-        my = grid.numberOfPointsPerDim[2]
-        X = repeat(pointsalongdim(grid,1),1,my)
-        Y = permutedims(repeat(pointsalongdim(grid,2),1,mx))
-        plot_surface(X,Y,reshape(gridfunction,mx,my));
-    else
-        error(string("Plot not implemented for dimension ", string(numberOfDimensions(grid))))
-    end
-end

diff -r f01b70b81e95 -r 6b6d921e8f05 TODO.txt
--- a/TODO.txt	Thu Feb 07 16:34:55 2019 +0100
+++ b/TODO.txt	Thu Feb 21 14:19:25 2019 +0100
@@ -4,4 +4,9 @@
 Kolla att vi gör boundschecks överallt och att de är markerade med @boundscheck
 Kolla att vi har @inline på rätt ställen
 
-Ändra namn på variabler och funktioner så att det följer style-guide
\ No newline at end of file
+Ändra namn på variabler och funktioner så att det följer style-guide
+
+Profilera
+
+Konvertera till paket
+Skriv tester

diff -r f01b70b81e95 -r 6b6d921e8f05 diffOp.jl
--- a/diffOp.jl	Thu Feb 07 16:34:55 2019 +0100
+++ b/diffOp.jl	Thu Feb 21 14:19:25 2019 +0100
@@ -44,6 +44,53 @@
     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))
+    for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2)) # TBD: Is this the right way, the right size?
+        for j ∈ tileaxs[2], i ∈ tileaxs[1]
+            I = ri[i,j]
+            u[i,j] = 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)
@@ -62,23 +109,17 @@
 
 # u = L*v
 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
-    h = Grid.spacings(L.grid)[1]
-    uᵢ = L.a * apply(L.op, h, v, i)
+    uᵢ = L.a * apply(L.op, L.grid.spacing[1], v, i)
     return uᵢ
 end
 
-using UnsafeArrays
-
-function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
-    h = Grid.spacings(L.grid)
-
+@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 = uview(v, :, Int(I[2]))
-    @inbounds uᵢ  = apply(L.op, h[1], vx , I[1])
+    @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 = uview(v, Int(I[1]), :)
-    @inbounds uᵢ += apply(L.op, h[2], vy, I[2])
-
+    @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
 

diff -r f01b70b81e95 -r 6b6d921e8f05 index.jl
--- a/index.jl	Thu Feb 07 16:34:55 2019 +0100
+++ b/index.jl	Thu Feb 21 14:19:25 2019 +0100
@@ -38,3 +38,25 @@
 end
 
 IndexTuple(t::Vararg{Tuple{T, DataType}}) where T<:Integer = Index.(t)
+
+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

diff -r f01b70b81e95 -r 6b6d921e8f05 sbpD2.jl
--- a/sbpD2.jl	Thu Feb 07 16:34:55 2019 +0100
+++ b/sbpD2.jl	Thu Feb 21 14:19:25 2019 +0100
@@ -2,16 +2,16 @@
 
 # Apply for different regions Lower/Interior/Upper or Unknown region
 @inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Lower})
-    return @inbounds apply(op.closureStencils[Int(i)], v, Int(i))/h^2
+    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 apply(op.innerStencil, v, Int(i))/h^2
+    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 Int(op.parity)*apply(flip(op.closureStencils[N-Int(i)+1]), v, Int(i))/h^2
+    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})

diff -r f01b70b81e95 -r 6b6d921e8f05 sbpPlot.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/sbpPlot.jl	Thu Feb 21 14:19:25 2019 +0100
@@ -0,0 +1,17 @@
+module sbpPlot
+using PyPlot, PyCall
+
+function plotgridfunction(grid::EquidistantGrid, gridfunction)
+    if dimension(grid) == 1
+        plot(pointsalongdim(grid,1), gridfunction, linewidth=2.0)
+    elseif dimension(grid) == 2
+        mx = grid.size[1]
+        my = grid.size[2]
+        X = repeat(pointsalongdim(grid,1),1,my)
+        Y = permutedims(repeat(pointsalongdim(grid,2),1,mx))
+        plot_surface(X,Y,reshape(gridfunction,mx,my));
+    else
+        error(string("Plot not implemented for dimension ", string(dimension(grid))))
+    end
+end
+end

diff -r f01b70b81e95 -r 6b6d921e8f05 stencil.jl
--- a/stencil.jl	Thu Feb 07 16:34:55 2019 +0100
+++ b/stencil.jl	Thu Feb 21 14:19:25 2019 +0100
@@ -14,19 +14,25 @@
 end
 
 # Provides index into the Stencil based on offset for the root element
-function Base.getindex(s::Stencil, i::Int)
+@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 function apply(s::Stencil, v::AbstractVector, i::Int)
-    w = zero(eltype(v))
-    for j ∈ s.range[1]:s.range[2]
-        @inbounds weight = s[j]
-        w += weight*v[i+j]
+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