--- a/Notes.md	Fri Aug 25 08:49:07 2023 +0200
+++ b/Notes.md	Fri Aug 25 15:41:56 2023 +0200
@@ -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	Fri Aug 25 15:41:56 2023 +0200
@@ -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	Fri Aug 25 15:41:56 2023 +0200
@@ -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	Fri Aug 25 15:41:56 2023 +0200
@@ -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 Aug 25 08:49:07 2023 +0200
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Fri Aug 25 15:41:56 2023 +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{T}) where T = (-one(T))*tm
 
 """
     InflatedTensor{T,R,D} <: LazyTensor{T,R,D}

--- a/src/SbpOperators/SbpOperators.jl	Fri Aug 25 08:49:07 2023 +0200
+++ b/src/SbpOperators/SbpOperators.jl	Fri Aug 25 15:41:56 2023 +0200
@@ -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 Aug 25 08:49:07 2023 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Fri Aug 25 15:41:56 2023 +0200
@@ -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 Aug 25 08:49:07 2023 +0200
+++ b/src/Sbplib.jl	Fri Aug 25 15:41:56 2023 +0200
@@ -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	Fri Aug 25 15:41:56 2023 +0200
@@ -0,0 +1,18 @@
+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) = x^2+y^2
+    @test DirichletCondition(g,id_n) isa BoundaryCondition{Float64}
+    @test NeumannCondition(f,id_n) isa BoundaryCondition{<:Function}
+
+   @test g*ones(11,1) ≈ discretize_data(grid_2d,DirichletCondition(g,id_n))
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/BoundaryConditions/sat_test.jl	Fri Aug 25 15:41:56 2023 +0200
@@ -0,0 +1,40 @@
+using Test
+
+
+using Sbplib.BoundaryConditions
+using Sbplib.Grids
+using Sbplib.RegionIndices
+using Sbplib.LazyTensors
+
+grid = equidistant_grid(11, 0.0, 1.0)
+(id_l,id_r) = boundary_identifiers(grid)
+struct MockOp
+end
+
+function BoundaryConditions.sat_tensors(op::MockOp, grid, bc::DirichletCondition)
+    sz = size(grid)
+    m = sz[1]
+    ind = (region(bc.id) == Lower()) ? 1 : m
+    e = zeros(m);
+    e[ind] = 1.
+    eᵀ = ones(Float64,m,0);
+    e[ind] = 1.
+    c_tensor = LazyTensors.DiagonalTensor(e)
+    p_tensor = DenseTensor(eᵀ, (1,), (2,))
+    closure(u) = c_tensor*u
+    function penalty(g)
+        @show g
+        return p_tensor*g
+    end
+    return closure, penalty
+end
+
+
+# @testset "sat" begin
+#     g = ConstantBoundaryData(2.0)
+#     dc = DirichletCondition(g,id_l)
+#     op = MockOp()
+#     f = sat(op, grid, dc)
+#     u = eval_on(grid, x-> -1/2 + x^2)
+#     @show f(0.,u)
+# end