--- a/AbstractGrid.jl	Fri Jan 25 15:20:40 2019 +0100
+++ b/AbstractGrid.jl	Thu Feb 21 16:27:28 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
 

--- a/EquidistantGrid.jl	Fri Jan 25 15:20:40 2019 +0100
+++ b/EquidistantGrid.jl	Thu Feb 21 16:27:28 2019 +0100
@@ -5,93 +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
 
-# 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)
-    dx̄ = (grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1)
-
-    points = Vector{typeof(dx̄)}(undef, numberOfPoints(grid))
-    # Compute the points based on their Cartesian indices and the signed
-    # grid spacings
-    cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim)
-    for i ∈ 1:numberOfPoints(grid)
-        ci = Tuple(cartesianIndices[i]) .-1
-        points[i] = grid.limit_lower .+ dx̄.*ci
-    end
-
-    # TBD: Keep? this? How do we want to represent points in 1D?
-    if numberOfDimensions(grid) == 1
-        points = broadcast(x -> x[1], points)
-    end
-    return points
+    # 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<=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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/TODO.txt	Thu Feb 21 16:27:28 2019 +0100
@@ -0,0 +1,12 @@
+# TODO
+
+Kolla att vi har @inbounds och @propagate_inbounds på rätt ställen
+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
+
+Profilera
+
+Konvertera till paket
+Skriv tester

--- a/diffOp.jl	Fri Jan 25 15:20:40 2019 +0100
+++ b/diffOp.jl	Thu Feb 21 16:27:28 2019 +0100
@@ -1,6 +1,7 @@
 abstract type DiffOp end
 
-function apply!(D::DiffOp, u::AbstractVector, v::AbstractVector)
+# 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
 
@@ -32,50 +33,98 @@
     error("not implemented")
 end
 
-# Differential operator for a*d^2/dx^2
-struct Laplace{Dim,T<:Real,N,M,K} <: DiffOp
+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))
+    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)
+    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
 
-# u = L*v
-function apply!(L::Laplace{1}, u::AbstractVector, v::AbstractVector)
-    h = Grid.spacings(L.grid)[1]
-    apply!(L.op, u, v, h)
-    u .= L.a * u
-    return nothing
+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{2}, u::AbstractVector, v::AbstractVector)
-    u .= 0*u
-    h = Grid.spacings(L.grid)
-
-    li = LinearIndices(L.grid.numberOfPointsPerDim)
-    n_x, n_y = L.grid.numberOfPointsPerDim
-
-
-    # For each x
-    temp = zeros(eltype(u), n_y)
-    for i ∈ 1:n_x
-
-        v_i = view(v, li[i,:])
-        apply!(L.op, temp, v_i, h[2])
+function apply(L::Laplace{1}, v::AbstractVector, i::Int)
+    uᵢ = L.a * apply(L.op, L.grid.spacing[1], v, i)
+    return uᵢ
+end
 
-        u[li[i,:]] += temp
-    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
 
-    # For each y
-    temp = zeros(eltype(u), n_x)
-    for i ∈ 1:n_y
-        v_i = view(v, li[:,i])
-        apply!(L.op, temp, v_i, h[1])
-
-        u[li[:,i]] += temp
-    end
-
-    u .= L.a*u
-
-    return nothing
+# 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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/index.jl	Thu Feb 21 16:27:28 2019 +0100
@@ -0,0 +1,62 @@
+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)
+
+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/plotDerivative.jl	Fri Jan 25 15:20:40 2019 +0100
+++ b/plotDerivative.jl	Thu Feb 21 16:27:28 2019 +0100
@@ -2,7 +2,7 @@
 op =sbp.readOperator("d2_4th.txt","h_4th.txt")
 Laplace = sbp.Laplace(g,1.0,op)
 
-init(x) = sin(x)
+init(x) = cos(x)
 v = sbp.Grid.evalOn(g,init)
 u = zeros(length(v))
 

--- a/plotDerivative2d.jl	Fri Jan 25 15:20:40 2019 +0100
+++ b/plotDerivative2d.jl	Thu Feb 21 16:27:28 2019 +0100
@@ -9,8 +9,8 @@
 
 sbp.apply!(Laplace,u,v)
 
-@show u
-@show u'*u
+#@show u
+#@show u'*u
 
 sbp.Grid.plotgridfunction(g,u)
 

--- a/sbp.jl	Fri Jan 25 15:20:40 2019 +0100
+++ b/sbp.jl	Thu Feb 21 16:27:28 2019 +0100
@@ -1,4 +1,5 @@
 module sbp
+include("index.jl")
 include("grid.jl")
 include("stencil.jl")
 include("sbpD2.jl")

--- a/sbpD2.jl	Fri Jan 25 15:20:40 2019 +0100
+++ b/sbpD2.jl	Thu Feb 21 16:27:28 2019 +0100
@@ -1,24 +1,40 @@
 abstract type ConstantStencilOperator end
 
-function apply!(op::ConstantStencilOperator, u::AbstractVector, v::AbstractVector, h::Real)
-    N = length(v)
-    cSize = closureSize(op)
+# 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
 
-    for i ∈ range(1; length=cSize)
-        u[i] = apply(op.closureStencils[i], v, i)/h^2
-    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)
 
-    innerStart = 1 + cSize
-    innerEnd = N - cSize
-    for i ∈ range(innerStart, stop=innerEnd)
-        u[i] = apply(op.innerStencil, v, i)/h^2
-    end
+    i = Int(index)
 
-    for i ∈ range(innerEnd+1, length=cSize)
-        u[i] = Int(op.parity)*apply(flip(op.closureStencils[N-i+1]), v, i)/h^2
+    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
 
-    return nothing
+
+# 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
@@ -27,11 +43,11 @@
 end
 
 struct D2{T,N,M,K} <: ConstantStencilOperator
-    quadratureClosure::Vector{T}
+    quadratureClosure::NTuple{M,T}
     innerStencil::Stencil{T,N}
-    closureStencils::NTuple{M, Stencil{T,K}}
-    eClosure::Vector{T}
-    dClosure::Vector{T}
+    closureStencils::NTuple{M,Stencil{T,K}}
+    eClosure::Stencil{T,M}
+    dClosure::Stencil{T,M}
     parity::Parity
 end
 
@@ -44,29 +60,33 @@
     h = readSectionedFile(Hfn)
 
     # Create inner stencil
-    innerStencilWeights = stringToVector(Float64, d["inner_stencil"][1])
+    innerStencilWeights = stringToTuple(Float64, d["inner_stencil"][1])
     width = length(innerStencilWeights)
     r = (-div(width,2), div(width,2))
 
-    innerStencil = Stencil(r, Tuple(innerStencilWeights))
+    innerStencil = Stencil(r, innerStencilWeights)
 
     # Create boundary stencils
     boundarySize = length(d["boundary_stencils"])
-    closureStencils = Vector{Stencil}()
+    closureStencils = Vector{typeof(innerStencil)}() # TBD: is the the right way to get the correct type?
 
     for i ∈ 1:boundarySize
-        stencilWeights = stringToVector(Float64, d["boundary_stencils"][i])
+        stencilWeights = stringToTuple(Float64, d["boundary_stencils"][i])
         width = length(stencilWeights)
         r = (1-i,width-i)
-        closureStencils = (closureStencils..., Stencil(r, Tuple(stencilWeights)))
+        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(
-        stringToVector(Float64, h["closure"][1]),
+        quadratureClosure,
         innerStencil,
         closureStencils,
-        stringToVector(Float64, d["e"][1]),
-        stringToVector(Float64, d["d1"][1]),
+        eClosure,
+        dClosure,
         even
     )
 
@@ -74,6 +94,23 @@
 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}}()
@@ -98,6 +135,19 @@
     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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/sbpPlot.jl	Thu Feb 21 16:27:28 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

--- a/stencil.jl	Fri Jan 25 15:20:40 2019 +0100
+++ b/stencil.jl	Thu Feb 21 16:27:28 2019 +0100
@@ -1,6 +1,11 @@
 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)
@@ -9,19 +14,25 @@
 end
 
 # Provides index into the Stencil based on offset for the root element
-function Base.getindex(s::Stencil, i::Int)
-    # TBD: Rearrange to mark with @boundscheck?
-    if s.range[1] <= i <= s.range[2]
-        return s.weights[1 + i - s.range[1]]
-    else
-        return 0
+@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
 
-function apply(s::Stencil, v::AbstractVector, i::Int)
-    w = zero(eltype(v))
-    for j ∈ s.range[1]:s.range[2]
-        w += s[j]*v[i+j] # TBD: Make this without boundschecks?
+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