--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/benchmarks/laplace_benchmark.jl	Thu Jan 27 11:00:31 2022 +0100
@@ -0,0 +1,42 @@
+using Sbplib.Grids, Sbplib.SbpOperators, Sbplib.LazyTensors, Sbplib.RegionIndices
+using Profile, BenchmarkTools
+
+function apply_laplace!(f, u, L, inds)
+    for I in inds
+        f[I] = (L*u)[I]
+    end
+end
+
+apply_laplace!(f, u, L) = apply_laplace!(f, u, L, eachindex(u))
+
+region_indices(L, N, ::Lower) = map(x->Index{Lower}(x),1:closure_size(L))
+region_indices(L, N, ::Interior) = map(x->Index{Interior}(x),closure_size(L)+1:N-closure_size(L))
+region_indices(L, N, ::Upper) = map(x->Index{Upper}(x),N-closure_size(L)+1:N)
+
+function get_region_indices(L,N)
+    ind_lower = region_indices(L, N, Lower())
+    ind_interior = region_indices(L, N, Interior())
+    ind_upper = region_indices(L, N, Upper())
+    return (ind_lower, ind_interior, ind_upper)
+end
+
+function apply_laplace_regions!(f, u, L, region_inds)
+    apply_laplace!(f, u, L, region_inds[1])
+    apply_laplace!(f, u, L, region_inds[2])
+    apply_laplace!(f, u, L, region_inds[3])
+end
+
+# Domain
+N = 4001
+g = EquidistantGrid(N,0.,1.)
+
+# Operators
+L = Laplace(g,sbp_operators_path()*"standard_diagonal.toml"; order=4)
+u = evalOn(g,x->x^2)
+f = similar(u)
+
+apply_laplace!(f,u,L) #ensure compilation
+@btime apply_laplace!(f,u,L)
+rinds = get_region_indices(L,N)
+apply_laplace_regions!(f,u,L,rinds) #ensure compilation
+@btime apply_laplace_regions!(f,u,L,rinds)