--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/laplace_benchmark.jl	Wed Jan 19 13:15:45 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)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/laplace_benchmark.m	Wed Jan 19 13:15:45 2022 +0100
@@ -0,0 +1,14 @@
+m = 4001;
+ops = sbp.D2Standard(m,{0,1},4);
+D2 = ops.D2;
+u = linspace(0,1,m)';
+f = zeros(size(u));
+
+nsample = 10000;
+ts = zeros(nsample,1);
+
+for i = 1:nsample
+    tic; f = D2*u; t = toc;
+    ts(i) = t;
+end
+min(ts)
\ No newline at end of file

--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Fri Jul 02 14:23:33 2021 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Wed Jan 19 13:15:45 2022 +0100
@@ -22,6 +22,7 @@
     e::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary restriction operators.
     d::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Normal derivative operators
     H_boundary::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary quadrature operators # TODO: Boundary inner product?
+    order::Int
 end
 export Laplace
 
@@ -48,7 +49,7 @@
     d_pairs = ntuple(i -> ids[i] => normal_derivative(grid,d_closure_stencil,ids[i]),n_ids)
     Hᵧ_pairs = ntuple(i -> ids[i] => inner_product(boundary_grid(grid,ids[i]),H_closure_stencils),n_ids)
 
-    return Laplace(Δ, H, H⁻¹, StaticDict(e_pairs), StaticDict(d_pairs), StaticDict(Hᵧ_pairs))
+    return Laplace(Δ, H, H⁻¹, StaticDict(e_pairs), StaticDict(d_pairs), StaticDict(Hᵧ_pairs), order)
 end
 
 # TODO: Consider pretty printing of the following form
@@ -58,6 +59,14 @@
 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D)
 LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...)
 
+"""
+    closure_size(L::Lapalace)
+
+Returns the inner product operator associated with `L`
+
+"""
+closure_size(L::Laplace) = closure_size(L.D)
+export closure_size
 
 """
     inner_product(L::Lapalace)
@@ -124,6 +133,54 @@
 boundary_quadrature(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.H_boundary[ids[i]],N)
 export boundary_quadrature
 
+abstract type BoundaryConditionType end
+struct NeumannBC <: BoundaryConditionType end
+struct DirichletBC <: BoundaryConditionType end
+export NeumannBC
+
+boundary_condition(L, id, ::NeumannBC) = neumann_bc(L, id)
+boundary_condition(L, id, ::DirichletBC) = dirichlet_bc(L, id)
+export boundary_condition
+
+function neumann_bc(L::Laplace, id::BoundaryIdentifier)
+    H_inv = inverse_inner_product(L)
+    e = boundary_restriction(L,id)
+    d = normal_derivative(L,id)
+    H_b = boundary_quadrature(L,id)
+    tau = H_inv∘e'∘H_b
+    closure = tau∘d
+    # TODO: Return penalty once we have implemented scalar scaling of the operators.
+    return closure
+end
+
+function dirichlet_bc(L::Laplace, id::BoundaryIdentifier)
+    error("Not implemented")
+    # H_inv = inverse_inner_product(L)
+    # e = boundary_restriction(L,id)
+    # d = normal_derivative(L,id)
+    # H_b = boundary_quadrature(L,id)
+    # gamma = borrowing_parameter(L)
+    # tuning = 1.2
+    # S = ScalingOperator(tuning * -1.0/gamma)
+    # tau = H_inv∘(S∘e' + d')∘H_b
+    # closure = tau∘e
+    # penalty = ScalingOperator(-1)∘tau
+    # return (closure, penalty)
+end
+
+# function borrowing_parameter(L)
+#     if L.order == 2
+#         return 0.4
+#     elseif L.order == 4
+#         return 0.2508
+#     elseif L.order == 6
+#         return 0.1878
+#     elseif L.order == 8
+#         return 0.0015
+#     elseif L.order == 10
+#         return 0.0351
+#     end
+# end
 
 """
     laplace(grid, inner_stencil, closure_stencils)

Binary file wave.gif has changed