diff -r d12ab8120d29 -r 101b88d13cf3 Manifest.toml
--- a/Manifest.toml	Wed Mar 23 12:43:03 2022 +0100
+++ b/Manifest.toml	Thu Apr 07 20:28:08 2022 +0200
@@ -1,13 +1,13 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.7.0"
+julia_version = "1.7.1"
 manifest_format = "2.0"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
+git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "3.3.1"
+version = "3.3.3"
 
 [[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

diff -r d12ab8120d29 -r 101b88d13cf3 Notes.md
--- a/Notes.md	Wed Mar 23 12:43:03 2022 +0100
+++ b/Notes.md	Thu Apr 07 20:28:08 2022 +0200
@@ -146,12 +146,52 @@
  - [ ] Is it ok to have "Constructors" for abstract types which create subtypes? For example a Grids() functions that gives different kind of grids based on input?
  - [ ] Figure out how to treat the borrowing parameters of operators. Include in into the struct? Expose via function dispatched on the operator type and grid?
 
+## Identifiers for regions
+The identifiers (`Upper`, `Lower`, `Interior`) used for region indecies should probabily be included in the grid module. This allows new grid types to come with their own regions.
+
 ## Regions and tensormappings
 - [ ] Use a trait to indicate if a LazyTensor uses indices with regions.
     The default should be that they do NOT.
         - [ ] What to name this trait? Can we call it IndexStyle but not export it to avoid conflicts with Base.IndexStyle?
  - [ ] Figure out repeated application of regioned LazyTensors. Maybe an instance of a tensor mapping needs to know the exact size of the range and domain for this to work?
 
+### Ideas for information sharing functions
+```julia
+using StaticArrays
+
+function regions(op::SecondDerivativeVariable)
+    t = ntuple(i->(Interior(),),range_dim(op))
+    return Base.setindex(t, (Lower(), Interior(), Upper()), derivative_direction(op))
+end
+
+function regionsizes(op::SecondDerivativeVariable)
+    sz = tuple.(range_size(op))
+
+    cl = closuresize(op)
+    return Base.setindex(sz, (cl, n-2cl, cl), derivative_direction(op))
+end
+
+
+g = EquidistantGrid((11,9), (0.,0.), (10.,8.)) # h = 1
+c = evalOn(g, (x,y)->x+y)
+
+D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,1)
+@test regions(D₂ᶜ) == (
+    (Lower(), Interior(), Upper()),
+    (Interior(),),
+)
+@test regionsizes(D₂ᶜ) == ((1,9,1),(9,))
+
+
+D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,2)
+@test regions(D₂ᶜ) == (
+    (Interior(),),
+    (Lower(), Interior(), Upper()),
+)
+@test regionsizes(D₂ᶜ) == ((11,),(1,7,1))
+```
+
+
 ## Boundschecking and dimension checking
 Does it make sense to have boundschecking only in getindex methods?
 This would mean no bounds checking in applys, however any indexing that they do would be boundschecked. The only loss would be readability of errors. But users aren't really supposed to call apply directly anyway.

diff -r d12ab8120d29 -r 101b88d13cf3 TODO.md
--- a/TODO.md	Wed Mar 23 12:43:03 2022 +0100
+++ b/TODO.md	Thu Apr 07 20:28:08 2022 +0200
@@ -1,14 +1,13 @@
 # TODO
 
-## Skämskudde
- - [ ] Ändra namn på variabler och funktioner så att det följer style-guide
- - [ ] Skriv tester
 
 ## Coding
+ - [ ] Ändra namn på variabler och funktioner så att det följer style-guide
  - [ ] Add new Laplace operator to DiffOps, probably named WaveEqOp(?!!?)
  - [ ] Create a struct that bundles the necessary Tensor operators for solving the wave equation.
  - [ ] Replace getindex hack for flattening tuples with flatten_tuple. (eg. `getindex.(range_size.(L.D2),1)`)
  - [ ] Use `@inferred` in a lot of tests.
+ - [ ] Replace `@inferred` tests with a benchmark suite that automatically tests for regressions.
  - [ ] Make sure we are setting tolerances in tests in a consistent way
  - [ ] Write down some coding guideline or checklist for code conventions. For example i,j,... for indices and I for multi-index
  - [ ] Clean up RegionIndices
@@ -19,12 +18,11 @@
  - [ ] Add possibility to create tensor mapping application with `()`, e.g `D1(v) <=> D1*v`?
  - [ ] Add custom pretty printing to LazyTensors/SbpOperators to enhance readability of e.g error messages.
        See (https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing)
+ - [ ] Samla noggrannhets- och SBP-ness-tester för alla operatorer på ett ställe
  - [ ] Move export statements to top of each module
 
-## Repo
- - [ ] Rename repo to Sbplib.jl
 
-# Wrap up tasks
+ - [ ] Gå igenom alla typ parametrar och kolla om de är motiverade. Både i signaturer och typer, tex D i VariableSecondDerivative. Kan vi använda promote istället?
  - [ ] 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

diff -r d12ab8120d29 -r 101b88d13cf3 src/SbpOperators/stencil.jl
--- a/src/SbpOperators/stencil.jl	Wed Mar 23 12:43:03 2022 +0100
+++ b/src/SbpOperators/stencil.jl	Thu Apr 07 20:28:08 2022 +0200
@@ -1,11 +1,12 @@
 export CenteredStencil
+export CenteredNestedStencil
 
 struct Stencil{T,N}
-    range::Tuple{Int,Int}
+    range::UnitRange{Int64}
     weights::NTuple{N,T}
 
-    function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T, N}
-        @assert range[2]-range[1]+1 == N
+    function Stencil(range::UnitRange,weights::NTuple{N,T}) where {T, N}
+        @assert length(range) == N
         new{T,N}(range,weights)
     end
 end
@@ -15,27 +16,30 @@
 
 Create a stencil with the given weights with element `center` as the center of the stencil.
 """
-function Stencil(weights::Vararg{T}; center::Int) where T # Type parameter T makes sure the weights are valid for the Stencil constuctors and throws an earlier, more readable, error
+function Stencil(weights...; center::Int)
+    weights = promote(weights...)
     N = length(weights)
-    range = (1, N) .- center
+    range = (1:N) .- center
 
     return Stencil(range, weights)
 end
 
-function Stencil{T}(s::Stencil) where T
-    return Stencil(s.range, T.(s.weights))
-end
+Stencil{T,N}(s::Stencil{S,N}) where {T,S,N} = Stencil(s.range, T.(s.weights))
+Stencil{T}(s::Stencil) where T = Stencil{T,length(s)}(s)
 
-Base.convert(::Type{Stencil{T}}, stencil) where T = Stencil{T}(stencil)
+Base.convert(::Type{Stencil{T1,N}}, s::Stencil{T2,N}) where {T1,T2,N} = Stencil{T1,N}(s)
+Base.convert(::Type{Stencil{T1}}, s::Stencil{T2,N}) where {T1,T2,N} = Stencil{T1,N}(s)
 
-function CenteredStencil(weights::Vararg)
+Base.promote_rule(::Type{Stencil{T1,N}}, ::Type{Stencil{T2,N}}) where {T1,T2,N} = Stencil{promote_type(T1,T2),N}
+
+function CenteredStencil(weights...)
     if iseven(length(weights))
         throw(ArgumentError("a centered stencil must have an odd number of weights."))
     end
 
     r = length(weights) ÷ 2
 
-    return Stencil((-r, r), weights)
+    return Stencil(-r:r, weights)
 end
 
 
@@ -48,7 +52,8 @@
     return Stencil(s.range, a.*s.weights)
 end
 
-Base.eltype(::Stencil{T}) where T = T
+Base.eltype(::Stencil{T,N}) where {T,N} = T
+Base.length(::Stencil{T,N}) where {T,N} = N
 
 function flip(s::Stencil)
     range = (-s.range[2], -s.range[1])
@@ -57,24 +62,103 @@
 
 # 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
+    @boundscheck if i ∉ s.range
         return zero(eltype(s))
     end
     return s.weights[1 + i - s.range[1]]
 end
 
-Base.@propagate_inbounds @inline function apply_stencil(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
-    w = s.weights[1]*v[i + s.range[1]]
-    @simd for k ∈ 2:N
-        w += s.weights[k]*v[i + s.range[1] + k-1]
+Base.@propagate_inbounds @inline function apply_stencil(s::Stencil, v::AbstractVector, i::Int)
+    w = zero(promote_type(eltype(s),eltype(v)))
+    @simd for k ∈ 1:length(s)
+        w += s.weights[k]*v[i + s.range[k]]
+    end
+
+    return w
+end
+
+Base.@propagate_inbounds @inline function apply_stencil_backwards(s::Stencil, v::AbstractVector, i::Int)
+    w = zero(promote_type(eltype(s),eltype(v)))
+    @simd for k ∈ length(s):-1:1
+        w += s.weights[k]*v[i - s.range[k]]
     end
     return w
 end
 
-Base.@propagate_inbounds @inline function apply_stencil_backwards(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
-    w = s.weights[N]*v[i - s.range[2]]
-    @simd for k ∈ N-1:-1:1
-        w += s.weights[k]*v[i - s.range[1] - k + 1]
-    end
-    return w
+
+struct NestedStencil{T,N,M}
+    s::Stencil{Stencil{T,N},M}
+end
+
+# Stencil input
+NestedStencil(s::Vararg{Stencil}; center) = NestedStencil(Stencil(s... ; center))
+CenteredNestedStencil(s::Vararg{Stencil}) = NestedStencil(CenteredStencil(s...))
+
+# Tuple input
+function NestedStencil(weights::Vararg{NTuple{N,Any}}; center) where N
+    inner_stencils = map(w -> Stencil(w...; center), weights)
+    return NestedStencil(Stencil(inner_stencils... ; center))
+end
+function CenteredNestedStencil(weights::Vararg{NTuple{N,Any}}) where N
+    inner_stencils = map(w->CenteredStencil(w...), weights)
+    return CenteredNestedStencil(inner_stencils...)
+end
+
+
+# Conversion
+function NestedStencil{T,N,M}(ns::NestedStencil{S,N,M}) where {T,S,N,M}
+    return NestedStencil(Stencil{Stencil{T}}(ns.s))
+end
+
+function NestedStencil{T}(ns::NestedStencil{S,N,M}) where {T,S,N,M}
+    NestedStencil{T,N,M}(ns)
+end
+
+function Base.convert(::Type{NestedStencil{T,N,M}}, s::NestedStencil{S,N,M}) where {T,S,N,M}
+    return NestedStencil{T,N,M}(s)
+end
+Base.convert(::Type{NestedStencil{T}}, stencil) where T = NestedStencil{T}(stencil)
+
+function Base.promote_rule(::Type{NestedStencil{T,N,M}}, ::Type{NestedStencil{S,N,M}}) where {T,S,N,M}
+    return NestedStencil{promote_type(T,S),N,M}
 end
+
+Base.eltype(::NestedStencil{T}) where T = T
+
+function scale(ns::NestedStencil, a)
+    range = ns.s.range
+    weights = ns.s.weights
+
+    return NestedStencil(Stencil(range, scale.(weights,a)))
+end
+
+function flip(ns::NestedStencil)
+    s_flip = flip(ns.s)
+    return NestedStencil(Stencil(s_flip.range, flip.(s_flip.weights)))
+end
+
+Base.getindex(ns::NestedStencil, i::Int) = ns.s[i]
+
+"Apply inner stencils to `c` and get a concrete stencil"
+Base.@propagate_inbounds function apply_inner_stencils(ns::NestedStencil, c::AbstractVector, i::Int)
+    weights = apply_stencil.(ns.s.weights, Ref(c), i)
+    return Stencil(ns.s.range, weights)
+end
+
+"Apply the whole nested stencil"
+Base.@propagate_inbounds function apply_stencil(ns::NestedStencil, c::AbstractVector, v::AbstractVector, i::Int)
+    s = apply_inner_stencils(ns,c,i)
+    return apply_stencil(s, v, i)
+end
+
+"Apply inner stencils backwards to `c` and get a concrete stencil"
+Base.@propagate_inbounds @inline function apply_inner_stencils_backwards(ns::NestedStencil, c::AbstractVector, i::Int)
+    weights = apply_stencil_backwards.(ns.s.weights, Ref(c), i)
+    return Stencil(ns.s.range, weights)
+end
+
+"Apply the whole nested stencil backwards"
+Base.@propagate_inbounds @inline function apply_stencil_backwards(ns::NestedStencil, c::AbstractVector, v::AbstractVector, i::Int)
+    s = apply_inner_stencils_backwards(ns,c,i)
+    return apply_stencil_backwards(s, v, i)
+end

diff -r d12ab8120d29 -r 101b88d13cf3 src/SbpOperators/volumeops/volume_operator.jl
--- a/src/SbpOperators/volumeops/volume_operator.jl	Wed Mar 23 12:43:03 2022 +0100
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Thu Apr 07 20:28:08 2022 +0200
@@ -24,7 +24,7 @@
 end
 
 """
-    VolumeOperator{T,N,M,K} <: TensorOperator{T,1}
+    VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1}
 Implements a one-dimensional constant coefficients volume operator
 """
 struct VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1}

diff -r d12ab8120d29 -r 101b88d13cf3 test/Manifest.toml
--- a/test/Manifest.toml	Wed Mar 23 12:43:03 2022 +0100
+++ b/test/Manifest.toml	Thu Apr 07 20:28:08 2022 +0200
@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.7.0"
+julia_version = "1.7.1"
 manifest_format = "2.0"
 
 [[deps.ArgTools]]
@@ -14,9 +14,9 @@
 
 [[deps.ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "4c26b4e9e91ca528ea212927326ece5918a04b47"
+git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
 uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.11.2"
+version = "1.14.0"
 
 [[deps.ChangesOfVariables]]
 deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
@@ -26,9 +26,9 @@
 
 [[deps.Compat]]
 deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
-git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582"
+git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "3.41.0"
+version = "3.42.0"
 
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -48,10 +48,10 @@
 uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 
 [[deps.DiffRules]]
-deps = ["LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
-git-tree-sha1 = "9bc5dac3c8b6706b58ad5ce24cffd9861f07c94f"
+deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
+git-tree-sha1 = "dd933c4ef7b4c270aacd4eb88fa64c147492acf0"
 uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
-version = "1.9.0"
+version = "1.10.0"
 
 [[deps.Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
@@ -78,9 +78,9 @@
 
 [[deps.InverseFunctions]]
 deps = ["Test"]
-git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
+git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
 uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
-version = "0.1.2"
+version = "0.1.3"
 
 [[deps.IrrationalConstants]]
 git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
@@ -89,9 +89,9 @@
 
 [[deps.JLLWrappers]]
 deps = ["Preferences"]
-git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
+git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.3.0"
+version = "1.4.1"
 
 [[deps.LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
@@ -118,9 +118,9 @@
 
 [[deps.LogExpFunctions]]
 deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1"
+git-tree-sha1 = "58f25e56b706f95125dcb796f39e1fb01d913a71"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.6"
+version = "0.3.10"
 
 [[deps.Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
@@ -140,9 +140,9 @@
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
 
 [[deps.NaNMath]]
-git-tree-sha1 = "f755f36b19a5116bb580de457cda0c140153f283"
+git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5"
 uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-version = "0.3.6"
+version = "1.0.0"
 
 [[deps.NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
@@ -167,9 +167,9 @@
 
 [[deps.Preferences]]
 deps = ["TOML"]
-git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
+git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
 uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.2.2"
+version = "1.2.5"
 
 [[deps.Printf]]
 deps = ["Unicode"]
@@ -185,9 +185,9 @@
 
 [[deps.Requires]]
 deps = ["UUIDs"]
-git-tree-sha1 = "8f82019e525f4d5c669692772a6f4b0a58b06a6a"
+git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
-version = "1.2.0"
+version = "1.3.0"
 
 [[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
@@ -208,9 +208,9 @@
 
 [[deps.SpecialFunctions]]
 deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "e08890d19787ec25029113e88c34ec20cac1c91e"
+git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.0.0"
+version = "2.1.4"
 
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
@@ -236,9 +236,9 @@
 
 [[deps.Tullio]]
 deps = ["ChainRulesCore", "DiffRules", "LinearAlgebra", "Requires"]
-git-tree-sha1 = "0288b7a395fc412952baf756fac94e4f28bfec65"
+git-tree-sha1 = "7830c974acc69437a3fee35dd7b510a74cbc862d"
 uuid = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
-version = "0.3.2"
+version = "0.3.3"
 
 [[deps.UUIDs]]
 deps = ["Random", "SHA"]

diff -r d12ab8120d29 -r 101b88d13cf3 test/SbpOperators/boundaryops/boundary_operator_test.jl
--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl	Wed Mar 23 12:43:03 2022 +0100
+++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl	Thu Apr 07 20:28:08 2022 +0200
@@ -9,7 +9,7 @@
 import Sbplib.SbpOperators.boundary_operator
 
 @testset "BoundaryOperator" begin
-    closure_stencil = Stencil((0,2), (2.,1.,3.))
+    closure_stencil = Stencil(2.,1.,3.; center = 1)
     g_1D = EquidistantGrid(11, 0.0, 1.0)
     g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
 

diff -r d12ab8120d29 -r 101b88d13cf3 test/SbpOperators/stencil_test.jl
--- a/test/SbpOperators/stencil_test.jl	Wed Mar 23 12:43:03 2022 +0100
+++ b/test/SbpOperators/stencil_test.jl	Thu Apr 07 20:28:08 2022 +0200
@@ -1,19 +1,25 @@
 using Test
 using Sbplib.SbpOperators
 import Sbplib.SbpOperators.Stencil
+import Sbplib.SbpOperators.NestedStencil
+import Sbplib.SbpOperators.scale
 
 @testset "Stencil" begin
-    s = Stencil((-2,2), (1.,2.,2.,3.,4.))
+    s = Stencil(-2:2, (1.,2.,2.,3.,4.))
     @test s isa Stencil{Float64, 5}
 
     @test eltype(s) == Float64
-    @test SbpOperators.scale(s, 2) == Stencil((-2,2), (2.,4.,4.,6.,8.))
+
+    @test length(s) == 5
+    @test length(Stencil(-1:2, (1,2,3,4))) == 4
+
+    @test SbpOperators.scale(s, 2) == Stencil(-2:2, (2.,4.,4.,6.,8.))
 
-    @test Stencil(1,2,3,4; center=1) == Stencil((0, 3),(1,2,3,4))
-    @test Stencil(1,2,3,4; center=2) == Stencil((-1, 2),(1,2,3,4))
-    @test Stencil(1,2,3,4; center=4) == Stencil((-3, 0),(1,2,3,4))
+    @test Stencil(1,2,3,4; center=1) == Stencil(0:3,(1,2,3,4))
+    @test Stencil(1,2,3,4; center=2) == Stencil(-1:2,(1,2,3,4))
+    @test Stencil(1,2,3,4; center=4) == Stencil(-3:0,(1,2,3,4))
 
-    @test CenteredStencil(1,2,3,4,5) == Stencil((-2, 2), (1,2,3,4,5))
+    @test CenteredStencil(1,2,3,4,5) == Stencil(-2:2, (1,2,3,4,5))
     @test_throws ArgumentError CenteredStencil(1,2,3,4)
 
     # Changing the type of the weights
@@ -24,8 +30,145 @@
 
     @testset "convert" begin
         @test convert(Stencil{Float64}, Stencil(1,2,3,4,5; center=2)) == Stencil(1.,2.,3.,4.,5.; center=2)
-        @test convert(Stencil{Float64}, CenteredStencil(1,2,3,4,5)) == CenteredStencil(1.,2.,3.,4.,5.)
-        @test convert(Stencil{Int}, Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1,2,3,4,5; center=2)
-        @test convert(Stencil{Rational}, Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1//1,2//1,3//1,4//1,5//1; center=2)
+        @test convert(Stencil{Float64,5}, CenteredStencil(1,2,3,4,5)) == CenteredStencil(1.,2.,3.,4.,5.)
+        @test convert(Stencil{Int,5}, Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1,2,3,4,5; center=2)
+        @test convert(Stencil{Rational,5}, Stencil(1.,2.,3.,4.,5.; center=2)) == Stencil(1//1,2//1,3//1,4//1,5//1; center=2)
+    end
+
+    @testset "promotion of weights" begin
+        @test Stencil(1.,2; center = 1) isa Stencil{Float64, 2}
+        @test Stencil(1,2//2; center = 1) isa Stencil{Rational{Int64}, 2}
+    end
+
+    @testset "promotion" begin
+        @test promote(Stencil(1,1;center=1), Stencil(2.,2.;center=2)) == (Stencil(1.,1.;center=1), Stencil(2.,2.;center=2))
+    end
+
+    @testset "type stability" begin
+        s_int = CenteredStencil(1,2,3)
+        s_float = CenteredStencil(1.,2.,3.)
+        v_int = rand(1:10,10);
+        v_float = rand(10);
+
+        @inferred SbpOperators.apply_stencil(s_int, v_int, 2)
+        @inferred SbpOperators.apply_stencil(s_float, v_float, 2)
+        @inferred SbpOperators.apply_stencil(s_int,  v_float, 2)
+        @inferred SbpOperators.apply_stencil(s_float, v_int, 2)
+
+        @inferred SbpOperators.apply_stencil_backwards(s_int, v_int, 5)
+        @inferred SbpOperators.apply_stencil_backwards(s_float, v_float, 5)
+        @inferred SbpOperators.apply_stencil_backwards(s_int,  v_float, 5)
+        @inferred SbpOperators.apply_stencil_backwards(s_float, v_int, 5)
     end
 end
+
+@testset "NestedStencil" begin
+
+    @testset "Constructors" begin
+        s1 = CenteredStencil(-1, 1, 0)
+        s2 = CenteredStencil(-1, 0, 1)
+        s3 = CenteredStencil( 0,-1, 1)
+
+        ns = NestedStencil(CenteredStencil(s1,s2,s3))
+        @test ns isa NestedStencil{Int,3}
+
+        @test CenteredNestedStencil(s1,s2,s3) == ns
+
+        @test NestedStencil(s1,s2,s3, center = 2) == ns
+        @test NestedStencil(s1,s2,s3, center = 1) == NestedStencil(Stencil(s1,s2,s3, center=1))
+
+        @test NestedStencil((-1,1,0),(-1,0,1),(0,-1,1), center=2) == ns
+        @test CenteredNestedStencil((-1,1,0),(-1,0,1),(0,-1,1)) == ns
+        @test NestedStencil((-1,1,0),(-1,0,1),(0,-1,1), center=1) == NestedStencil(Stencil(
+            Stencil(-1, 1, 0; center=1),
+            Stencil(-1, 0, 1; center=1),
+            Stencil( 0,-1, 1; center=1);
+            center=1
+        ))
+
+        @testset "Error handling" begin
+        end
+    end
+
+    @testset "scale" begin
+        ns = NestedStencil((-1,1,0),(-1,0,1),(0,-1,1), center=2)
+        @test SbpOperators.scale(ns, 2) == NestedStencil((-2,2,0),(-2,0,2),(0,-2,2), center=2)
+    end
+
+    @testset "conversion" begin
+        ns = NestedStencil((-1,1,0),(-1,0,1),(0,-1,1), center=2)
+        @test NestedStencil{Float64}(ns) == NestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-1.,1.), center=2)
+        @test NestedStencil{Rational}(ns) == NestedStencil((-1//1,1//1,0//1),(-1//1,0//1,1//1),(0//1,-1//1,1//1), center=2)
+
+        @test convert(NestedStencil{Float64}, ns) == NestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-1.,1.), center=2)
+        @test convert(NestedStencil{Rational}, ns) == NestedStencil((-1//1,1//1,0//1),(-1//1,0//1,1//1),(0//1,-1//1,1//1), center=2)
+    end
+
+    @testset "promotion of weights" begin
+        @test NestedStencil((-1,1,0),(-1.,0.,1.),(0,-1,1), center=2) isa NestedStencil{Float64,3,3}
+        @test NestedStencil((-1,1,0),(-1,0,1),(0//1,-1,1), center=2) isa NestedStencil{Rational{Int64},3,3}
+    end
+
+    @testset "promotion" begin
+        promote(
+            CenteredNestedStencil((-1,1,0),(-1,0,1),(0,-1,1)),
+            CenteredNestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-1.,1.))
+        ) == (
+            CenteredNestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-1.,1.)),
+            CenteredNestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-1.,1.))
+        )
+    end
+
+    @testset "apply" begin
+        c = [  1,  3,  6, 10, 15, 21, 28, 36, 45, 55]
+        v = [  2,  3,  5,  7, 11, 13, 17, 19, 23, 29]
+
+        # Centered
+        ns = NestedStencil((-1,1,0),(-1,0,1),(0,-2,2), center=2)
+        @test SbpOperators.apply_inner_stencils(ns, c, 4) == Stencil(4,9,10; center=2)
+        @test SbpOperators.apply_inner_stencils_backwards(ns, c, 4) == Stencil(-5,-9,-8; center=2)
+
+        @test SbpOperators.apply_stencil(ns, c, v, 4) == 4*5 + 9*7 + 10*11
+        @test SbpOperators.apply_stencil_backwards(ns, c, v, 4) == -8*5 - 9*7 - 5*11
+
+        # Non-centered
+        ns = NestedStencil((-1,1,0),(-1,0,1),(0,-1,1), center=1)
+        @test SbpOperators.apply_inner_stencils(ns, c, 4) == Stencil(5,11,6; center=1)
+        @test SbpOperators.apply_inner_stencils_backwards(ns, c, 4) == Stencil(-4,-7,-3; center=1)
+
+        @test SbpOperators.apply_stencil(ns, c, v, 4) == 5*7 + 11*11 + 6*13
+        @test SbpOperators.apply_stencil_backwards(ns, c, v, 4) == -3*3 - 7*5 - 4*7
+    end
+
+    @testset "type stability" begin
+        s_int = CenteredNestedStencil((1,2,3),(1,2,3),(1,2,3))
+        s_float = CenteredNestedStencil((1.,2.,3.),(1.,2.,3.),(1.,2.,3.))
+
+        v_int = rand(1:10,10);
+        v_float = rand(10);
+
+        c_int = rand(1:10,10);
+        c_float = rand(10);
+
+        @inferred SbpOperators.apply_stencil(s_int,   c_int, v_int,   2)
+        @inferred SbpOperators.apply_stencil(s_float, c_int, v_float, 2)
+        @inferred SbpOperators.apply_stencil(s_int,   c_int, v_float, 2)
+        @inferred SbpOperators.apply_stencil(s_float, c_int, v_int,   2)
+
+        @inferred SbpOperators.apply_stencil(s_int,   c_float, v_int,   2)
+        @inferred SbpOperators.apply_stencil(s_float, c_float, v_float, 2)
+        @inferred SbpOperators.apply_stencil(s_int,   c_float, v_float, 2)
+        @inferred SbpOperators.apply_stencil(s_float, c_float, v_int,   2)
+
+        @inferred SbpOperators.apply_stencil_backwards(s_int,   c_int, v_int,   2)
+        @inferred SbpOperators.apply_stencil_backwards(s_float, c_int, v_float, 2)
+        @inferred SbpOperators.apply_stencil_backwards(s_int,   c_int, v_float, 2)
+        @inferred SbpOperators.apply_stencil_backwards(s_float, c_int, v_int,   2)
+
+        @inferred SbpOperators.apply_stencil_backwards(s_int,   c_float, v_int,   2)
+        @inferred SbpOperators.apply_stencil_backwards(s_float, c_float, v_float, 2)
+        @inferred SbpOperators.apply_stencil_backwards(s_int,   c_float, v_float, 2)
+        @inferred SbpOperators.apply_stencil_backwards(s_float, c_float, v_int,   2)
+    end
+
+end