--- a/EquidistantGrid.jl	Thu Jan 24 17:46:57 2019 +0100
+++ b/EquidistantGrid.jl	Fri Jan 25 15:26:47 2019 +0100
@@ -77,29 +77,21 @@
 function pointsalongdim(grid::EquidistantGrid, dim::Integer)
     @assert dim<=numberOfDimensions(grid)
     @assert dim>0
-    points = range(grid.limits[1][dim],stop=grid.limits[2][dim],length=grid.numberOfPointsPerDim[dim])
+    points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.numberOfPointsPerDim[dim])
 end
 
 using PyPlot, PyCall
-#pygui(:qt)
-#using Plots; pyplot()
 
 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 = pointsalongdim(grid,1)
-        X = repeat(x,1,my)
-        y = pointsalongdim(grid,2)
-        Y = permutedims(repeat(y,1,mx))
-        plot_surface(X,Y,reshape(gridfunction,mx,my))
-        # fig = figure()
-        # ax = fig[:add_subplot](1,1,1, projection = "3d")
-        # ax[:plot_surface](X,Y,reshape(gridfunction,mx,my))
+        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
-    savefig("gridfunction")
 end

--- a/diffOp.jl	Thu Jan 24 17:46:57 2019 +0100
+++ b/diffOp.jl	Fri Jan 25 15:26:47 2019 +0100
@@ -47,10 +47,10 @@
     return u
 end
 
-struct Laplace{Dim,T<:Real} <: DiffOp
+struct Laplace{Dim,T<:Real,N,M,K} <: DiffOp
     grid::Grid.EquidistantGrid{Dim,T}
     a::T
-    op::D2{Float64}
+    op::D2{Float64,N,M,K}
 end
 
 # u = L*v

--- a/plotDerivative.jl	Thu Jan 24 17:46:57 2019 +0100
+++ b/plotDerivative.jl	Fri Jan 25 15:26:47 2019 +0100
@@ -1,6 +1,6 @@
-g = sbp.Grid.EquidistantGrid(200,(0, 2pi))
+g = sbp.Grid.EquidistantGrid((200,), (0.0,), (2pi,))
 op =sbp.readOperator("d2_4th.txt","h_4th.txt")
-Laplace = sbp.Laplace1D(g,1,op)
+Laplace = sbp.Laplace(g,1.0,op)
 
 init(x) = cos(x)
 v = sbp.Grid.evalOn(g,init)

--- a/plotDerivative2d.jl	Thu Jan 24 17:46:57 2019 +0100
+++ b/plotDerivative2d.jl	Fri Jan 25 15:26:47 2019 +0100
@@ -1,6 +1,6 @@
-g = sbp.Grid.EquidistantGrid((100,75),((0, 0), (2pi, 3/2*pi)))
+g = sbp.Grid.EquidistantGrid((100,75), (0.0, 0.0), (2pi, 3/2*pi))
 op = sbp.readOperator("d2_4th.txt","h_4th.txt")
-Laplace = sbp.Laplace2D(g,1,op)
+Laplace = sbp.Laplace(g, 1.0, op)
 
 init(x,y) = sin(x) + cos(y)
 v = sbp.Grid.evalOn(g,init)

--- a/sbpD2.jl	Thu Jan 24 17:46:57 2019 +0100
+++ b/sbpD2.jl	Fri Jan 25 15:26:47 2019 +0100
@@ -20,10 +20,10 @@
     even = 1
 end
 
-struct D2{T} <: ConstantStencilOperator
+struct D2{T,N,M,K} <: ConstantStencilOperator
     quadratureClosure::Vector{T}
-    innerStencil::Stencil{T}
-    closureStencils::Vector{Stencil{T}} # TBD: Should this be a tuple?
+    innerStencil::Stencil{T,N}
+    closureStencils::NTuple{M, Stencil{T,K}}
     eClosure::Vector{T}
     dClosure::Vector{T}
     parity::Parity
@@ -42,7 +42,7 @@
     width = length(innerStencilWeights)
     r = (-div(width,2), div(width,2))
 
-    innerStencil = Stencil(r, innerStencilWeights)
+    innerStencil = Stencil(r, Tuple(innerStencilWeights))
 
     # Create boundary stencils
     boundarySize = length(d["boundary_stencils"])
@@ -52,7 +52,7 @@
         stencilWeights = stringToVector(Float64, d["boundary_stencils"][i])
         width = length(stencilWeights)
         r = (1-i,width-i)
-        push!(closureStencils,Stencil(r, stencilWeights))
+        closureStencils = (closureStencils..., Stencil(r, Tuple(stencilWeights)))
     end
 
     d2 = D2(

--- a/stencil.jl	Thu Jan 24 17:46:57 2019 +0100
+++ b/stencil.jl	Fri Jan 25 15:26:47 2019 +0100
@@ -1,19 +1,11 @@
-struct Stencil{T<:Real}
-    range::NTuple{2,Int}
-    weights::Vector{T} # Should this be a tuple?? (Check type stability)
-
-    function Stencil(range, weights)
-        width = range[2]-range[1]+1
-        if width != length(weights)
-            error("The width and the number of weights must be the same")
-        end
-        new{eltype(weights)}(range, weights)
-    end
+struct Stencil{T<:Real,N}
+    range::Tuple{Int,Int}
+    weights::NTuple{N,T}
 end
 
 function flip(s::Stencil)
     range = (-s.range[2], -s.range[1])
-    s = Stencil(range, s.weights[end:-1:1])
+    return Stencil(range, reverse(s.weights))
 end
 
 # Provides index into the Stencil based on offset for the root element