Mercurial > repos > public > sbplib_julia

changeset 562:8f7919a9b398 feature/boundary_ops

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 30 Nov 2020 18:30:24 +0100
parents 884be64e82d9 (current diff) d1929491180b (diff)
children 212e266043dd
files src/LazyTensors/lazy_tensor_operations.jl src/SbpOperators/boundaryops/boundary_restriction.jl test/testSbpOperators.jl
diffstat 14 files changed, 326 insertions(+), 213 deletions(-) [+]
line wrap: on
line diff
--- a/TODO.md	Thu Nov 26 09:03:54 2020 +0100
+++ b/TODO.md	Mon Nov 30 18:30:24 2020 +0100
@@ -10,10 +10,10 @@
  - [ ] Create a struct that bundles the necessary Tensor operators for solving the wave equation.
  - [ ] Add a quick and simple way of running all tests for all subpackages.
  - [ ] Replace getindex hack for flatteing tuples with flatten_tuple.
- - [ ] Fix indexing signatures. We should make sure we are not too specific. For the "inbetween" layers we don't know what type of index is coming so we should use `I...` instead of `I::Vararg{Int,R}` or probably better `I::Vararg{Any,R}`
  - [ ] Use `@inferred` in a lot of tests.
  - [ ] Make sure we are setting tolerances in tests in a consistent way
  - [ ] Add check for correct domain sizes to lazy tensor operations using SizeMismatch
+ - [ ] Write down some coding guideline or checklist for code convetions. For example i,j,... för indecies and I for multi-index
 
 ## Repo
  - [ ] Add Vidar to the authors list
--- a/src/Grids/Grids.jl	Thu Nov 26 09:03:54 2020 +0100
+++ b/src/Grids/Grids.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -7,7 +7,7 @@
 abstract type BoundaryIdentifier end
 struct CartesianBoundary{Dim, R<:Region} <: BoundaryIdentifier end
 dim(::CartesianBoundary{Dim, R}) where {Dim, R} = Dim
-region(::CartesianBoundary{Dim, R}) where {Dim, R} = R
+region(::CartesianBoundary{Dim, R}) where {Dim, R} = R  #TODO: Should return R()
 
 export dim, region
 
--- a/src/LazyTensors/lazy_tensor_operations.jl	Thu Nov 26 09:03:54 2020 +0100
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -14,10 +14,7 @@
 # TODO: Do boundschecking on creation!
 export LazyTensorMappingApplication
 
-# TODO: Go through and remove unneccerary type parameters on functions
-Base.getindex(ta::LazyTensorMappingApplication{T,0}, I::Index) where T = apply(ta.t, ta.o, I)
-Base.getindex(ta::LazyTensorMappingApplication{T,R}, I::Vararg{Index,R}) where {T,R} = apply(ta.t, ta.o, I...)
-Base.getindex(ta::LazyTensorMappingApplication{T,R}, I::Vararg{Int,R}) where {T,R} = apply(ta.t, ta.o, Index{Unknown}.(I)...)
+Base.getindex(ta::LazyTensorMappingApplication{T,R}, I::Vararg{Any,R}) where {T,R} = apply(ta.t, ta.o, I...)
 Base.size(ta::LazyTensorMappingApplication) = range_size(ta.t)
 # TODO: What else is needed to implement the AbstractArray interface?
 
@@ -50,8 +47,8 @@
 Base.adjoint(tm::TensorMapping) = LazyTensorMappingTranspose(tm)
 Base.adjoint(tmt::LazyTensorMappingTranspose) = tmt.tm
 
-apply(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Index,D}) where {T,R,D} = apply_transpose(tmt.tm, v, I...)
-apply_transpose(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Index,R}) where {T,R,D} = apply(tmt.tm, v, I...)
+apply(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,D} = apply_transpose(tmt.tm, v, I...)
+apply_transpose(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmt.tm, v, I...)
 
 range_size(tmt::LazyTensorMappingTranspose) = domain_size(tmt.tm)
 domain_size(tmt::LazyTensorMappingTranspose) = range_size(tmt.tm)
@@ -67,11 +64,11 @@
 end
 # TODO: Boundschecking in constructor.
 
-apply(tmBinOp::LazyTensorMappingBinaryOperation{:+,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Index,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) + apply(tmBinOp.tm2, v, I...)
-apply(tmBinOp::LazyTensorMappingBinaryOperation{:-,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Index,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) - apply(tmBinOp.tm2, v, I...)
+apply(tmBinOp::LazyTensorMappingBinaryOperation{:+,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) + apply(tmBinOp.tm2, v, I...)
+apply(tmBinOp::LazyTensorMappingBinaryOperation{:-,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) - apply(tmBinOp.tm2, v, I...)
 
-range_size(tmBinOp::LazyTensorMappingBinaryOperation{Op,T,R,D}) where {Op,T,R,D} = range_size(tmBinOp.tm1)
-domain_size(tmBinOp::LazyTensorMappingBinaryOperation{Op,T,R,D}) where {Op,T,R,D} = domain_size(tmBinOp.tm1)
+range_size(tmBinOp::LazyTensorMappingBinaryOperation) = range_size(tmBinOp.tm1)
+domain_size(tmBinOp::LazyTensorMappingBinaryOperation) = domain_size(tmBinOp.tm1)
 
 Base.:+(tm1::TensorMapping{T,R,D}, tm2::TensorMapping{T,R,D}) where {T,R,D} = LazyTensorMappingBinaryOperation{:+,T,R,D}(tm1,tm2)
 Base.:-(tm1::TensorMapping{T,R,D}, tm2::TensorMapping{T,R,D}) where {T,R,D} = LazyTensorMappingBinaryOperation{:-,T,R,D}(tm1,tm2)
@@ -95,11 +92,11 @@
 range_size(tm::TensorMappingComposition) = range_size(tm.t1)
 domain_size(tm::TensorMappingComposition) = domain_size(tm.t2)
 
-function apply(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{T,D}, I::Vararg{S,R} where S) where {T,R,K,D}
+function apply(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,K,D}
     apply(c.t1, c.t2*v, I...)
 end
 
-function apply_transpose(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{T,R}, I::Vararg{S,D} where S) where {T,R,K,D}
+function apply_transpose(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,K,D}
     apply_transpose(c.t2, c.t1'*v, I...)
 end
 
@@ -132,7 +129,7 @@
 range_size(llm::LazyLinearMap) = size(llm.A)[[llm.range_indicies...]]
 domain_size(llm::LazyLinearMap) = size(llm.A)[[llm.domain_indicies...]]
 
-function apply(llm::LazyLinearMap{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Index,R}) where {T,R,D}
+function apply(llm::LazyLinearMap{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D}
     view_index = ntuple(i->:,ndims(llm.A))
     for i ∈ 1:R
         view_index = Base.setindex(view_index, Int(I[i]), llm.range_indicies[i])
@@ -141,7 +138,7 @@
     return sum(A_view.*v)
 end
 
-function apply_transpose(llm::LazyLinearMap{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Index,D}) where {T,R,D}
+function apply_transpose(llm::LazyLinearMap{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,D}
     apply(LazyLinearMap(llm.A, llm.domain_indicies, llm.range_indicies), v, I...)
 end
 
@@ -240,8 +237,6 @@
 # Resolve ambiguity between the two previous methods
 InflatedTensorMapping(I1::IdentityMapping{T}, I2::IdentityMapping{T}) where T = InflatedTensorMapping(I1,I2,IdentityMapping{T}())
 
-# TODO: Implement syntax and constructors for products of different combinations of InflatedTensorMapping and IdentityMapping
-
 # TODO: Implement some pretty printing in terms of ⊗. E.g InflatedTensorMapping(I(3),B,I(2)) -> I(3)⊗B⊗I(2)
 
 function range_size(itm::InflatedTensorMapping)
@@ -261,30 +256,56 @@
 end
 
 function apply(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D}
-    view_index, inner_index = split_index(itm, I...)
+    dim_before = range_dim(itm.before)
+    dim_domain = domain_dim(itm.tm)
+    dim_range = range_dim(itm.tm)
+    dim_after = range_dim(itm.after)
+
+    view_index, inner_index = split_index(Val(dim_before), Val(dim_domain), Val(dim_range), Val(dim_after), I...)
 
     v_inner = view(v, view_index...)
     return apply(itm.tm, v_inner, inner_index...)
 end
 
+function apply_transpose(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,D}
+    dim_before = range_dim(itm.before)
+    dim_domain = domain_dim(itm.tm)
+    dim_range = range_dim(itm.tm)
+    dim_after = range_dim(itm.after)
+
+    view_index, inner_index = split_index(Val(dim_before), Val(dim_range), Val(dim_domain), Val(dim_after), I...)
+
+    v_inner = view(v, view_index...)
+    return apply_transpose(itm.tm, v_inner, inner_index...)
+end
+
 
 """
-    split_index(...)
+    split_index(::Val{dim_before}, ::Val{dim_view}, ::Val{dim_index}, ::Val{dim_after}, I...)
 
-Splits the multi-index into two parts. One part for the view that the inner TensorMapping acts on, and one part for indexing the result
+Splits the multi-index `I` into two parts. One part which is expected to be
+used as a view, and one which is expected to be used as an index.
 Eg.
 ```
-(1,2,3,4) -> (1,:,:,4), (2,3)
+split_index(Val(1),Val(3),Val(2),Val(1),(1,2,3,4)) -> (1,:,:,:,4), (2,3)
 ```
+
+`dim_view` controls how many colons are in the view, and `dim_index` controls
+how many elements are extracted from the middle.
+`dim_before` and `dim_after` decides the length of the index parts before and after the colons in the view index.
+
+Arguments should satisfy `length(I) == dim_before+B_domain+dim_after`.
+
+The returned values satisfy
+ * `length(view_index) == dim_before + dim_view + dim_after`
+ * `length(I_middle) == dim_index`
 """
-function split_index(itm::InflatedTensorMapping{T,R,D}, I::Vararg{Any,R}) where {T,R,D}
-    I_before = slice_tuple(I, Val(1), Val(range_dim(itm.before)))
-    I_after = slice_tuple(I, Val(R-range_dim(itm.after)+1), Val(R))
+function split_index(::Val{dim_before}, ::Val{dim_view}, ::Val{dim_index}, ::Val{dim_after}, I...) where {dim_before,dim_view, dim_index,dim_after}
+    I_before, I_middle, I_after = split_tuple(I, Val(dim_before), Val(dim_index))
 
-    view_index = (I_before..., ntuple((i)->:,domain_dim(itm.tm))..., I_after...)
-    inner_index = slice_tuple(I, Val(range_dim(itm.before)+1), Val(R-range_dim(itm.after)))
+    view_index = (I_before..., ntuple((i)->:, dim_view)..., I_after...)
 
-    return (view_index, inner_index)
+    return view_index, I_middle
 end
 
 # TODO: Can this be replaced by something more elegant while still being type stable? 2020-10-21
@@ -302,6 +323,32 @@
 end
 
 """
+    split_tuple(t::Tuple{...}, ::Val{M}) where {N,M}
+
+Split the tuple `t` into two parts. the first part is `M` long.
+E.g
+```
+split_tuple((1,2,3,4),Val(3)) -> (1,2,3), (4,)
+```
+"""
+function split_tuple(t::NTuple{N},::Val{M}) where {N,M}
+    return slice_tuple(t,Val(1), Val(M)), slice_tuple(t,Val(M+1), Val(N))
+end
+
+"""
+    split_tuple(t::Tuple{...},::Val{M},::Val{K}) where {N,M,K}
+
+Same as `split_tuple(t::NTuple{N},::Val{M})` but splits the tuple in three parts. With the first
+two parts having lenght `M` and `K`.
+"""
+function split_tuple(t::NTuple{N},::Val{M},::Val{K}) where {N,M,K}
+    p1, tail = split_tuple(t, Val(M))
+    p2, p3 = split_tuple(tail, Val(K))
+    return p1,p2,p3
+end
+
+
+"""
     flatten_tuple(t)
 
 Takes a nested tuple and flattens the whole structure
--- a/src/RegionIndices/RegionIndices.jl	Thu Nov 26 09:03:54 2020 +0100
+++ b/src/RegionIndices/RegionIndices.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -4,9 +4,8 @@
 struct Interior <: Region end
 struct Lower    <: Region end
 struct Upper    <: Region end
-struct Unknown  <: Region end
 
-export Region, Interior, Lower, Upper, Unknown
+export Region, Interior, Lower, Upper
 
 struct Index{R<:Region, T<:Integer}
     i::T
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Thu Nov 26 09:03:54 2020 +0100
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -41,24 +41,18 @@
 LazyTensors.range_size(e::BoundaryRestriction) = ()
 LazyTensors.domain_size(e::BoundaryRestriction) = e.size
 
-# TODO: Currently not working.
-# We need to handle getindex for LazyTensorMappingApplication such that we pass more #indices than the
-# range size of the TensorMapping. Or we need to be able to handle the case where we dont pass any index, for
-# 0-dimensional tensormappings.
+# TODO: Should we support indexing into the 0-dimensional lazyarray? This is
+# supported for arrays with linear index style (i.e for e.g
+# u = fill(1), u[] and u[1] are both valid.) This currently not supported by
+# LazyTensorMappingApplication.
 " Restricts a grid function v on a grid of size m to the scalar element v[1]"
-function LazyTensors.apply(e::BoundaryRestriction{T,M,Lower}, v::AbstractVector{T}, i::Index{Lower}) where {T,M}
-    @boundscheck if Int(i)!=1
-        throw(BoundsError())
-    end
-    apply_stencil(e.stencil,v,Int(i))
+function LazyTensors.apply(e::BoundaryRestriction{T,M,Lower}, v::AbstractVector{T}) where {T,M}
+    apply_stencil(e.stencil,v,1)
 end
 
 " Restricts a grid function v on a grid of size m to the scalar element v[m]"
-function LazyTensors.apply(e::BoundaryRestriction{T,M,Upper}, v::AbstractVector{T}, i::Index{Upper}) where {T,M}
-    @boundscheck if Int(i) != e.size[1]
-        throw(BoundsError())
-    end
-    apply_stencil_backwards(e.stencil,v,Int(i))
+function LazyTensors.apply(e::BoundaryRestriction{T,M,Upper}, v::AbstractVector{T}) where {T,M}
+    apply_stencil_backwards(e.stencil,v,e.size[1])
 end
 
 " Transpose of a restriction is an inflation or prolongation.
--- a/src/SbpOperators/constantstenciloperator.jl	Thu Nov 26 09:03:54 2020 +0100
+++ b/src/SbpOperators/constantstenciloperator.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -14,7 +14,7 @@
     return @inbounds h_inv*h_inv*Int(op.parity)*apply_stencil_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))
 end
 
-@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, index::Index{Unknown})
+@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i)
     N = length(v)
     r = getregion(Int(index), closuresize(op), N)
     i = Index(Int(index), r)
@@ -26,10 +26,9 @@
 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Upper}, N::Integer) where T = v*h*op.quadratureClosure[N-Int(i)+1]
 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Interior}, N::Integer) where T = v*h
 
-function apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, index::Index{Unknown}, N::Integer) where T
-    r = getregion(Int(index), closuresize(op), N)
-    i = Index(Int(index), r)
-    return apply_quadrature(op, h, v, i, N)
+function apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i, N::Integer) where T
+    r = getregion(i, closuresize(op), N)
+    return apply_quadrature(op, h, v, Index(i, r), N)
 end
 export apply_quadrature
 
@@ -38,10 +37,9 @@
 apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Index{Upper}, N::Integer) where T = h_inv*v/op.quadratureClosure[N-Int(i)+1]
 apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Index{Interior}, N::Integer) where T = v*h_inv
 
-function apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, index::Index{Unknown}, N::Integer) where T
-    r = getregion(Int(index), closuresize(op), N)
-    i = Index(Int(index), r)
-    return apply_inverse_quadrature(op, h_inv, v, i, N)
+function apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i, N::Integer) where T
+    r = getregion(i, closuresize(op), N)
+    return apply_inverse_quadrature(op, h_inv, v, Index(i, r), N)
 end
 
 export apply_inverse_quadrature
@@ -62,14 +60,14 @@
 
 export apply_normal_derivative_transpose
 
-function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, i::Index, N::Integer, ::Type{Lower})
+function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, i, N::Integer, ::Type{Lower})
     @boundscheck if !(0<length(Int(i)) <= N)
         throw(BoundsError())
     end
     h_inv*op.dClosure[Int(i)-1]*v
 end
 
-function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, i::Index, N::Integer, ::Type{Upper})
+function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, i, N::Integer, ::Type{Upper})
     @boundscheck if !(0<length(Int(i)) <= N)
         throw(BoundsError())
     end
--- a/src/SbpOperators/laplace/laplace.jl	Thu Nov 26 09:03:54 2020 +0100
+++ b/src/SbpOperators/laplace/laplace.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -23,24 +23,24 @@
 LazyTensors.range_size(L::Laplace) = getindex.(range_size.(L.D2),1)
 LazyTensors.domain_size(L::Laplace) = getindex.(domain_size.(L.D2),1)
 
-function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim}
+function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {T,Dim}
     error("not implemented")
 end
 
 # u = L*v
-function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, I::Index) where T
-    @inbounds u = LazyTensors.apply(L.D2[1],v,I)
+function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, i) where T
+    @inbounds u = LazyTensors.apply(L.D2[1],v,i)
     return u
 end
 
-function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::Index, J::Index) where T
+function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, i, j) where T
     # 2nd x-derivative
-    @inbounds vx = view(v, :, Int(J))
-    @inbounds uᵢ = LazyTensors.apply(L.D2[1], vx , I)
+    @inbounds vx = view(v, :, Int(j))
+    @inbounds uᵢ = LazyTensors.apply(L.D2[1], vx , i)
 
     # 2nd y-derivative
-    @inbounds vy = view(v, Int(I), :)
-    @inbounds uᵢ += LazyTensors.apply(L.D2[2], vy , J)
+    @inbounds vy = view(v, Int(i), :)
+    @inbounds uᵢ += LazyTensors.apply(L.D2[2], vy , j)
 
     return uᵢ
 end
--- a/src/SbpOperators/laplace/secondderivative.jl	Thu Nov 26 09:03:54 2020 +0100
+++ b/src/SbpOperators/laplace/secondderivative.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -20,29 +20,24 @@
 LazyTensors.range_size(D2::SecondDerivative) = D2.size
 LazyTensors.domain_size(D2::SecondDerivative) = D2.size
 
-#TODO: The 1D tensor mappings should not have to dispatch on 1D tuples if we write LazyTensor.apply for vararg right?!?!
-#      Currently have to index the Tuple{Index} in each method in order to call the stencil methods which is ugly.
-#      I thought I::Vararg{Index,R} fell back to just Index for R = 1
+# Apply for different regions Lower/Interior/Upper or Unknown region
+function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, i::Index{Lower}) where T
+    return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.closureStencils[Int(i)], v, Int(i))
+end
 
-# Apply for different regions Lower/Interior/Upper or Unknown region
-function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Lower}) where T
-    return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.closureStencils[Int(I)], v, Int(I))
+function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, i::Index{Interior}) where T
+    return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.innerStencil, v, Int(i))
 end
 
-function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Interior}) where T
-    return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.innerStencil, v, Int(I))
+function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, i::Index{Upper}) where T
+    N = length(v) # TODO: Use domain_size here instead? N = domain_size(D2,size(v))
+    return @inbounds D2.h_inv*D2.h_inv*apply_stencil_backwards(D2.closureStencils[N-Int(i)+1], v, Int(i))
 end
 
-function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Upper}) where T
-    N = length(v) # TODO: Use domain_size here instead? N = domain_size(D2,size(v))
-    return @inbounds D2.h_inv*D2.h_inv*apply_stencil_backwards(D2.closureStencils[N-Int(I)+1], v, Int(I))
-end
-
-function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, index::Index{Unknown}) where T
+function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, i) where T
     N = length(v)  # TODO: Use domain_size here instead?
-    r = getregion(Int(index), closuresize(D2), N)
-    I = Index(Int(index), r)
-    return LazyTensors.apply(D2, v, I)
+    r = getregion(i, closuresize(D2), N)
+    return LazyTensors.apply(D2, v, Index(i, r))
 end
 
 closuresize(D2::SecondDerivative{T,N,M,K}) where {T<:Real,N,M,K} = M
--- a/src/SbpOperators/quadrature/diagonal_inner_product.jl	Thu Nov 26 09:03:54 2020 +0100
+++ b/src/SbpOperators/quadrature/diagonal_inner_product.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -17,30 +17,25 @@
 LazyTensors.range_size(H::DiagonalInnerProduct) = H.size
 LazyTensors.domain_size(H::DiagonalInnerProduct) = H.size
 
-function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index) where T
-    return @inbounds apply(H, v, I)
+function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Lower}) where T
+    return @inbounds H.h*H.quadratureClosure[Int(i)]*v[Int(i)]
 end
 
-function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index{Lower}) where T
-    return @inbounds H.h*H.quadratureClosure[Int(I)]*v[Int(I)]
-end
-
-function LazyTensors.apply(H::DiagonalInnerProduct{T},v::AbstractVector{T}, I::Index{Upper}) where T
+function LazyTensors.apply(H::DiagonalInnerProduct{T},v::AbstractVector{T}, i::Index{Upper}) where T
     N = length(v);
-    return @inbounds H.h*H.quadratureClosure[N-Int(I)+1]*v[Int(I)]
+    return @inbounds H.h*H.quadratureClosure[N-Int(i)+1]*v[Int(i)]
 end
 
-function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index{Interior}) where T
-    return @inbounds H.h*v[Int(I)]
+function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Interior}) where T
+    return @inbounds H.h*v[Int(i)]
 end
 
-function LazyTensors.apply(H::DiagonalInnerProduct{T},  v::AbstractVector{T}, index::Index{Unknown}) where T
+function LazyTensors.apply(H::DiagonalInnerProduct{T},  v::AbstractVector{T}, i) where T
     N = length(v);
-    r = getregion(Int(index), closuresize(H), N)
-    i = Index(Int(index), r)
-    return LazyTensors.apply(H, v, i)
+    r = getregion(i, closuresize(H), N)
+    return LazyTensors.apply(H, v, Index(i, r))
 end
 
-LazyTensors.apply_transpose(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index) where T = LazyTensors.apply(H,v,I)
+LazyTensors.apply_transpose(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, i) where T = LazyTensors.apply(H,v,i)
 
 closuresize(H::DiagonalInnerProduct{T,M}) where {T,M} = M
--- a/src/SbpOperators/quadrature/inverse_diagonal_inner_product.jl	Thu Nov 26 09:03:54 2020 +0100
+++ b/src/SbpOperators/quadrature/inverse_diagonal_inner_product.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -18,27 +18,26 @@
 LazyTensors.domain_size(Hi::InverseDiagonalInnerProduct) = Hi.size
 
 
-function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index{Lower}) where T
-    return @inbounds Hi.h_inv*Hi.inverseQuadratureClosure[Int(I)]*v[Int(I)]
+function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Lower}) where T
+    return @inbounds Hi.h_inv*Hi.inverseQuadratureClosure[Int(i)]*v[Int(i)]
 end
 
-function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index{Upper}) where T
+function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Upper}) where T
     N = length(v);
-    return @inbounds Hi.h_inv*Hi.inverseQuadratureClosure[N-Int(I)+1]*v[Int(I)]
+    return @inbounds Hi.h_inv*Hi.inverseQuadratureClosure[N-Int(i)+1]*v[Int(i)]
 end
 
-function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index{Interior}) where T
-    return @inbounds Hi.h_inv*v[Int(I)]
+function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Interior}) where T
+    return @inbounds Hi.h_inv*v[Int(i)]
 end
 
-function LazyTensors.apply(Hi::InverseDiagonalInnerProduct,  v::AbstractVector{T}, index::Index{Unknown}) where T
+function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T},  v::AbstractVector{T}, i) where T
     N = length(v);
-    r = getregion(Int(index), closuresize(Hi), N)
-    i = Index(Int(index), r)
-    return LazyTensors.apply(Hi, v, i)
+    r = getregion(i, closuresize(Hi), N)
+    return LazyTensors.apply(Hi, v, Index(i, r))
 end
 
-LazyTensors.apply_transpose(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index) where T = LazyTensors.apply(Hi,v,I)
+LazyTensors.apply_transpose(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, i) where T = LazyTensors.apply(Hi,v,i)
 
 
 closuresize(Hi::InverseDiagonalInnerProduct{T,M}) where {T,M} =  M
--- a/src/SbpOperators/quadrature/inverse_quadrature.jl	Thu Nov 26 09:03:54 2020 +0100
+++ b/src/SbpOperators/quadrature/inverse_quadrature.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -25,23 +25,23 @@
 
 LazyTensors.domain_size(Qi::InverseQuadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
 
-function LazyTensors.apply(Qi::InverseQuadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim}
+function LazyTensors.apply(Qi::InverseQuadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {T,Dim}
     error("not implemented")
 end
 
-@inline function LazyTensors.apply(Qi::InverseQuadrature{1,T}, v::AbstractVector{T}, I::Index) where T
-    @inbounds q = apply(Qi.Hi[1], v , I)
+@inline function LazyTensors.apply(Qi::InverseQuadrature{1,T}, v::AbstractVector{T}, i) where T
+    @inbounds q = apply(Qi.Hi[1], v , i)
     return q
 end
 
-@inline function LazyTensors.apply(Qi::InverseQuadrature{2,T}, v::AbstractArray{T,2}, I::Index, J::Index) where T
+@inline function LazyTensors.apply(Qi::InverseQuadrature{2,T}, v::AbstractArray{T,2}, i,j) where T
     # InverseQuadrature in x direction
-    @inbounds vx = view(v, :, Int(J))
-    @inbounds qx_inv = apply(Qi.Hi[1], vx , I)
+    @inbounds vx = view(v, :, Int(j))
+    @inbounds qx_inv = apply(Qi.Hi[1], vx , i)
     # InverseQuadrature in y-direction
-    @inbounds vy = view(v, Int(I), :)
-    @inbounds qy_inv = apply(Qi.Hi[2], vy, J)
+    @inbounds vy = view(v, Int(i), :)
+    @inbounds qy_inv = apply(Qi.Hi[2], vy, j)
     return qx_inv*qy_inv
 end
 
-LazyTensors.apply_transpose(Qi::InverseQuadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {Dim,T} = LazyTensors.apply(Qi,v,I...)
+LazyTensors.apply_transpose(Qi::InverseQuadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {Dim,T} = LazyTensors.apply(Qi,v,I...)
--- a/src/SbpOperators/quadrature/quadrature.jl	Thu Nov 26 09:03:54 2020 +0100
+++ b/src/SbpOperators/quadrature/quadrature.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -22,23 +22,23 @@
 LazyTensors.range_size(H::Quadrature) = getindex.(range_size.(H.H),1)
 LazyTensors.domain_size(H::Quadrature) = getindex.(domain_size.(H.H),1)
 
-function LazyTensors.apply(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim}
+function LazyTensors.apply(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {T,Dim}
     error("not implemented")
 end
 
-function LazyTensors.apply(Q::Quadrature{1,T}, v::AbstractVector{T}, I::Index) where T
-    @inbounds q = apply(Q.H[1], v , I)
+function LazyTensors.apply(Q::Quadrature{1,T}, v::AbstractVector{T}, i) where T
+    @inbounds q = apply(Q.H[1], v , i)
     return q
 end
 
-function LazyTensors.apply(Q::Quadrature{2,T}, v::AbstractArray{T,2}, I::Index, J::Index) where T
+function LazyTensors.apply(Q::Quadrature{2,T}, v::AbstractArray{T,2}, i, j) where T
     # Quadrature in x direction
-    @inbounds vx = view(v, :, Int(J))
-    @inbounds qx = apply(Q.H[1], vx , I)
+    @inbounds vx = view(v, :, Int(j))
+    @inbounds qx = apply(Q.H[1], vx , i)
     # Quadrature in y-direction
-    @inbounds vy = view(v, Int(I), :)
-    @inbounds qy = apply(Q.H[2], vy, J)
+    @inbounds vy = view(v, Int(i), :)
+    @inbounds qy = apply(Q.H[2], vy, j)
     return qx*qy
 end
 
-LazyTensors.apply_transpose(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {Dim,T} = LazyTensors.apply(Q,v,I...)
+LazyTensors.apply_transpose(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {Dim,T} = LazyTensors.apply(Q,v,I...)
--- a/test/testLazyTensors.jl	Thu Nov 26 09:03:54 2020 +0100
+++ b/test/testLazyTensors.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -8,10 +8,10 @@
 
 @testset "Generic Mapping methods" begin
     struct DummyMapping{T,R,D} <: TensorMapping{T,R,D} end
-    LazyTensors.apply(m::DummyMapping{T,R,D}, v, i::NTuple{R,Index{<:Region}}) where {T,R,D} = :apply
+    LazyTensors.apply(m::DummyMapping{T,R,D}, v, I::Vararg{Any,R}) where {T,R,D} = :apply
     @test range_dim(DummyMapping{Int,2,3}()) == 2
     @test domain_dim(DummyMapping{Int,2,3}()) == 3
-    @test apply(DummyMapping{Int,2,3}(), zeros(Int, (0,0,0)),(Index{Unknown}(0),Index{Unknown}(0))) == :apply
+    @test apply(DummyMapping{Int,2,3}(), zeros(Int, (0,0,0)),0,0) == :apply
     @test eltype(DummyMapping{Int,2,3}()) == Int
     @test eltype(DummyMapping{Float64,2,3}()) == Float64
 end
@@ -19,19 +19,18 @@
 @testset "Mapping transpose" begin
     struct DummyMapping{T,R,D} <: TensorMapping{T,R,D} end
 
-    LazyTensors.apply(m::DummyMapping{T,R,D}, v, I::Vararg{Index{<:Region},R}) where {T,R,D} = :apply
-    LazyTensors.apply_transpose(m::DummyMapping{T,R,D}, v, I::Vararg{Index{<:Region},D}) where {T,R,D} = :apply_transpose
+    LazyTensors.apply(m::DummyMapping{T,R}, v, I::Vararg{Any,R}) where {T,R} = :apply
+    LazyTensors.apply_transpose(m::DummyMapping{T,R,D}, v, I::Vararg{Any,D}) where {T,R,D} = :apply_transpose
 
-    LazyTensors.range_size(m::DummyMapping{T,R,D}) where {T,R,D} = :range_size
-    LazyTensors.domain_size(m::DummyMapping{T,R,D}) where {T,R,D} = :domain_size
+    LazyTensors.range_size(m::DummyMapping) = :range_size
+    LazyTensors.domain_size(m::DummyMapping) = :domain_size
 
     m = DummyMapping{Float64,2,3}()
-    I = Index{Unknown}(0)
     @test m' isa TensorMapping{Float64, 3,2}
     @test m'' == m
-    @test apply(m',zeros(Float64,(0,0)), I, I, I) == :apply_transpose
-    @test apply(m'',zeros(Float64,(0,0,0)), I, I) == :apply
-    @test apply_transpose(m', zeros(Float64,(0,0,0)), I, I) == :apply
+    @test apply(m',zeros(Float64,(0,0)), 0, 0, 0) == :apply_transpose
+    @test apply(m'',zeros(Float64,(0,0,0)), 0, 0) == :apply
+    @test apply_transpose(m', zeros(Float64,(0,0,0)), 0, 0) == :apply
 
     @test range_size(m') == :domain_size
     @test domain_size(m') == :range_size
@@ -42,7 +41,7 @@
         domain_size::NTuple{D,Int}
     end
 
-    LazyTensors.apply(m::SizeDoublingMapping{T,R,D}, v, i::Vararg{Index{<:Region},R}) where {T,R,D} = (:apply,v,i)
+    LazyTensors.apply(m::SizeDoublingMapping{T,R}, v, i::Vararg{Any,R}) where {T,R} = (:apply,v,i)
     LazyTensors.range_size(m::SizeDoublingMapping) = 2 .* m.domain_size
     LazyTensors.domain_size(m::SizeDoublingMapping) = m.domain_size
 
@@ -51,15 +50,11 @@
     v = [0,1,2]
     @test m*v isa AbstractVector{Int}
     @test size(m*v) == 2 .*size(v)
-    @test (m*v)[Index{Upper}(0)] == (:apply,v,(Index{Upper}(0),))
-    @test (m*v)[0] == (:apply,v,(Index{Unknown}(0),))
+    @test (m*v)[0] == (:apply,v,(0,))
     @test m*m*v isa AbstractVector{Int}
-    @test (m*m*v)[Index{Upper}(1)] == (:apply,m*v,(Index{Upper}(1),))
-    @test (m*m*v)[1] == (:apply,m*v,(Index{Unknown}(1),))
-    @test (m*m*v)[Index{Interior}(3)] == (:apply,m*v,(Index{Interior}(3),))
-    @test (m*m*v)[3] == (:apply,m*v,(Index{Unknown}(3),))
-    @test (m*m*v)[Index{Lower}(6)] == (:apply,m*v,(Index{Lower}(6),))
-    @test (m*m*v)[6] == (:apply,m*v,(Index{Unknown}(6),))
+    @test (m*m*v)[1] == (:apply,m*v,(1,))
+    @test (m*m*v)[3] == (:apply,m*v,(3,))
+    @test (m*m*v)[6] == (:apply,m*v,(6,))
     @test_broken BoundsError == (m*m*v)[0]
     @test_broken BoundsError == (m*m*v)[7]
     @test_throws MethodError m*m
@@ -70,16 +65,15 @@
 
     m = SizeDoublingMapping{Float64, 2, 2}((3,3))
     v = ones(3,3)
-    I = (Index{Lower}(1),Index{Interior}(2));
     @test size(m*v) == 2 .*size(v)
-    @test (m*v)[I] == (:apply,v,I)
+    @test (m*v)[1,2] == (:apply,v,(1,2))
 
     struct ScalingOperator{T,D} <: TensorMapping{T,D,D}
         λ::T
         size::NTuple{D,Int}
     end
 
-    LazyTensors.apply(m::ScalingOperator{T,D}, v, I::Vararg{Index,D}) where {T,D} = m.λ*v[I]
+    LazyTensors.apply(m::ScalingOperator{T,D}, v, I::Vararg{Any,D}) where {T,D} = m.λ*v[I...]
     LazyTensors.range_size(m::ScalingOperator) = m.size
     LazyTensors.domain_size(m::ScalingOperator) = m.size
 
@@ -91,8 +85,7 @@
     m = ScalingOperator{Int,2}(2,(2,2))
     v = [[1 2];[3 4]]
     @test m*v == [[2 4];[6 8]]
-    I = (Index{Upper}(2),Index{Lower}(1))
-    @test (m*v)[I] == 6
+    @test (m*v)[2,1] == 6
 end
 
 @testset "TensorMapping binary operations" begin
@@ -102,7 +95,7 @@
         domain_size::NTuple{D,Int}
     end
 
-    LazyTensors.apply(m::ScalarMapping{T,R,D}, v, I::Vararg{Index{<:Region}}) where {T,R,D} = m.λ*v[I...]
+    LazyTensors.apply(m::ScalarMapping{T,R}, v, I::Vararg{Any,R}) where {T,R} = m.λ*v[I...]
     LazyTensors.range_size(m::ScalarMapping) = m.domain_size
     LazyTensors.domain_size(m::ScalarMapping) = m.range_size
 
@@ -340,61 +333,115 @@
     B = LazyLinearMap(B̃,(1,2),(3,))
     C = LazyLinearMap(C̃,(1,),(2,3))
 
-    @test InflatedTensorMapping(I(3,2), A, I(4)) isa TensorMapping{Float64, 4, 4}
-    @test InflatedTensorMapping(I(3,2), B, I(4)) isa TensorMapping{Float64, 5, 4}
-    @test InflatedTensorMapping(I(3), C, I(2,3)) isa TensorMapping{Float64, 4, 5}
-    @test InflatedTensorMapping(C, I(2,3)) isa TensorMapping{Float64, 3, 4}
-    @test InflatedTensorMapping(I(3), C) isa TensorMapping{Float64, 2, 3}
-    @test InflatedTensorMapping(I(3), I(2,3)) isa TensorMapping{Float64, 3, 3}
-
-    @test range_size(InflatedTensorMapping(I(3,2), A, I(4))) == (3,2,4,4)
-    @test domain_size(InflatedTensorMapping(I(3,2), A, I(4))) == (3,2,2,4)
-
-    @test range_size(InflatedTensorMapping(I(3,2), B, I(4))) == (3,2,4,2,4)
-    @test domain_size(InflatedTensorMapping(I(3,2), B, I(4))) == (3,2,3,4)
-
-    @test range_size(InflatedTensorMapping(I(3), C, I(2,3))) == (3,4,2,3)
-    @test domain_size(InflatedTensorMapping(I(3), C, I(2,3))) == (3,2,3,2,3)
-
-    @inferred range_size(InflatedTensorMapping(I(3,2), A, I(4))) == (3,2,4,4)
-    @inferred domain_size(InflatedTensorMapping(I(3,2), A, I(4))) == (3,2,2,4)
+    @testset "Constructors" begin
+        @test InflatedTensorMapping(I(3,2), A, I(4)) isa TensorMapping{Float64, 4, 4}
+        @test InflatedTensorMapping(I(3,2), B, I(4)) isa TensorMapping{Float64, 5, 4}
+        @test InflatedTensorMapping(I(3), C, I(2,3)) isa TensorMapping{Float64, 4, 5}
+        @test InflatedTensorMapping(C, I(2,3)) isa TensorMapping{Float64, 3, 4}
+        @test InflatedTensorMapping(I(3), C) isa TensorMapping{Float64, 2, 3}
+        @test InflatedTensorMapping(I(3), I(2,3)) isa TensorMapping{Float64, 3, 3}
+    end
 
-    # Test InflatedTensorMapping mapping w. before and after
-    tm = InflatedTensorMapping(I(3,2), A, I(4))
-    v = rand(domain_size(tm)...)
-    @tullio IAIv[a,b,c,d] := Ã[c,i]*v[a,b,i,d]
-    @test tm*v ≈ IAIv rtol=1e-14
-    @inferred LazyTensors.split_index(tm,1,1,1,1)
+    @testset "Range and domain size" begin
+        @test range_size(InflatedTensorMapping(I(3,2), A, I(4))) == (3,2,4,4)
+        @test domain_size(InflatedTensorMapping(I(3,2), A, I(4))) == (3,2,2,4)
 
-    # Test InflatedTensorMapping mapping w. before
-    tm = InflatedTensorMapping(I(3,2), A)
-    v = rand(domain_size(tm)...)
-    @tullio IAIv[a,b,c] := Ã[c,i]*v[a,b,i]
-    @test tm*v ≈ IAIv rtol=1e-14
-    @inferred LazyTensors.split_index(tm,1,1,1)
+        @test range_size(InflatedTensorMapping(I(3,2), B, I(4))) == (3,2,4,2,4)
+        @test domain_size(InflatedTensorMapping(I(3,2), B, I(4))) == (3,2,3,4)
 
-    # Test InflatedTensorMapping mapping w. after
-    tm = InflatedTensorMapping(A,I(4))
-    v = rand(domain_size(tm)...)
-    @tullio IAIv[c,d] := Ã[c,i]*v[i,d]
-    @test tm*v ≈ IAIv rtol=1e-14
-    @inferred LazyTensors.split_index(tm,1,1)
+        @test range_size(InflatedTensorMapping(I(3), C, I(2,3))) == (3,4,2,3)
+        @test domain_size(InflatedTensorMapping(I(3), C, I(2,3))) == (3,2,3,2,3)
 
-    struct ScalingOperator{T,D} <: TensorMapping{T,D,D}
-        λ::T
-        size::NTuple{D,Int}
+        @inferred range_size(InflatedTensorMapping(I(3,2), A, I(4))) == (3,2,4,4)
+        @inferred domain_size(InflatedTensorMapping(I(3,2), A, I(4))) == (3,2,2,4)
     end
 
-    LazyTensors.apply(m::ScalingOperator{T,D}, v, I::Vararg{Index,D}) where {T,D} = m.λ*v[I]
-    LazyTensors.range_size(m::ScalingOperator) = m.size
-    LazyTensors.domain_size(m::ScalingOperator) = m.size
+    @testset "Application" begin
+        # Testing regular application and transposed application with inflation "before", "after" and "before and after".
+        # The inflated tensor mappings are chosen to preserve, reduce and increase the dimension of the result compared to the input.
+        tests = [
+            (
+                InflatedTensorMapping(I(3,2), A, I(4)),
+                (v-> @tullio res[a,b,c,d] := Ã[c,i]*v[a,b,i,d]), # Expected result of apply
+                (v-> @tullio res[a,b,c,d] := Ã[i,c]*v[a,b,i,d]), # Expected result of apply_transpose
+            ),
+            (
+                InflatedTensorMapping(I(3,2), B, I(4)),
+                (v-> @tullio res[a,b,c,d,e] := B̃[c,d,i]*v[a,b,i,e]),
+                (v-> @tullio res[a,b,c,d] := B̃[i,j,c]*v[a,b,i,j,d]),
+            ),
+            (
+                InflatedTensorMapping(I(3,2), C, I(4)),
+                (v-> @tullio res[a,b,c,d] := C̃[c,i,j]*v[a,b,i,j,d]),
+                (v-> @tullio res[a,b,c,d,e] := C̃[i,c,d]*v[a,b,i,e]),
+            ),
+            (
+                InflatedTensorMapping(I(3,2), A),
+                (v-> @tullio res[a,b,c] := Ã[c,i]*v[a,b,i]),
+                (v-> @tullio res[a,b,c] := Ã[i,c]*v[a,b,i]),
+            ),
+            (
+                InflatedTensorMapping(I(3,2), B),
+                (v-> @tullio res[a,b,c,d] := B̃[c,d,i]*v[a,b,i]),
+                (v-> @tullio res[a,b,c] := B̃[i,j,c]*v[a,b,i,j]),
+            ),
+            (
+                InflatedTensorMapping(I(3,2), C),
+                (v-> @tullio res[a,b,c] := C̃[c,i,j]*v[a,b,i,j]),
+                (v-> @tullio res[a,b,c,d] := C̃[i,c,d]*v[a,b,i]),
+            ),
+            (
+                InflatedTensorMapping(A,I(4)),
+                (v-> @tullio res[a,b] := Ã[a,i]*v[i,b]),
+                (v-> @tullio res[a,b] := Ã[i,a]*v[i,b]),
+            ),
+            (
+                InflatedTensorMapping(B,I(4)),
+                (v-> @tullio res[a,b,c] := B̃[a,b,i]*v[i,c]),
+                (v-> @tullio res[a,b] := B̃[i,j,a]*v[i,j,b]),
+            ),
+            (
+                InflatedTensorMapping(C,I(4)),
+                (v-> @tullio res[a,b] := C̃[a,i,j]*v[i,j,b]),
+                (v-> @tullio res[a,b,c] := C̃[i,a,b]*v[i,c]),
+            ),
+        ]
 
-    tm = InflatedTensorMapping(I(2,3),ScalingOperator(2.0, (3,2)),I(3,4))
-    v = rand(domain_size(tm)...)
+        @testset "apply" begin
+            for i ∈ 1:length(tests)
+                tm = tests[i][1]
+                v = rand(domain_size(tm)...)
+                true_value = tests[i][2](v)
+                @test tm*v ≈ true_value rtol=1e-14
+            end
+        end
+
+        @testset "apply_transpose" begin
+            for i ∈ 1:length(tests)
+                tm = tests[i][1]
+                v = rand(range_size(tm)...)
+                true_value = tests[i][3](v)
+                @test tm'*v ≈ true_value rtol=1e-14
+            end
+        end
 
-    @inferred LazyTensors.split_index(tm,1,2,3,2,2,4)
-    @inferred apply(tm,v,Index{Unknown}.((1,2,3,2,2,4))...)
-    @inferred (tm*v)[1,2,3,2,2,4]
+        @testset "Inference of application" begin
+            struct ScalingOperator{T,D} <: TensorMapping{T,D,D}
+                λ::T
+                size::NTuple{D,Int}
+            end
+
+            LazyTensors.apply(m::ScalingOperator{T,D}, v, I::Vararg{Any,D}) where {T,D} = m.λ*v[I...]
+            LazyTensors.range_size(m::ScalingOperator) = m.size
+            LazyTensors.domain_size(m::ScalingOperator) = m.size
+
+            tm = InflatedTensorMapping(I(2,3),ScalingOperator(2.0, (3,2)),I(3,4))
+            v = rand(domain_size(tm)...)
+
+            @inferred apply(tm,v,1,2,3,2,2,4)
+            @inferred (tm*v)[1,2,3,2,2,4]
+        end
+    end
 
     @testset "InflatedTensorMapping of InflatedTensorMapping" begin
         A = ScalingOperator(2.0,(2,3))
@@ -405,7 +452,20 @@
 
         @test InflatedTensorMapping(I(2), I(2), I(2)) isa InflatedTensorMapping # The constructor should always return its type.
     end
+end
 
+@testset "split_index" begin
+    @test LazyTensors.split_index(Val(2),Val(1),Val(2),Val(2),1,2,3,4,5,6) == ((1,2,:,5,6),(3,4))
+    @test LazyTensors.split_index(Val(2),Val(3),Val(2),Val(2),1,2,3,4,5,6) == ((1,2,:,:,:,5,6),(3,4))
+    @test LazyTensors.split_index(Val(3),Val(1),Val(1),Val(2),1,2,3,4,5,6) == ((1,2,3,:,5,6),(4,))
+    @test LazyTensors.split_index(Val(3),Val(2),Val(1),Val(2),1,2,3,4,5,6) == ((1,2,3,:,:,5,6),(4,))
+    @test LazyTensors.split_index(Val(1),Val(1),Val(2),Val(3),1,2,3,4,5,6) == ((1,:,4,5,6),(2,3))
+    @test LazyTensors.split_index(Val(1),Val(2),Val(2),Val(3),1,2,3,4,5,6) == ((1,:,:,4,5,6),(2,3))
+
+    @test LazyTensors.split_index(Val(0),Val(1),Val(3),Val(3),1,2,3,4,5,6) == ((:,4,5,6),(1,2,3))
+    @test LazyTensors.split_index(Val(3),Val(1),Val(3),Val(0),1,2,3,4,5,6) == ((1,2,3,:),(4,5,6))
+
+    @inferred LazyTensors.split_index(Val(2),Val(3),Val(2),Val(2),1,2,3,2,2,4)
 end
 
 @testset "slice_tuple" begin
@@ -415,6 +475,32 @@
     @test LazyTensors.slice_tuple((1,2,3,4,5,6),Val(4), Val(6)) == (4,5,6)
 end
 
+@testset "split_tuple" begin
+    @testset "2 parts" begin
+        @test LazyTensors.split_tuple((),Val(0)) == ((),())
+        @test LazyTensors.split_tuple((1,),Val(0)) == ((),(1,))
+        @test LazyTensors.split_tuple((1,),Val(1)) == ((1,),())
+
+        @test LazyTensors.split_tuple((1,2,3,4),Val(0)) == ((),(1,2,3,4))
+        @test LazyTensors.split_tuple((1,2,3,4),Val(1)) == ((1,),(2,3,4))
+        @test LazyTensors.split_tuple((1,2,3,4),Val(2)) == ((1,2),(3,4))
+        @test LazyTensors.split_tuple((1,2,3,4),Val(3)) == ((1,2,3),(4,))
+        @test LazyTensors.split_tuple((1,2,3,4),Val(4)) == ((1,2,3,4),())
+
+        @inferred LazyTensors.split_tuple((1,2,3,4),Val(3))
+    end
+
+    @testset "3 parts" begin
+        @test LazyTensors.split_tuple((),Val(0),Val(0)) == ((),(),())
+        @test LazyTensors.split_tuple((1,2,3),Val(1), Val(1)) == ((1,),(2,),(3,))
+
+        @test LazyTensors.split_tuple((1,2,3,4,5,6),Val(1),Val(2)) == ((1,),(2,3),(4,5,6))
+        @test LazyTensors.split_tuple((1,2,3,4,5,6),Val(3),Val(2)) == ((1,2,3),(4,5),(6,))
+
+        @inferred LazyTensors.split_tuple((1,2,3,4,5,6),Val(3),Val(2))
+    end
+end
+
 @testset "flatten_tuple" begin
     @test LazyTensors.flatten_tuple((1,)) == (1,)
     @test LazyTensors.flatten_tuple((1,2,3,4,5,6)) == (1,2,3,4,5,6)
@@ -430,7 +516,7 @@
         size::NTuple{D,Int}
     end
 
-    LazyTensors.apply(m::ScalingOperator{T,D}, v, I::Vararg{Index,D}) where {T,D} = m.λ*v[I]
+    LazyTensors.apply(m::ScalingOperator{T,D}, v, I::Vararg{Any,D}) where {T,D} = m.λ*v[I...]
     LazyTensors.range_size(m::ScalingOperator) = m.size
     LazyTensors.domain_size(m::ScalingOperator) = m.size
 
--- a/test/testSbpOperators.jl	Thu Nov 26 09:03:54 2020 +0100
+++ b/test/testSbpOperators.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -184,8 +184,8 @@
     v = evalOn(g,x->1+x^2)
     u = fill(3.124)
 
-    @test (e_l*v)[Index{Lower}(1)] == v[1]
-    @test (e_r*v)[Index{Upper}(4)] == v[end]
+    @test (e_l*v)[] == v[1]
+    @test (e_r*v)[] == v[end]
     @test e_l'*u == [u[], 0, 0, 0]
     @test e_r'*u == [0, 0, 0, u[]]
     @test_throws BoundsError (e_l*v)[Index{Lower}(3)]
@@ -218,19 +218,19 @@
     @test range_size(e_s) == (4,)
     @test range_size(e_n) == (4,)
 
-    I_w = [(Index{Lower}(1),Index{Lower}(1)),
-           (Index{Lower}(1),Index{Interior}(2)),
-           (Index{Lower}(1),Index{Interior}(3)),
-           (Index{Lower}(1),Index{Interior}(4)),
-           (Index{Lower}(1),Index{Upper}(5))]
+    I_w = [(Index{Lower}(1),),
+           (Index{Interior}(2),),
+           (Index{Interior}(3),),
+           (Index{Interior}(4),),
+           (Index{Upper}(5),)]
     v_w = [10,7,4,1.0,1];
     for i = 1:length(I_w)
-        @test_broken (e_w*v)[I_w[i]...] == v_w[i];
+        @test (e_w*v)[I_w[i]...] == v_w[i];
     end
-    @test_broken e_w*v == [10,7,4,1.0,1]
-    @test_broken e_e*v == [13,10,7,4,4.0]
-    @test_broken e_s*v == [10,11,12,13.0]
-    @test_broken e_n*v == [1,2,3,4.0]
+    @test e_w*v == [10,7,4,1.0,1]
+    @test e_e*v == [13,10,7,4,4.0]
+    @test e_s*v == [10,11,12,13.0]
+    @test e_n*v == [1,2,3,4.0]
 
     g_x = [1,2,3,4.0]
     g_y = [5,4,3,2,1.0]
@@ -247,10 +247,10 @@
     G_n = zeros(Float64, (4,5))
     G_n[:,5] = g_x
 
-    @test_broken e_w'*g_y == G_w
-    @test_broken e_e'*g_y == G_e
-    @test_broken e_s'*g_x == G_s
-    @test_broken e_n'*g_x == G_n
+    @test e_w'*g_y == G_w
+    @test e_e'*g_y == G_e
+    @test e_s'*g_x == G_s
+    @test e_n'*g_x == G_n
 end
 #
 # @testset "NormalDerivative" begin