--- a/src/LazyTensors/lazy_tensor_operations.jl	Sun May 08 11:35:22 2022 +0200
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Tue May 10 20:24:20 2022 +0200
@@ -269,6 +269,29 @@
 LazyOuterProduct(tms::Vararg{LazyTensor}) = foldl(LazyOuterProduct, tms)
 
 
+
+"""
+    inflate(tm::LazyTensor, sz, dir)
+
+Inflate `tm` such that it gets the size `sz` in all directions except `dir`.
+Here `sz[dir]` is ignored and replaced with the range and domains size of
+`tm`.
+
+An example of when this operation is useful is when extending a one
+dimensional difference operator `D` to a 2D grid of a ceratin size. In that
+case we could have
+
+```julia
+Dx = inflate(D, (10,10), 1)
+Dy = inflate(D, (10,10), 2)
+```
+"""
+function inflate(tm::LazyTensor, sz, dir)
+    Is = IdentityTensor{eltype(tm)}.(sz)
+    parts = Base.setindex(Is, tm, dir)
+    return foldl(⊗, parts)
+end
+
 function check_domain_size(tm::LazyTensor, sz)
     if domain_size(tm) != sz
         throw(DomainSizeMismatch(tm,sz))

--- a/src/SbpOperators/boundaryops/boundary_operator.jl	Sun May 08 11:35:22 2022 +0200
+++ b/src/SbpOperators/boundaryops/boundary_operator.jl	Tue May 10 20:24:20 2022 +0200
@@ -12,20 +12,16 @@
 function boundary_operator(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary)
     #TODO:Check that dim(boundary) <= Dim?
 
-    # Create 1D boundary operator
-    r = region(boundary)
     d = dim(boundary)
-    op = BoundaryOperator(restrict(grid, d), closure_stencil, r)
+    op = BoundaryOperator(restrict(grid, d), closure_stencil, region(boundary))
 
     # Create 1D IdentityTensors for each coordinate direction
     one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid)))
     Is = IdentityTensor{eltype(grid)}.(size.(one_d_grids))
 
-    # Formulate the correct outer product sequence of the identity mappings and
-    # the boundary operator
-    parts = Base.setindex(Is, op, d)
-    return foldl(⊗, parts)
+    return LazyTensors.inflate(op, size(grid), d)
 end
+# TBD: Should the inflation happen here or should we remove this method and do it at the caller instead?
 
 """
     BoundaryOperator{T,R,N} <: LazyTensor{T,0,1}

--- a/src/SbpOperators/volumeops/volume_operator.jl	Sun May 08 11:35:22 2022 +0200
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Tue May 10 20:24:20 2022 +0200
@@ -12,16 +12,10 @@
 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 IdentityTensors for each coordinate direction
-    one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid)))
-    Is = IdentityTensor{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)
-    return foldl(⊗, parts)
+    return LazyTensors.inflate(op, size(grid), direction)
 end
+# TBD: Should the inflation happen here or should we remove this method and do it at the caller instead?
 
 """
     VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1}
@@ -59,3 +53,4 @@
     r = getregion(i, closure_size(op), op.size[1])
     return LazyTensors.apply(op, v, Index(i, r))
 end
+# TODO: Move this to LazyTensors when we have the region communication down.

--- a/test/LazyTensors/lazy_tensor_operations_test.jl	Sun May 08 11:35:22 2022 +0200
+++ b/test/LazyTensors/lazy_tensor_operations_test.jl	Tue May 10 20:24:20 2022 +0200
@@ -366,3 +366,17 @@
         @test I1⊗Ã⊗I2 == InflatedTensor(I1, Ã, I2)
     end
 end
+
+@testset "inflate" begin
+    I = LazyTensors.inflate(IdentityTensor(),(3,4,5,6), 2)
+    @test I isa LazyTensor{Float64, 3,3}
+    @test range_size(I) == (3,5,6)
+    @test domain_size(I) == (3,5,6)
+
+    @test LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 1) == InflatedTensor(IdentityTensor{Float64}(),ScalingTensor(1., (4,)),IdentityTensor(4,5,6))
+    @test LazyTensors.inflate(ScalingTensor(2., (1,)),(3,4,5,6), 2) == InflatedTensor(IdentityTensor(3),ScalingTensor(2., (1,)),IdentityTensor(5,6))
+    @test LazyTensors.inflate(ScalingTensor(3., (6,)),(3,4,5,6), 4) == InflatedTensor(IdentityTensor(3,4,5),ScalingTensor(3., (6,)),IdentityTensor{Float64}())
+
+    @test_throws BoundsError LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 0)
+    @test_throws BoundsError LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 5)
+end