--- a/Notes.md	Tue Jan 23 20:39:06 2024 +0100
+++ b/Notes.md	Tue Jan 23 20:48:25 2024 +0100
@@ -1,5 +1,23 @@
 # Notes
 
+## Boundary Conditions and SATs
+
+Types for boundary conditions:
+
+ * abstract type `BoundaryData`
+ * abstract type `BoundaryCondition{T<:BoundaryData}`
+ * concrete types `ConstantBoundaryData <: BoundaryData` and similar
+ * concrete types `NeumannCondition{BD<:BoundaryData} <: BoundaryCondition{BD}` and similar
+The concrete `BoundaryData` subtypes are "thin types" wrapping the boundary data, and are used to indicate how the boundary data should be used in e.g. sat routines. The concrete `BoundaryCondition{BD}` subtypes are used for assembling the tensors used to construct e.g. a SAT.
+
+SAT methods:
+There are multiple options for what the SAT methods could return.
+* (Current) a function which returns a `LazyTensorApplication`, e.g. `f = sat(grid,op,bc)`. The the resulting `LazyTensorApplication` can then be added to scheme i.e. `scheme = op*u + f(u)`.  Depdending on the type of data in the BC, e.g. time-depdendent etc one can return f(u,t).
+* `LazyTensor`s `closure, penalty = sat(grid,op,bc)` like in the matlab version. Probably the most general one. Up to the user to make use of the returned `LazyTensor`s. One can for example then easily include the closures to the operator and have eg. `D = (op + closure)*u`.
+* A `LazyTensor` for closure, and a `LazyArray` for `penalty*data`. Mix of the above.
+* Same as first but of  the  form sat = `sat_op*(L*u-g)`. This is how one typically would write the SAT in the litterature. The function `sat_tensors` would return `sat_op` and `L`. Need to get compositions working before we can implement this approach.
+
+
 ## Reading operators
 
 Jonatan's suggestion is to add methods to `Laplace`, `SecondDerivative` and

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/BoundaryConditions/BoundaryConditions.jl	Tue Jan 23 20:48:25 2024 +0100
@@ -0,0 +1,21 @@
+module BoundaryConditions
+
+
+export BoundaryCondition
+export discretize_data
+export data
+export id
+
+export NeumannCondition
+export DirichletCondition
+
+export sat
+export sat_tensors
+
+using Sbplib.Grids
+using Sbplib.LazyTensors
+
+include("boundary_condition.jl")
+include("sat.jl")
+
+end # module

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/BoundaryConditions/boundary_condition.jl	Tue Jan 23 20:48:25 2024 +0100
@@ -0,0 +1,47 @@
+"""
+    BoundaryCondition
+
+A type for implementing data needed in order to impose a boundary condition.
+Subtypes refer to perticular types of boundary conditions, e.g. Neumann conditions.
+"""
+abstract type BoundaryCondition{T1,T2} end
+
+"""
+    id(::BoundaryCondition)
+
+The boundary identifier of the BoundaryCondition.
+Must be implemented by subtypes.
+"""
+function id end
+
+"""
+    data(::BoundaryCondition)
+
+If implemented, the data associated with the BoundaryCondition
+"""
+function data end
+
+"""
+    discretize(grid, bc::BoundaryCondition)
+
+The grid function obtained from discretizing the `bc` data on the boundary grid
+    specified the by bc `id`.
+"""
+function discretize_data(grid, bc::BoundaryCondition)
+    return eval_on(boundary_grid(grid, id(bc)), data(bc))
+end
+
+struct DirichletCondition{T1,T2} <: BoundaryCondition{T1,T2}
+    data::T1
+    id::T2
+end
+id(bc::DirichletCondition) = bc.id
+data(bc::DirichletCondition) = bc.data
+
+struct NeumannCondition{T1,T2} <: BoundaryCondition{T1,T2}
+    data::T1
+    id::T2
+end
+id(bc::NeumannCondition) = bc.id
+data(bc::NeumannCondition) = bc.data
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/BoundaryConditions/sat.jl	Tue Jan 23 20:48:25 2024 +0100
@@ -0,0 +1,24 @@
+"""
+    sat_tensors(op, grid, bc::BoundaryCondition, params...)
+
+The tensor and boundary operator used to construct a simultaneous-approximation-term
+for imposing `bc` related to `op`.
+
+For `sat_op, L  = sat_tensors(...)` then `SAT = sat_op*(L*u - g)`  where `g` 
+is the boundary data.
+"""
+function sat_tensors end
+
+
+"""
+    sat(op, grid, bc::BoundaryCondition, params...)
+
+Simultaneous-Approximation-Term for a general `bc` to `op`. 
+Returns a function `SAT(u,g)` weakly imposing `bc` when added to `op*u`.
+
+`op` must implement the function `sat_tensors`.
+"""
+function sat(op, grid, bc::BoundaryCondition, params...)
+    sat_op, L = sat_tensors(op, grid, bc, params...)
+    return SAT(u, g) = sat_op*(L*u - g)
+end
\ No newline at end of file

--- a/src/LazyTensors/lazy_tensor_operations.jl	Tue Jan 23 20:39:06 2024 +0100
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Tue Jan 23 20:48:25 2024 +0100
@@ -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{T}) where T = (-one(T))*tm
 
 """
     InflatedTensor{T,R,D} <: LazyTensor{T,R,D}

--- a/src/SbpOperators/SbpOperators.jl	Tue Jan 23 20:39:06 2024 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Tue Jan 23 20:48:25 2024 +0100
@@ -26,6 +26,7 @@
 using Sbplib.RegionIndices
 using Sbplib.LazyTensors
 using Sbplib.Grids
+using Sbplib.BoundaryConditions
 
 @enum Parity begin
     odd = -1

--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Tue Jan 23 20:39:06 2024 +0100
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Tue Jan 23 20:48:25 2024 +0100
@@ -52,3 +52,25 @@
     return Δ
 end
 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
+
+# TODO: Add sat_tensor for Diirichlet condition
+
+"""
+sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
+
+The operators required to construct the SAT for imposing Neumann condition
+
+
+See also: [`sat`,`NeumannCondition`](@ref).
+"""
+function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
+    id = bc.id
+    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)
+
+    sat_op = H⁻¹∘e'∘Hᵧ
+    return sat_op, d
+end

--- a/src/Sbplib.jl	Tue Jan 23 20:39:06 2024 +0100
+++ b/src/Sbplib.jl	Tue Jan 23 20:48:25 2024 +0100
@@ -4,6 +4,7 @@
 include("RegionIndices/RegionIndices.jl")
 include("LazyTensors/LazyTensors.jl")
 include("Grids/Grids.jl")
+include("BoundaryConditions/BoundaryConditions.jl")
 include("SbpOperators/SbpOperators.jl")
 include("DiffOps/DiffOps.jl")
 
@@ -12,5 +13,6 @@
 export Grids
 export SbpOperators
 export DiffOps
+export BoundaryConditions
 
 end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/BoundaryConditions/boundary_condition_test.jl	Tue Jan 23 20:48:25 2024 +0100
@@ -0,0 +1,25 @@
+using Test
+
+using Sbplib.BoundaryConditions
+using Sbplib.Grids
+
+@testset "BoundaryCondition" begin
+    grid_1d = equidistant_grid(11, 0.0, 1.0)
+    grid_2d = equidistant_grid((11,15), (0.0, 0.0), (1.0,1.0))
+    grid_3d = equidistant_grid((11,15,13), (0.0, 0.0, 0.0), (1.0,1.0, 1.0))
+    (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
+    @test DirichletCondition(g,id_l) isa BoundaryCondition{Float64}
+    @test DirichletCondition(g,id_n) isa BoundaryCondition{Float64}
+    @test NeumannCondition(f,id_b) isa BoundaryCondition{<:Function}
+
+    @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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/BoundaryConditions/sat_test.jl	Tue Jan 23 20:48:25 2024 +0100
@@ -0,0 +1,80 @@
+using Test
+
+
+using Sbplib.BoundaryConditions
+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 BoundaryConditions.sat_tensors(op::MockOp, g::Grid, bc::DirichletCondition)
+    e = boundary_restriction(g, stencil_set, id(bc))
+    L = e
+    sat_op = e'
+    return sat_op, L
+end
+
+function BoundaryConditions.sat_tensors(op::MockOp, g::Grid, bc::DirichletCondition, a)
+    e = boundary_restriction(g, stencil_set, id(bc))
+    L = a*e
+    sat_op = e'
+    return sat_op, L
+end
+
+function BoundaryConditions.sat_tensors(op::MockOp, g::Grid, bc::NeumannCondition)
+    e = boundary_restriction(g, stencil_set, id(bc))
+    d = normal_derivative(g, stencil_set, id(bc))
+    L = d
+    sat_op = e'
+    return sat_op, L
+end
+
+@testset "sat" begin
+    op = MockOp()
+    @testset "1D" begin
+        grid  = equidistant_grid(11, 0., 1.)
+        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((11,13), (0.,0.), (1.,1.))
+        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, 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	Tue Jan 23 20:39:06 2024 +0100
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Tue Jan 23 20:48:25 2024 +0100
@@ -3,6 +3,7 @@
 using Sbplib.SbpOperators
 using Sbplib.Grids
 using Sbplib.LazyTensors
+using Sbplib.BoundaryConditions
 
 @testset "Laplace" begin
     # Default stencils (4th order)
@@ -88,3 +89,49 @@
     end
 end
 
+@testset "sat_tensors" begin
+    operator_path = sbp_operators_path()*"standard_diagonal.toml"
+    stencil_set = read_stencil_set(operator_path; order=4)
+    g = equidistant_grid((101,102), (-1.,-1.), (1.,1.))
+    W,E,S,N = boundary_identifiers(g)
+    
+    u = eval_on(g, (x,y) -> sin(x+y))
+    uWx = eval_on(boundary_grid(g,W), (x,y) -> -cos(x+y))
+    uEx = eval_on(boundary_grid(g,E), (x,y) -> cos(x+y))
+    uSy = eval_on(boundary_grid(g,S), (x,y) -> -cos(x+y))
+    uNy = eval_on(boundary_grid(g,N), (x,y) -> cos(x+y))
+    
+    
+    v = eval_on(g, (x,y) -> cos(x+y))
+    vW = eval_on(boundary_grid(g,W), (x,y) -> cos(x+y))
+    vE = eval_on(boundary_grid(g,E), (x,y) -> cos(x+y))
+    vS = eval_on(boundary_grid(g,S), (x,y) -> cos(x+y))
+    vN = eval_on(boundary_grid(g,N), (x,y) -> cos(x+y))
+
+
+    @testset "Neumann" begin
+        Δ = Laplace(g, stencil_set)
+        H = inner_product(g, stencil_set)
+        HW = inner_product(boundary_grid(g,W), stencil_set)
+        HE = inner_product(boundary_grid(g,E), stencil_set)
+        HS = inner_product(boundary_grid(g,S), stencil_set)
+        HN = inner_product(boundary_grid(g,N), stencil_set)
+        
+        ncW = NeumannCondition(0., W)
+        ncE = NeumannCondition(0., E)
+        ncS = NeumannCondition(0., S)
+        ncN = NeumannCondition(0., N)
+        
+        SATW = foldl(∘,sat_tensors(Δ, g, ncW))
+        SATE = foldl(∘,sat_tensors(Δ, g, ncE))
+        SATS = foldl(∘,sat_tensors(Δ, g, ncS))
+        SATN = foldl(∘,sat_tensors(Δ, g, ncN))
+        
+
+        @test sum((H*SATW*u).*v) ≈ sum((HW*uWx).*vW) rtol = 1e-6
+        @test sum((H*SATE*u).*v) ≈ sum((HE*uEx).*vE) rtol = 1e-6
+        @test sum((H*SATS*u).*v) ≈ sum((HS*uSy).*vS) rtol = 1e-6
+        @test sum((H*SATN*u).*v) ≈ sum((HN*uNy).*vN) rtol = 1e-6
+    end
+end
+