--- a/src/LazyTensors/LazyTensors.jl	Thu Feb 02 22:53:11 2023 +0100
+++ b/src/LazyTensors/LazyTensors.jl	Fri Feb 03 22:14:47 2023 +0100
@@ -1,5 +1,7 @@
 module LazyTensors
 
+using Sbplib.RegionIndices
+
 export TensorApplication
 export TensorTranspose
 export TensorComposition

--- a/src/LazyTensors/lazy_tensor_operations.jl	Thu Feb 02 22:53:11 2023 +0100
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Fri Feb 03 22:14:47 2023 +0100
@@ -324,3 +324,30 @@
     print(io, "RangeSizeMismatch: ")
     print(io, "range size $(range_size(err.tm)) of LazyTensor not matching size $(err.sz)")
 end
+
+
+function apply_with_region(op, v, boundary_width::Integer, dim_size::Integer, i)
+    if 0 < i <= boundary_width
+        return LazyTensors.apply(op,v,Index(i,Lower))
+    elseif boundary_width < i <= dim_size-boundary_width
+        return LazyTensors.apply(op,v,Index(i,Interior))
+    elseif dim_size-boundary_width < i <= dim_size
+        return LazyTensors.apply(op,v,Index(i,Upper))
+    else
+        error("Bounds error") # TODO: Make this more standard
+    end
+end
+# TBD: Can these methods be merge by having a function as an arguement instead?
+# TODO: Add inference test that show where things break and how this rewrite fixes it.
+function apply_transpose_with_region(op, v, boundary_width::Integer, dim_size::Integer, i)
+    if 0 < i <= boundary_width
+        return LazyTensors.apply_transpose(op,v,Index(i,Lower))
+    elseif boundary_width < i <= dim_size-boundary_width
+        return LazyTensors.apply_transpose(op,v,Index(i,Interior))
+    elseif dim_size-boundary_width < i <= dim_size
+        return LazyTensors.apply_transpose(op,v,Index(i,Upper))
+    else
+        error("Bounds error") # TODO: Make this more standard
+    end
+end
+

--- a/src/RegionIndices/RegionIndices.jl	Thu Feb 02 22:53:11 2023 +0100
+++ b/src/RegionIndices/RegionIndices.jl	Fri Feb 03 22:14:47 2023 +0100
@@ -65,6 +65,7 @@
         error("Bounds error") # TODO: Make this more standard
     end
 end
+# 2022-02-21: Using the return values of getregion cause type inference to give up in ceratin cases for example H*H*v
 
 export getregion
 

--- a/src/SbpOperators/boundaryops/boundary_operator.jl	Thu Feb 02 22:53:11 2023 +0100
+++ b/src/SbpOperators/boundaryops/boundary_operator.jl	Fri Feb 03 22:14:47 2023 +0100
@@ -54,6 +54,5 @@
 end
 
 function LazyTensors.apply_transpose(op::BoundaryOperator, v::AbstractArray{<:Any,0}, i)
-    r = getregion(i, closure_size(op), op.size)
-    apply_transpose(op, v, Index(i,r))
+    return LazyTensors.apply_transpose_with_region(op, v, closure_size(op), op.size[1], i)
 end

--- a/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl	Thu Feb 02 22:53:11 2023 +0100
+++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl	Fri Feb 03 22:14:47 2023 +0100
@@ -41,8 +41,7 @@
 end
 
 function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i)
-    r = getregion(i, closure_size(op), op.size[1])
-    return LazyTensors.apply(op, v, Index(i, r))
+    return LazyTensors.apply_with_region(op, v, closure_size(op), op.size[1], i)
 end
 
 LazyTensors.apply_transpose(op::ConstantInteriorScalingOperator, v, i) = apply(op, v, i)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/inference_trouble.txt	Fri Feb 03 22:14:47 2023 +0100
@@ -0,0 +1,47 @@
+Innan ändringarna på den här branchen:
+
+begin
+    using Sbplib
+    using Sbplib.Grids
+    using Sbplib.SbpOperators
+
+    g = EquidistantGrid((10,10),(0.,0.), (1.,1.))
+    v = evalOn(g, (x,y)->x^2+y^2+1)
+    H = inner_product(g, 1., [1/2])
+end
+
+# Not type stable
+LazyTensors.apply(H.t1, H.t2*v, 1,2) 
+@code_warntype LazyTensors.apply(H.t1, H.t2*v, 1,2)
+@code_warntype LazyTensors.apply(H.t1.tm, view(H.t2*v,:,1), 2)
+
+# Nedan är halvdåliga
+@code_warntype LazyTensors.apply(H.t1.tm, view(H.t2*v,1,:), 2)
+@code_warntype LazyTensors.apply(H.t1.tm, view(v,1,:), 2)
+@code_warntype LazyTensors.apply(H.t1.tm, view(v,:,1), 2)
+@code_warntype LazyTensors.apply(H.t1.tm, v[:,1], 2)
+
+
+
+
+
+
+
+
+
+begin
+    using Sbplib
+    using Sbplib.Grids
+    using Sbplib.SbpOperators
+    import Sbplib.SbpOperators: Stencil
+    using Sbplib.RegionIndices
+
+    g = EquidistantGrid(10,0., 1.)
+    v = evalOn(g, (x)->x^2+1)
+    H = inner_product(g, 1., [1/2])
+    V = SbpOperators.VolumeOperator(g, Stencil(1.,center=1), (Stencil(1/2,center=1),), SbpOperators.even)
+    b = SbpOperators.BoundaryOperator(g, Stencil(1/2,center=1), Lower())
+end
+
+@code_warntype LazyTensors.apply(H, H*v, 2)
+@code_warntype LazyTensors.apply(V, V*v, 2)

--- a/src/SbpOperators/volumeops/volume_operator.jl	Thu Feb 02 22:53:11 2023 +0100
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Fri Feb 03 22:14:47 2023 +0100
@@ -32,7 +32,6 @@
 end
 
 function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i)
-    r = getregion(i, closure_size(op), op.size[1])
-    return LazyTensors.apply(op, v, Index(i, r))
+    return LazyTensors.apply_with_region(op, v, closure_size(op), op.size[1], i)
 end
 # TODO: Move this to LazyTensors when we have the region communication down.