--- a/Notes.md	Sat May 25 16:02:55 2024 -0700
+++ b/Notes.md	Tue Jun 25 15:32:19 2024 +0200
@@ -1,5 +1,26 @@
 # Notes
 
+## How to dispatch for different operators
+We have a problem in how dispatch for different operators work.
+ * We want to keep the types simple and flat (Awkward to forward `apply`)
+ * We want to dispatch SATs on the parameters of the continuous operator. (a * div for example)
+ * We want to allow keeping the same stencil_set across different calls. (maybe not so bad for the user to be responsible)
+
+Could remove the current opset idea and introduce a description of continuous operators
+ ```julia
+abstract type DifferentialOperator end
+
+struct Laplace <: DifferentialOperator end
+struct Advection <: DifferentialOperator
+    v
+end
+
+difference_operator(::Laplace, grid, stencil_set) = ... # Returns a plain LazyTensor. Replaces the current `laplace()` function.
+sat_tensors(::Laplace, grid, stencil_set, bc) = ...
+
+sat(::DifferentialOperator, grid, stencil_set, bc) = ...
+ ```
+
 ## Reading operators
 
 Jonatan's suggestion is to add methods to `Laplace`, `SecondDerivative` and

--- a/Project.toml	Sat May 25 16:02:55 2024 -0700
+++ b/Project.toml	Tue Jun 25 15:32:19 2024 +0200
@@ -1,7 +1,7 @@
 name = "Sbplib"
 uuid = "5a373a26-915f-4769-bcab-bf03835de17b"
 authors = ["Jonatan Werpers <jonatan@werpers.com>", "Vidar Stiernström <vidar.stiernstrom@it.uu.se>, and contributors"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

--- a/TODO.md	Sat May 25 16:02:55 2024 -0700
+++ b/TODO.md	Tue Jun 25 15:32:19 2024 +0200
@@ -21,6 +21,10 @@
        See (https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing)
  - [ ] Samla noggrannhets- och SBP-ness-tester för alla operatorer på ett ställe
  - [ ] Move export statements to top of each module
+ - [ ] Implement apply_transpose for
+      - [ ] ElementwiseTensorOperation
+      - [ ] VolumeOperator
+      - [ ] Laplace
 
 
  - [ ] Gå igenom alla typ parametrar och kolla om de är motiverade. Både i signaturer och typer, tex D i VariableSecondDerivative. Kan vi använda promote istället?

--- a/benchmark/benchmarks.jl	Sat May 25 16:02:55 2024 -0700
+++ b/benchmark/benchmarks.jl	Tue Jun 25 15:32:19 2024 +0200
@@ -15,7 +15,7 @@
 ll(d) = ntuple(i->0., d)
 lu(d) = ntuple(i->1., d)
 
-g1 = equidistant_grid(ll(1)[1], lu(1)[1], sz(1)[1])
+g1 = equidistant_grid(ll(1)[1], lu(1)[1], sz(1)...)
 g2 = equidistant_grid(ll(2), lu(2), sz(2)...)
 g3 = equidistant_grid(ll(3), lu(3), sz(3)...)
 

--- a/src/Grids/grid.jl	Sat May 25 16:02:55 2024 -0700
+++ b/src/Grids/grid.jl	Tue Jun 25 15:32:19 2024 +0200
@@ -74,7 +74,7 @@
     end
 end
 
-Base.size(cv) = size(cv.v)
+Base.size(cv::ArrayComponentView) = size(cv.v)
 Base.getindex(cv::ArrayComponentView, i::Int) = cv.v[i][cv.component_index...]
 Base.getindex(cv::ArrayComponentView, I::Vararg{Int}) = cv.v[I...][cv.component_index...]
 IndexStyle(::Type{<:ArrayComponentView{<:Any,<:Any,AT}}) where AT = IndexStyle(AT)

--- a/src/LazyTensors/lazy_tensor_operations.jl	Sat May 25 16:02:55 2024 -0700
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Tue Jun 25 15:32:19 2024 +0200
@@ -121,6 +121,7 @@
 
 Base.:*(a::T, tm::LazyTensor{T}) where T = TensorComposition(ScalingTensor{T,range_dim(tm)}(a,range_size(tm)), tm)
 Base.:*(tm::LazyTensor{T}, a::T) where T = a*tm
+Base.:-(tm::LazyTensor) = (-one(eltype(tm)))*tm
 
 """
     InflatedTensor{T,R,D} <: LazyTensor{T,R,D}

--- a/src/SbpOperators/SbpOperators.jl	Sat May 25 16:02:55 2024 -0700
+++ b/src/SbpOperators/SbpOperators.jl	Tue Jun 25 15:32:19 2024 +0200
@@ -11,7 +11,6 @@
 export sbp_operators_path
 
 # Operators
-export boundary_quadrature
 export boundary_restriction
 export inner_product
 export inverse_inner_product
@@ -22,20 +21,34 @@
 export second_derivative
 export second_derivative_variable
 export undivided_skewed04
-
-using Sbplib.RegionIndices
-using Sbplib.LazyTensors
-using Sbplib.Grids
+export closure_size
 
 @enum Parity begin
     odd = -1
     even = 1
 end
 
-export closure_size
 
+# Boundary conditions
+export BoundaryCondition
+export NeumannCondition
+export DirichletCondition
+export discretize_data
+export boundary_data
+export boundary
+export sat
+export sat_tensors
+
+# Using
+using Sbplib.RegionIndices
+using Sbplib.LazyTensors
+using Sbplib.Grids
+
+# Includes
 include("stencil.jl")
 include("stencil_set.jl")
+include("boundary_conditions/boundary_condition.jl")
+include("boundary_conditions/sat.jl")
 include("volumeops/volume_operator.jl")
 include("volumeops/stencil_operator_distinct_closures.jl")
 include("volumeops/constant_interior_scaling_operator.jl")

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/boundary_conditions/boundary_condition.jl	Tue Jun 25 15:32:19 2024 +0200
@@ -0,0 +1,62 @@
+"""
+    BoundaryCondition
+
+Description of a boundary condition. Implementations describe the kind of
+boundary condition, what boundary the condition applies to, and any associated
+data. Should implement [`boundary`](@ref) and may implement
+[`boundary_data`](@ref) if applicable.
+
+For examples see [`DirichletCondition`](@ref) and [`NeumannCondition`](@ref)
+"""
+abstract type BoundaryCondition end
+
+"""
+    boundary(::BoundaryCondition)
+
+The boundary identifier of the BoundaryCondition.
+"""
+function boundary end
+
+"""
+    boundary_data(::BoundaryCondition)
+
+If implemented, the data associated with the BoundaryCondition.
+"""
+function boundary_data end
+
+"""
+    discretize_data(grid, bc::BoundaryCondition)
+
+The data of `bc` as a lazily evaluated grid function on the boundary grid
+specified by `boundary(bc)`.
+"""
+function discretize_data(grid, bc::BoundaryCondition)
+    return eval_on(boundary_grid(grid, boundary(bc)), boundary_data(bc))
+end
+
+"""
+    DirichletCondition{DT,BID}
+
+A Dirichlet condition with `data::DT` on the boundary
+specified by the boundary identifier `BID`.
+"""
+struct DirichletCondition{DT,BID} <: BoundaryCondition
+    data::DT
+    boundary::BID
+end
+boundary_data(bc::DirichletCondition) = bc.data
+boundary(bc::DirichletCondition) = bc.boundary
+
+"""
+    NeumannCondition{DT,BID}
+
+A Neumann condition with `data::DT` on the boundary
+specified by the boundary identifier `BID`.
+"""
+struct NeumannCondition{DT,BID} <: BoundaryCondition
+    data::DT
+    boundary::BID
+end
+boundary_data(bc::NeumannCondition) = bc.data
+boundary(bc::NeumannCondition) = bc.boundary
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/boundary_conditions/sat.jl	Tue Jun 25 15:32:19 2024 +0200
@@ -0,0 +1,25 @@
+"""
+    sat_tensors(op, grid, bc::BoundaryCondition; kwargs...)
+
+The penalty tensor and boundary operator used to construct a
+simultaneous-approximation-term for imposing `bc` related to `op`.
+
+For `penalty_tensor, L  = sat_tensors(...)` then `SAT(u,g) =
+penalty_tensor*(L*u - g)`  where `g` is the boundary data.
+"""
+function sat_tensors end
+
+
+"""
+    sat(op, grid, bc::BoundaryCondition; kwargs...)
+
+Simultaneous-Approximation-Term for a general `bc` to `op`. Returns a function
+`SAT(u,g)` weakly imposing `bc` when added to `op*u`.
+
+Internally `sat_tensors(op, grid, bc; ...)` is called to construct the
+necessary parts for the SAT.
+"""
+function sat(op, grid, bc::BoundaryCondition; kwargs...)
+    penalty_tensor, L = sat_tensors(op, grid, bc; kwargs...)
+    return SAT(u, g) = penalty_tensor*(L*u - g)
+end

--- a/src/SbpOperators/operators/standard_diagonal.toml	Sat May 25 16:02:55 2024 -0700
+++ b/src/SbpOperators/operators/standard_diagonal.toml	Tue Jun 25 15:32:19 2024 +0200
@@ -34,11 +34,14 @@
     {s = ["1", "-2", "1"], c = 1},
 ]
 
+D2.positivity = {theta_M = "0.3636363636", theta_R = "1.000000538455350", m_b = "2"}
+
 D2variable.inner_stencil = [["1/2", "1/2", "0"],[ "-1/2", "-1", "-1/2"],["0", "1/2", "1/2"]]
 D2variable.closure_stencils = [
         {s = [["2", "-1", "0"],["-3", "1",   "0"],["1","0","0"]], c = 1},
 ]
 
+
 [[stencil_set]]
 
 order = 4
@@ -65,6 +68,8 @@
     {s = [ "-1/49",     "0",   "59/49", "-118/49", "64/49", "-4/49"], c = 4},
 ]
 
+D2.positivity = {theta_M = "0.2505765857", theta_R = "0.577587500088313", m_b = "4"}
+
 D2variable.inner_stencil = [
     ["-1/8",   "1/6", "-1/8",   "0",    "0"  ],
     [ "1/6",   "1/2",  "1/2",  "1/6",   "0"  ],
@@ -136,7 +141,6 @@
 ]
 
 
-
 [[stencil_set]]
 
 order = 6
@@ -150,3 +154,5 @@
 
 e.closure = ["1"]
 d1.closure = ["-25/12", "4", "-3", "4/3", "-1/4"]
+
+D2.positivity = {theta_M = "0.1878687080", theta_R = "0.3697", m_b = "7"}

--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Sat May 25 16:02:55 2024 -0700
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Tue Jun 25 15:32:19 2024 +0200
@@ -52,3 +52,78 @@
     return Δ
 end
 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
+
+"""
+    sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
+
+The operators required to construct the SAT for imposing a Dirichlet
+condition. `H_tuning` and `R_tuning` are used to specify the strength of the
+penalty.
+
+See also: [`sat`](@ref),[`DirichletCondition`](@ref), [`positivity_decomposition`](@ref).
+"""
+function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.)
+    id = boundary(bc)
+    set  = Δ.stencil_set
+    H⁻¹ = inverse_inner_product(g,set)
+    Hᵧ = inner_product(boundary_grid(g, id), set)
+    e = boundary_restriction(g, set, id)
+    d = normal_derivative(g, set, id)
+    B = positivity_decomposition(Δ, g, bc; H_tuning, R_tuning)
+    penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ
+    return penalty_tensor, e
+end
+
+"""
+    sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
+
+The operators required to construct the SAT for imposing a Neumann condition.
+
+See also: [`sat`](@ref), [`NeumannCondition`](@ref).
+"""
+function sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
+    id = boundary(bc)
+    set  = Δ.stencil_set
+    H⁻¹ = inverse_inner_product(g,set)
+    Hᵧ = inner_product(boundary_grid(g, id), set)
+    e = boundary_restriction(g, set, id)
+    d = normal_derivative(g, set, id)
+
+    penalty_tensor = -H⁻¹∘e'∘Hᵧ
+    return penalty_tensor, d
+end
+
+"""
+    positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
+
+Constructs the scalar `B` such that `d' - 1/2*B*e'` is symmetric positive
+definite with respect to the boundary quadrature. Here `d` is the normal
+derivative and `e` is the boundary restriction operator. `B` can then be used
+to form a symmetric and energy stable penalty for a Dirichlet condition. The
+parameters `H_tuning` and `R_tuning` are used to specify the strength of the
+penalty and must be greater than 1. For details we refer to
+https://doi.org/10.1016/j.jcp.2020.109294
+"""
+function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
+    @assert(H_tuning ≥ 1.)
+    @assert(R_tuning ≥ 1.)
+    Nτ_H, τ_R = positivity_limits(Δ,g,bc)
+    return H_tuning*Nτ_H + R_tuning*τ_R
+end
+
+# TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then
+# change bc::BoundaryCondition to id::BoundaryIdentifier
+function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition)
+    h = spacing(g)
+    θ_H = parse_scalar(Δ.stencil_set["H"]["closure"][1])
+    θ_R = parse_scalar(Δ.stencil_set["D2"]["positivity"]["theta_R"])
+
+    τ_H = 1/(h*θ_H)
+    τ_R = 1/(h*θ_R)
+    return τ_H, τ_R
+end
+
+function positivity_limits(Δ::Laplace, g::TensorGrid, bc::DirichletCondition)
+    τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(boundary(bc))], bc)
+    return τ_H*ndims(g), τ_R
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/boundary_conditions/boundary_condition_test.jl	Tue Jun 25 15:32:19 2024 +0200
@@ -0,0 +1,39 @@
+using Test
+
+using Sbplib.Grids
+using Sbplib.RegionIndices
+using Sbplib.SbpOperators
+
+@testset "BoundaryCondition" begin
+    grid_1d = equidistant_grid(0.0, 1.0, 11)
+    grid_2d = equidistant_grid((0.0, 0.0), (1.0,1.0), 11, 15)
+    grid_3d = equidistant_grid((0.0, 0.0, 0.0), (1.0,1.0, 1.0), 11, 15, 13)
+    (id_l,_) = boundary_identifiers(grid_1d)
+    (_,_,_,id_n) = boundary_identifiers(grid_2d)
+    (_,_,_,_,id_b,_) = boundary_identifiers(grid_3d)
+
+    g = 3.14
+    f(x,y,z) = x^2+y^2+z^2
+    @testset "Constructors" begin
+        @test DirichletCondition(g,id_l) isa DirichletCondition{Float64,Lower}
+        @test NeumannCondition(f,id_b) isa NeumannCondition{<:Function,CartesianBoundary{3,Lower}}
+    end
+
+    @testset "boundary" begin
+        @test boundary(DirichletCondition(g,id_l)) == id_l
+        @test boundary(NeumannCondition(f,id_b)) == id_b
+    end
+
+    @testset "boundary_data" begin
+        @test boundary_data(DirichletCondition(g,id_l)) == g
+        @test boundary_data(NeumannCondition(f,id_b)) == f
+    end
+
+    @testset "discretize_data" begin
+        @test fill(g) ≈ discretize_data(grid_1d,DirichletCondition(g,id_l))
+        @test g*ones(11,1) ≈ discretize_data(grid_2d,DirichletCondition(g,id_n))
+        X = repeat(0:1/10:1, inner = (1,15))
+        Y = repeat(0:1/14:1, outer = (1,11))
+        @test map((x,y)->f(x,y,0), X,Y') ≈ discretize_data(grid_3d,NeumannCondition(f,id_b))
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/boundary_conditions/sat_test.jl	Tue Jun 25 15:32:19 2024 +0200
@@ -0,0 +1,71 @@
+using Test
+
+using Sbplib.Grids
+using Sbplib.LazyTensors
+using Sbplib.SbpOperators
+
+stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4)
+
+struct MockOp end
+
+function SbpOperators.sat_tensors(op::MockOp, g::Grid, bc::DirichletCondition; a = 1.)
+    e = boundary_restriction(g, stencil_set, boundary(bc))
+    L = a*e
+    sat_op = e'
+    return sat_op, L
+end
+
+function SbpOperators.sat_tensors(op::MockOp, g::Grid, bc::NeumannCondition)
+    e = boundary_restriction(g, stencil_set, boundary(bc))
+    d = normal_derivative(g, stencil_set, boundary(bc))
+    L = d
+    sat_op = e'
+    return sat_op, L
+end
+
+@testset "sat" begin
+    op = MockOp()
+    @testset "1D" begin
+        grid  = equidistant_grid(0., 1., 11)
+        l, r = boundary_identifiers(grid)
+        u = eval_on(grid, x-> 1. + 2x^2)
+        dc = DirichletCondition(1.0, l)
+        g_l = discretize_data(grid, dc)
+        SAT_l = sat(op, grid, dc)
+        @test SAT_l(u, g_l) ≈ zeros((size(grid))) atol = 1e-13
+        
+        nc = NeumannCondition(4.0, r)
+        g_r = discretize_data(grid, nc)
+        SAT_r = sat(op, grid, nc)
+        @test SAT_r(u, g_r) ≈ zeros((size(grid))) atol = 1e-13
+    end
+    @testset "2D" begin
+        grid  = equidistant_grid((0.,0.), (1.,1.), 11, 13)
+        W, E, S, N = boundary_identifiers(grid)
+        u = eval_on(grid, (x,y) -> x+y^2)
+
+        dc_W = DirichletCondition(1.0, W)
+        SAT_W = sat(op, grid, dc_W)
+        g_W = discretize_data(grid, dc_W)
+        r_W = zeros(size(grid))
+        r_W[1,:] .= map(y -> (y^2-1.), range(0., 1., length=13))
+        @test SAT_W(u, g_W) ≈ r_W atol = 1e-13
+
+        dc_E = DirichletCondition(2, E)
+        SAT_E = sat(op, grid, dc_E; a = 2.)
+        g_E = discretize_data(grid, dc_E)
+        r_E = zeros(size(grid))
+        r_E[end,:] .= map(y -> (2*(1. + y^2)-2.), range(0., 1., length=13))
+        @test SAT_E(u, g_E) ≈ r_E atol = 1e-13
+
+        nc_S = NeumannCondition(.0, S)
+        SAT_S = sat(op, grid, nc_S)
+        g_S = discretize_data(grid, nc_S)
+        @test SAT_S(u, g_S) ≈ zeros(size(grid)) atol = 1e-13
+
+        nc_N = NeumannCondition(2.0, N)
+        SAT_N = sat(op, grid, nc_N)
+        g_N = discretize_data(grid, nc_N)
+        @test SAT_N(u, g_N) ≈ zeros(size(grid)) atol = 1e-13
+    end
+end

--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Sat May 25 16:02:55 2024 -0700
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Tue Jun 25 15:32:19 2024 +0200
@@ -88,3 +88,66 @@
     end
 end
 
+@testset "sat_tensors" begin
+    # TODO: The following tests should be implemented
+    #       1. Symmetry D'H == H'D (test_broken below)
+    #       2. Test eigenvalues of and/or solution to Poisson
+    #       3. Test tuning of Dirichlet conditions
+    #
+    #       These tests are likely easiest to implement once
+    #       we have support for generating matrices from tensors.
+
+    operator_path = sbp_operators_path()*"standard_diagonal.toml"
+    orders = (2,4)
+    tols = (5e-2,5e-4)
+    sz = (201,401)
+    g = equidistant_grid((0.,0.), (1.,1.), sz...)
+    
+    # Verify implementation of sat_tensors by testing accuracy and symmetry (TODO) 
+    # of the operator D = Δ + SAT, where SAT is the tensor composition of the 
+    # operators from sat_tensor. Note that SAT*u should approximate 0 for the 
+    # conditions chosen.
+
+    @testset "Dirichlet" begin
+        for (o, tol) ∈ zip(orders,tols)
+            stencil_set = read_stencil_set(operator_path; order=o)
+            Δ = Laplace(g, stencil_set)
+            H = inner_product(g, stencil_set)
+            u = collect(eval_on(g, (x,y) -> sin(π*x)sin(2*π*y)))
+            Δu = collect(eval_on(g, (x,y) -> -5*π^2*sin(π*x)sin(2*π*y)))
+            D = Δ 
+            for id ∈ boundary_identifiers(g)
+                D = D + foldl(∘, sat_tensors(Δ, g, DirichletCondition(0., id)))
+            end
+            e = D*u .- Δu
+            # Accuracy
+            @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
+            # Symmetry
+            r = randn(size(u))
+            @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
+        end
+    end
+
+    @testset "Neumann" begin
+        @testset "Dirichlet" begin
+            for (o, tol) ∈ zip(orders,tols)
+                stencil_set = read_stencil_set(operator_path; order=o)
+                Δ = Laplace(g, stencil_set)
+                H = inner_product(g, stencil_set)
+                u = collect(eval_on(g, (x,y) -> cos(π*x)cos(2*π*y)))
+                Δu = collect(eval_on(g, (x,y) -> -5*π^2*cos(π*x)cos(2*π*y)))
+                D = Δ 
+                for id ∈ boundary_identifiers(g)
+                    D = D + foldl(∘, sat_tensors(Δ, g, NeumannCondition(0., id)))
+                end
+                e = D*u .- Δu
+                # Accuracy
+                @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
+                # Symmetry
+                r = randn(size(u))
+                @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
+            end
+        end
+    end
+end
+