--- a/Notes.md	Thu Jun 08 15:52:22 2023 +0200
+++ b/Notes.md	Wed Jul 26 21:35:50 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	Wed Jul 26 21:35:50 2023 +0200
@@ -0,0 +1,27 @@
+module BoundaryConditions
+
+export BoundaryData
+export ConstantBoundaryData
+export SpaceDependentBoundaryData
+export TimeDependentBoundaryData
+export SpaceTimeDependentBoundaryData
+export ZeroBoundaryData
+export discretize
+
+export BoundaryCondition
+export NeumannCondition
+export DirichletCondition
+export RobinCondition
+
+export data
+
+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	Wed Jul 26 21:35:50 2023 +0200
@@ -0,0 +1,121 @@
+# TODO:  Should BoundaryData just be used for traits
+# of the BoundaryConditions? Seems like one then  could move the
+# the  boundary data value val directly to BoundaryCondition
+# Not sure how one would do this tho.
+"""
+    BoundaryData
+
+A type for storing boundary data, e.g. constant, space-dependent, time-dependent etc.
+Subtypes of `BoundaryData` should store the boundary data in a field `val`. The exception
+to this is ZeroBoundaryData.
+"""
+abstract type BoundaryData end
+
+"""
+    ConstantBoundaryData
+
+`val` is a scalar value of type T
+"""
+struct ConstantBoundaryData{T<:Number} <: BoundaryData
+    val::T
+end
+
+"""
+    SpaceDependentBoundaryData
+
+`val` is a function of dimensionality equal to the boundary
+"""
+struct SpaceDependentBoundaryData{T<:Function} <: BoundaryData
+    val::T
+end
+
+"""
+    TimeDependentBoundaryData
+
+`val` is a scalar function val(t)
+"""
+struct TimeDependentBoundaryData{T<:Function} <: BoundaryData
+    val::T
+end
+
+"""
+    SpaceTimeDependentBoundaryData
+
+`val` is a timedependent function returning the spacedependent
+    boundary data at a specific time. For instance, if f(t,x)
+    is the function describing the spacetimedependent  boundary data then
+    val(t*) returns the function g(x) = f(t*,x...)
+"""
+struct SpaceTimeDependentBoundaryData{T<:Function} <: BoundaryData
+    val::T
+
+    function SpaceTimeDependentBoundaryData(f::Function)
+        val(t) = (args...) -> f(t,args...)
+        return new{typeof(val)}(val)
+    end
+end
+
+"""
+    ZeroBoundaryData
+"""
+struct ZeroBoundaryData <: BoundaryData 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`.
+"""
+# TODO: Is the return type of discretize really a good interface
+# for the boundary data?
+# Moreover, instead of explicitly converting to a LazyArray here
+# should we defer this to eval_on (and extend eval_on for scalars as well)?
+# I.e. if eval_on returns a LazyArray, the boundary data is lazy. Otherwise
+# it is preallocated.
+
+function discretize(bd::ConstantBoundaryData, boundary_grid)
+    return t -> LazyTensors.LazyConstantArray(bd.val, size(boundary_grid))
+end
+
+function discretize(bd::TimeDependentBoundaryData, boundary_grid)
+    return t -> LazyTensors.LazyConstantArray(bd.val(t), size(boundary_grid))
+end
+
+function discretize(bd::SpaceDependentBoundaryData, boundary_grid)
+    return t -> eval_on(boundary_grid, bd.val)
+end
+
+function discretize(bd::SpaceTimeDependentBoundaryData, boundary_grid)
+    return t -> eval_on(boundary_grid, bd.val(t))
+end
+
+function discretize(::ZeroBoundaryData, boundary_grid)
+    return t -> LazyTensors.LazyConstantArray(zero(eltype(boundary_grid)), size(boundary_grid))
+end
+
+"""
+    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<:BoundaryData} end
+
+"""
+    data(::BoundaryCondition)
+
+Returns the data stored by the `BoundaryCondition`.
+"""
+data(bc::BoundaryCondition) = bc.data
+ 
+
+struct NeumannCondition{BD<:BoundaryData} <: BoundaryCondition{BD}
+    data::BD
+    id::BoundaryIdentifier 
+end
+
+struct DirichletCondition{BD<:BoundaryData} <: BoundaryCondition{BD}
+    data::BD
+    id::BoundaryIdentifier
+end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/BoundaryConditions/sat.jl	Wed Jul 26 21:35:50 2023 +0200
@@ -0,0 +1,41 @@
+"""
+    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...)
+    closure, penalty = sat_tensors(op, grid, bc, params...)
+    g = discretize(data(bc),boundary_grid(grid, bc.id))
+    return (t,u) -> closure(u) + penalty(g(t))
+end
+
+
+"""
+    sat(op, grid, bc::BoundaryCondition{ZeroBoundaryData}, 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 a homogenous boundary condition, when added to `op*u`.
+
+`op` must implement the function `sat_tensors`. `f` is then constructed as 
+`f(t,u) = closure(u)`.
+"""
+function sat(op, grid, bc::BoundaryCondition{ZeroBoundaryData}, params...)
+    closure = sat_tensors(op, grid, bc, params...)
+    return (t,u) -> closure(u)
+end
\ No newline at end of file

--- a/src/LazyTensors/lazy_tensor_operations.jl	Thu Jun 08 15:52:22 2023 +0200
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Wed Jul 26 21:35:50 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	Thu Jun 08 15:52:22 2023 +0200
+++ b/src/SbpOperators/SbpOperators.jl	Wed Jul 26 21:35:50 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	Thu Jun 08 15:52:22 2023 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Wed Jul 26 21:35:50 2023 +0200
@@ -52,3 +52,25 @@
     return Δ
 end
 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
+
+
+"""
+sat_tensors(Δ::Laplace, g::TensorGrid, bc::NeumannCondition)
+
+Returns anonymous functions for construction the `LazyTensorApplication`s
+recuired in order to impose a Neumann boundary 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)
+
+    closure(u) = H⁻¹*e'*Hᵧ*d*u
+    penalty(g) = -H⁻¹*e'*Hᵧ*g
+    return closure, penalty
+end

--- a/src/Sbplib.jl	Thu Jun 08 15:52:22 2023 +0200
+++ b/src/Sbplib.jl	Wed Jul 26 21:35:50 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	Wed Jul 26 21:35:50 2023 +0200
@@ -0,0 +1,99 @@
+using Test
+
+using Sbplib.BoundaryConditions
+using Sbplib.Grids
+
+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)
+
+# @testset "BoundaryData" begin
+    
+#     @testset "ConstantBoundaryData" begin
+#         c = float(pi)
+#         @test ConstantBoundaryData(c) isa BoundaryData
+#         g_1D = discretize(ConstantBoundaryData(c),boundary_grid(grid_1D, id_l))
+#         g_2D = discretize(ConstantBoundaryData(c),boundary_grid(grid_2D, id_n))
+#         @test g_1D isa Function
+#         @test g_2D isa Function
+#         @test g_1D(0.) == fill(c)
+#         @test g_2D(2.) == c*ones(11)
+#         @test_throws MethodError g_1D(0.,0.)
+#         @test_throws MethodError g_2D(0.,0.)
+#     end
+
+#     @testset "TimeDependentBoundaryData" begin
+#         f(t) = 1. /(t+0.1)
+#         @test TimeDependentBoundaryData(f) isa BoundaryData
+#         g_1D = discretize(TimeDependentBoundaryData(f),boundary_grid(grid_1D, id_l))
+#         g_2D = discretize(TimeDependentBoundaryData(f),boundary_grid(grid_2D, id_n))
+#         @test g_1D isa Function
+#         @test g_2D isa Function
+#         @test g_1D(0.) == f(0.)*fill(1)
+#         @test g_2D(2.) == f(2.)*ones(11)
+#         @test_throws MethodError g_1D(0.,0.)
+#         @test_throws MethodError g_2D(0.,0.)
+#     end
+    
+#     #TBD: Is it reasoanble to have SpaceDependentBoundaryData for 1D-grids? It would then be a constant
+#     #     which then may be represented by ConstantBoundaryData.
+#     @testset "SpaceDependentBoundaryData" begin
+#         f0(x) = 2
+#         f1(x,y) = x.^2
+#         f2(x,y,z) = x.^2 - y
+#         @test SpaceDependentBoundaryData(f1) isa BoundaryData
+#         g_1D = discretize(SpaceDependentBoundaryData(f0),boundary_grid(grid_1D, id_l))
+#         g_2D = discretize(SpaceDependentBoundaryData(f1),boundary_grid(grid_2D, id_n))
+#         g_3D = discretize(SpaceDependentBoundaryData(f2),boundary_grid(grid_3D, id_n))
+#         @test g_1D isa Function
+#         @test g_2D isa Function
+#         @test g_3D isa Function
+#         @test g_1D(1.) == fill(f0()) # Does not work since eval_on for f0 returns ().
+#         @test g_2D(2.) ≈ f1.(range(0., 1., 11)) rtol=1e-14
+#         @test g_3D(0.) ≈ eval_on(boundary_grid(grid_3D, id_n),f2) rtol=1e-14
+#         @test_throws MethodError g_1D(0.,0.)
+#         @test_throws MethodError g_2D(0.,0.)
+#         @test_throws MethodError g_3D(0.,0.)
+#     end
+    
+#     # TBD: Include tests for 1D-grids? See TBD above
+#     @testset "SpaceTimeDependentBoundaryData" begin
+#         fx1(x) = x.^2
+#         fx2(x,y) = x.^2 - y
+#         ft(t) = exp(t)
+#         ftx1(t,x) = ft(t)*fx1(x)
+#         ftx2(t,x,y) = ft(t)*fx2(x,y)
+#         @test SpaceTimeDependentBoundaryData(ftx1) isa BoundaryData
+#         g_2D = discretize(SpaceTimeDependentBoundaryData(ftx1),boundary_grid(grid_2D, id_n))
+#         g_3D = discretize(SpaceTimeDependentBoundaryData(ftx2),boundary_grid(grid_3D, id_b))
+#         @test g_2D isa Function
+#         @test g_3D isa Function
+#         @test g_2D(2.) ≈ ft(2.)*fx1.(range(0., 1., 11)) rtol=1e-14
+#         @test g_3D(3.14) ≈ ft(3.14)*eval_on(boundary_grid(grid_3D, id_b),fx2) rtol=1e-14
+#         @test_throws MethodError g_2D(0.,0.)
+#         @test_throws MethodError g_3D(0.,0.)
+#     end
+
+#     @testset "ZeroBoundaryData" begin
+#         @test ZeroBoundaryData() isa BoundaryData
+#         g_2D = discretize(ZeroBoundaryData(), boundary_grid(grid_2D, id_n))
+#         g_3D = discretize(ZeroBoundaryData(), boundary_grid(grid_3D, id_b))
+#         @test g_2D isa Function
+#         @test g_3D isa Function
+#         @test g_2D(2.) ≈ 0.0*range(0., 1., 11) rtol=1e-14
+#         f(x,y,z) = 0
+#         @test g_3D(3.14) ≈ 0.0*eval_on(boundary_grid(grid_3D, id_b), f) rtol=1e-14
+#         @test_throws MethodError g_2D(0.,0.)
+#         @test_throws MethodError g_3D(0.,0.)
+#     end
+# end
+
+# @testset "BoundaryCondition" begin
+#     g = ConstantBoundaryData(1.0)
+#     NeumannCondition(g,id_n) isa BoundaryCondition{ConstantBoundaryData}
+#     DirichletCondition(g,id_n) isa BoundaryCondition{ConstantBoundaryData}
+#     @test data(NeumannCondition(g,id_n)) == g
+# end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/BoundaryConditions/sat_test.jl	Wed Jul 26 21:35:50 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