--- a/Notes.md	Fri Dec 01 11:48:04 2023 +0100
+++ b/Notes.md	Sat Dec 23 23:03:13 2023 +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	Sat Dec 23 23:03:13 2023 +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	Sat Dec 23 23:03:13 2023 +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{T} 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(::BoundaryData, boundary_grid)
+
+Returns an anonymous time-dependent function f, such that f(t) is
+a `LazyArray` holding the `BoundaryData` discretized on `boundary_grid`.
+"""
+function discretize_data(grid, bc::BoundaryCondition)
+    return eval_on(boundary_grid(grid, id(bc)), data(bc))
+end
+
+struct DirichletCondition{T} <: BoundaryCondition{T}
+    data::T
+    id::BoundaryIdentifier
+end
+id(bc::DirichletCondition) = bc.id
+data(bc::DirichletCondition) = bc.data
+
+struct NeumannCondition{T} <: BoundaryCondition{T}
+    data::T
+    id::BoundaryIdentifier 
+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	Sat Dec 23 23:03:13 2023 +0100
@@ -0,0 +1,24 @@
+"""
+    sat_tensors(op, grid, bc::BoundaryCondition, params...)
+
+Returns the functions `closure(u)` and `penalty(g)` used to construct a SAT for the
+`LazyTensor` operator `op` on `grid` associated with the boundary condition `bc`,
+where g is the discretized data of `bc`.
+"""
+function sat_tensors end
+
+
+"""
+    sat(op, grid, bc::BoundaryCondition, params...)
+
+Simultaneous-Approximation-Term for a general `BoundaryCondition` `bc` to `LazyTensor` `op`. 
+The function returns a function `f`, where f(t,u)` returns a `LazyTensorApplication`
+weakly imposing the boundary condition at time `t`, when added to `op*u`.
+
+`op` must implement the function `sat_tensors`. `f` is then constructed as 
+`f(t,u) = closure(u) + `penalty(g(t))`.
+"""
+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	Fri Dec 01 11:48:04 2023 +0100
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Sat Dec 23 23:03:13 2023 +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	Fri Dec 01 11:48:04 2023 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Sat Dec 23 23:03:13 2023 +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	Fri Dec 01 11:48:04 2023 +0100
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Sat Dec 23 23:03:13 2023 +0100
@@ -52,3 +52,24 @@
     return Δ
 end
 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
+
+
+"""
+sat_tensors(Δ::Laplace, g::TensorGrid, bc::NeumannCondition)
+
+Returns the LazyTensors required to impose a Neumann condition
+SAT = sat_op(d*u - g)
+
+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_tensor = H⁻¹∘e'∘Hᵧ
+    return sat_tensor, d
+end

--- a/src/Sbplib.jl	Fri Dec 01 11:48:04 2023 +0100
+++ b/src/Sbplib.jl	Sat Dec 23 23:03:13 2023 +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	Sat Dec 23 23:03:13 2023 +0100
@@ -0,0 +1,21 @@
+using Test
+
+using Sbplib.BoundaryConditions
+using Sbplib.Grids
+
+@testset "BoundaryCondition" begin
+    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_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_n) isa BoundaryCondition{Float64}
+    @test NeumannCondition(f,id_b) isa BoundaryCondition{<:Function}
+
+    @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	Sat Dec 23 23:03:13 2023 +0100
@@ -0,0 +1,66 @@
+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::TensorGrid, 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::TensorGrid, 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::TensorGrid, 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()
+    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_S = sat(op, grid, nc_N)
+    g_N = discretize_data(grid, nc_N)
+    @test SAT_S(u, g_N) ≈ zeros(size(grid)) atol = 1e-13
+end