--- a/src/SbpOperators/volumeops/volume_operator.jl	Fri Jul 02 14:23:33 2021 +0200
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Wed Jan 19 14:44:24 2022 +0100
@@ -1,21 +1,22 @@
 """
-    volume_operator(grid,inner_stencil,closure_stencils,parity,direction)
+    volume_operator(grid, inner_stencil, closure_stencils, parity, direction)
+
 Creates a volume operator on a `Dim`-dimensional grid acting along the
-specified coordinate `direction`. The action of the operator is determined by the
-stencils `inner_stencil` and `closure_stencils`.
-When `Dim=1`, the corresponding `VolumeOperator` tensor mapping is returned.
-When `Dim>1`, the `VolumeOperator` `op` is inflated by the outer product
-of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`,
-the boundary restriction operator in the y-direction direction is `Ix⊗op⊗Iz`.
+specified coordinate `direction`. The action of the operator is determined by
+the stencils `inner_stencil` and `closure_stencils`. When `Dim=1`, the
+corresponding `VolumeOperator` tensor mapping is returned. When `Dim>1`, the
+returned operator is the appropriate outer product of a one-dimensional
+operators and `IdentityMapping`s, e.g for `Dim=3` the volume operator in the
+y-direction is `I⊗op⊗I`.
 """
-function volume_operator(grid::EquidistantGrid{Dim,T}, inner_stencil::Stencil{T}, closure_stencils::NTuple{M,Stencil{T}}, parity, direction) where {Dim,T,M}
+function volume_operator(grid::EquidistantGrid, inner_stencil, closure_stencils, parity, direction)
     #TODO: Check that direction <= Dim?
 
     # Create 1D volume operator in along coordinate direction
     op = VolumeOperator(restrict(grid, direction), inner_stencil, closure_stencils, parity)
     # Create 1D IdentityMappings for each coordinate direction
-    one_d_grids = restrict.(Ref(grid), Tuple(1:Dim))
-    Is = IdentityMapping{T}.(size.(one_d_grids))
+    one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid)))
+    Is = IdentityMapping{eltype(grid)}.(size.(one_d_grids))
     # Formulate the correct outer product sequence of the identity mappings and
     # the volume operator
     parts = Base.setindex(Is, op, direction)
@@ -34,7 +35,7 @@
 end
 
 function VolumeOperator(grid::EquidistantGrid{1}, inner_stencil, closure_stencils, parity)
-    return VolumeOperator(inner_stencil, closure_stencils, size(grid), parity)
+    return VolumeOperator(inner_stencil, Tuple(closure_stencils), size(grid), parity)
 end
 
 closure_size(::VolumeOperator{T,N,M}) where {T,N,M} = M