diff -r e6d8fd5e8268 -r 069e58fb3829 src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl
--- a/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl	Mon Feb 07 20:31:02 2022 +0100
+++ b/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl	Mon Feb 07 20:38:07 2022 +0100
@@ -30,35 +30,80 @@
 
 Implements the one-dimensional second derivative with variable coefficients.
 """
-struct SecondDerivativeVariable{T,N,M,K,TArray<:AbstractVector} <: TensorMapping{T,1,1}
+struct SecondDerivativeVariable{Dir,T,D,N,M,K,TArray<:AbstractArray} <: TensorMapping{T,D,D}
     inner_stencil::NestedStencil{T,N}
     closure_stencils::NTuple{M,NestedStencil{T,K}}
-    size::NTuple{1,Int}
+    size::NTuple{D,Int}
     coefficient::TArray
+
+    function SecondDerivativeVariable{Dir, D}(inner_stencil::NestedStencil{T,N}, closure_stencils::NTuple{M,NestedStencil{T,K}}, size::NTuple{D,Int}, coefficient::TArray) where {Dir,T,D,N,M,K,TArray<:AbstractArray}
+        return new{Dir,T,D,N,M,K,TArray}(inner_stencil,closure_stencils,size, coefficient)
+    end
+end
+
+function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, inner_stencil, closure_stencils, dir)
+    return SecondDerivativeVariable{dir, dimension(grid)}(inner_stencil, Tuple(closure_stencils), size(grid), coeff)
 end
 
 function SecondDerivativeVariable(grid::EquidistantGrid{1}, coeff::AbstractVector, inner_stencil, closure_stencils)
-    return SecondDerivativeVariable(inner_stencil, Tuple(closure_stencils), size(grid), coeff)
+    return SecondDerivativeVariable(grid, coeff, inner_stencil, closure_stencils, 1)
 end
 
+derivative_direction(::SecondDerivativeVariable{Dir}) where {Dir} = Dir
+
 closure_size(::SecondDerivativeVariable{T,N,M}) where {T,N,M} = M
 
 LazyTensors.range_size(op::SecondDerivativeVariable) = op.size
 LazyTensors.domain_size(op::SecondDerivativeVariable) = op.size
 
-function LazyTensors.apply(op::SecondDerivativeVariable{T}, v::AbstractVector{T}, i::Index{Lower}) where T
-    return @inbounds apply_stencil(op.closure_stencils[Int(i)], op.coefficient, v, Int(i))
+
+function derivative_view(op, a, I)
+    d = derivative_direction(op)
+
+    Iview = Base.setindex(I,:,d)
+    return @view a[Iview...]
+
+    # D = domain_dim(op)
+    # Iₗ, _, Iᵣ = split_tuple(I, Val(d-1), Val(1),  Val(D-d))
+    # return @view a[Iₗ..., :, Iᵣ...]
+end
+
+function apply_lower(op::SecondDerivativeVariable, v, I...)
+    ṽ = derivative_view(op, v, I)
+    c̃ = derivative_view(op, op.coefficient, I)
+
+    i = I[derivative_direction(op)]
+    return @inbounds apply_stencil(op.closure_stencils[i], c̃, ṽ, i)
 end
 
-function LazyTensors.apply(op::SecondDerivativeVariable{T}, v::AbstractVector{T}, i::Index{Interior}) where T
-    return apply_stencil(op.inner_stencil, op.coefficient, v, Int(i))
+function apply_interior(op::SecondDerivativeVariable, v, I...)
+    ṽ = derivative_view(op, v, I)
+    c̃ = derivative_view(op, op.coefficient, I)
+
+    i = I[derivative_direction(op)]
+    return apply_stencil(op.inner_stencil, c̃, ṽ, i)
 end
 
-function LazyTensors.apply(op::SecondDerivativeVariable{T}, v::AbstractVector{T}, i::Index{Upper}) where T
-    return @inbounds apply_stencil_backwards(op.closure_stencils[op.size[1]-Int(i)+1], op.coefficient, v, Int(i))
+function apply_upper(op::SecondDerivativeVariable, v, I...)
+    ṽ = derivative_view(op, v, I)
+    c̃ = derivative_view(op, op.coefficient, I)
+
+    i = I[derivative_direction(op)]
+    return @inbounds apply_stencil_backwards(op.closure_stencils[op.size[1]-i+1], c̃, ṽ, i)
 end
 
-function LazyTensors.apply(op::SecondDerivativeVariable{T}, v::AbstractVector{T}, i) where T
+function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractVector, I::Vararg{Index})
+    if I[derivative_direction(op)] isa Index{Lower}
+        return apply_lower(op, v, Int.(I)...)
+    elseif I[derivative_direction(op)] isa Index{Upper}
+        return apply_upper(op, v, Int.(I)...)
+    else
+        return apply_interior(op, v, Int.(I)...)
+    end
+end
+
+function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractVector, I...)
+    i = I[derivative_direction(op)]
     r = getregion(i, closure_size(op), op.size[1])
     return LazyTensors.apply(op, v, Index(i, r))
 end

diff -r e6d8fd5e8268 -r 069e58fb3829 test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl
--- a/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl	Mon Feb 07 20:31:02 2022 +0100
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl	Mon Feb 07 20:38:07 2022 +0100
@@ -8,17 +8,15 @@
 
 
 @testset "SecondDerivativeVariable" begin
-    g = EquidistantGrid(11, 0., 1.)
-    c = [  1,  3,  6, 10, 15, 21, 28, 36, 45, 55, 66]
-
     interior_stencil = CenteredNestedStencil((1/2, 1/2, 0.),(-1/2, -1., -1/2),( 0., 1/2, 1/2))
     closure_stencils = [
         NestedStencil(( 2.,  -1., 0.),(-3., 1.,  0.), (1., 0., 0.), center = 1),
     ]
 
     @testset "1D" begin
+        g = EquidistantGrid(11, 0., 1.)
+        c = [  1,  3,  6, 10, 15, 21, 28, 36, 45, 55, 66]
         @testset "Constructors" begin
-            @test SecondDerivativeVariable(interior_stencil,Tuple(closure_stencils), (4,), c) isa TensorMapping
             @test SecondDerivativeVariable(g, c, interior_stencil, closure_stencils) isa TensorMapping
         end
 
@@ -44,6 +42,21 @@
             @test apply_to_functions(v=x->1., c=x->-x) == zeros(11)
             @test apply_to_functions(v=x->x, c=x-> 1.) == zeros(11)
             @test apply_to_functions(v=x->x, c=x->-x) == -ones(11)
+            @test apply_to_functions(v=x->x^2, c=x->1.) == 2ones(11)
+        end
+    end
+
+    @testset "2D" begin
+        g = EquidistantGrid((11,9), (0.,0.), (10.,8.)) # h = 1
+        c = evalOn(g, (x,y)->x+y)
+        @testset "Constructors" begin
+            @test SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,1) isa TensorMapping
+        end
+
+        @testset "sizes" begin
+        end
+
+        @testset "application" begin
         end
     end
 end