--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/wave_eq.jl	Thu Jan 20 21:51:53 2022 +0100
@@ -0,0 +1,98 @@
+using Sbplib.Grids, Sbplib.SbpOperators, Sbplib.LazyTensors, Sbplib.RegionIndices
+using OrdinaryDiffEq, Plots, Printf, Base.Threads
+
+function apply_tm!(f,u,tm,ind)
+    for I in ind
+        @inbounds f[I] = (tm*u)[I]
+    end
+end
+
+function apply_tm_all_regions!(f,u,tm,rinds)
+    apply_tm!(f,u,tm,rinds[1])
+    apply_tm!(f,u,tm,rinds[2])
+    apply_tm!(f,u,tm,rinds[3])
+end
+
+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 wave_eq_sim(alg,T,CFL)
+    # Domain
+    N = 101
+    g = EquidistantGrid(N,0.,1.)
+    dx = min(spacing(g)...)
+
+    # Spatial discretization
+    Δ = Laplace(g,sbp_operators_path()*"standard_diagonal.toml"; order=4)
+    (id_l, id_r) = boundary_identifiers(g)
+    SAT_l = boundary_condition(Δ,id_l,NeumannBC())
+    SAT_r = boundary_condition(Δ,id_r,NeumannBC())
+    tm = (Δ + SAT_l + SAT_r)
+
+    # RHS function
+    rinds = get_region_indices(Δ,N)
+    function f(du,u,p,t)
+        du[1:N] .= u[N+1:end]
+        apply_tm_all_regions!(view(du,N+1:2*N), view(u,1:N), tm, rinds)
+    end
+    # Initial condition
+    sigma = 0.1
+    ic_u = x->1/(sigma*sqrt(2*pi))*exp(-1/2*((x-0.5)^2/sigma^2))
+    ic_u_t = x->0
+    w0 = [evalOn(g,ic_u);
+          evalOn(g,ic_u_t)]
+    # Setup ODE and solve
+    tspan = (0.,T)
+    prob = ODEProblem(f,w0,tspan)
+    sol = solve(prob, alg, dt=CFL*dx, saveat=0.05)
+
+    # Plotting
+    x = [x[1] for x in points(g)]
+    anim = @animate for i ∈ eachindex(sol.t)
+        u_i = sol.u[i]
+        plot(x, u_i[1:N], ylims = (0,4), lw=3,ls=:dash,label="",title=@sprintf("u at t = %.3f", sol.t[i]))
+    end
+    gif(anim, "wave.gif", fps = 15)
+end
+
+wave_eq_sim(CarpenterKennedy2N54(),1.,0.25)
+
+# function boundary_condition(L,id; )
+#     neumann_closure
+#     neumann_penalty
+# end
+#
+# function neumann_bc(L,id)
+#     e = boundary_restriction.(L,ids)
+#     d = normal_derivative(L,ids)
+#     return (e'∘d,)
+# end
+#
+# function (closure, penalty) neumann_bc(L,id,g::Function)
+#     e = boundary_restriction.(L,ids)
+#     d = normal_derivative.(L,ids)
+#     return e'∘d
+# end
+#
+# function dirichlet_closure(L,id)
+#     e = boundary_restriction.(L,ids)
+#     d = normal_derivative.(L,ids)
+#     return ...
+# end
+#
+# function SAT(L,ids,BCs)
+#     BC = BCs(L,ids[1])
+#     for i = 2:length(ids)
+#         BC = BC+ BCs(L,ids[1])
+#     end
+#     H_inv = inverse_inner_product(L)
+#     return H_inv∘BC
+# end