--- a/Notes.md	Wed Oct 12 10:44:44 2022 +0200
+++ b/Notes.md	Wed Oct 12 10:45:47 2022 +0200
@@ -1,5 +1,23 @@
 # Notes
 
+## Boundary Conditions and SATs
+
+Types for boundary conditions:
+
+ * abstract type `BoundaryDataType`
+ * abstract type `BoundaryCondition{T<:BoundaryDataType}`
+ * concrete types `ConstantBoundaryData <: BoundaryDataType` and similar
+ * concrete types `NeumannCondition{BDT<:BoundaryDataType} <: BoundaryCondition{BDT}` and similar
+The concrete `BoundaryDataType` 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{BDT}` 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)`. This is how one typically would write the SAT in the litterature. 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.
+
+It is not clear if any of these are preferable as it currently stands.
+
 ## 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 Oct 12 10:45:47 2022 +0200
@@ -0,0 +1,24 @@
+module BoundaryConditions
+
+export BoundaryDataType
+export ConstantBoundaryData
+export SpaceDependentBoundaryData
+export TimeDependentBoundaryData
+export SpaceDependentBoundaryData
+
+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 Oct 12 10:45:47 2022 +0200
@@ -0,0 +1,52 @@
+"""
+    BoundaryDataType
+
+A type for storing boundary data, e.g. constant, space-dependent, time-dependent etc.
+Subtypes of `BoundaryDataType` should store the boundary data in a field `val`, i.e.
+`struct MyBoundaryDataType{T} <: BoundaryDataType val::T end`.
+"""
+abstract type BoundaryDataType end
+
+struct ConstantBoundaryData{T} <: BoundaryDataType
+    val::T
+end
+
+struct SpaceDependentBoundaryData{T} <: BoundaryDataType
+    val::T
+end
+
+struct TimeDependentBoundaryData{T} <: BoundaryDataType
+    val::T
+end
+
+struct SpaceTimeDependentBoundaryData{T} <: BoundaryDataType
+    val::T
+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.
+"""
+# TODO: Parametrize the boundary id as well?
+abstract type BoundaryCondition{T<:BoundaryDataType} end
+
+data(bc::BoundaryCondition) = bc.data.val
+
+struct NeumannCondition{BDT<:BoundaryDataType} <: BoundaryCondition{BDT}
+    id::BoundaryIdentifier
+    data::BDT
+end
+
+struct DirichletCondition{BDT<:BoundaryDataType} <: BoundaryCondition{BDT}
+    id::BoundaryIdentifier
+    data::BDT
+end
+
+struct RobinCondition{BDT<:BoundaryDataType,T<:Real} <: BoundaryCondition{BDT}
+    id::BoundaryIdentifier
+    data::BDT
+    α::T
+    β::T
+end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/BoundaryConditions/sat.jl	Wed Oct 12 10:45:47 2022 +0200
@@ -0,0 +1,62 @@
+"""
+    sat_tensors(op, grid, bc::BoundaryCondition{T}, params...)
+
+Returns the `LazyTensor`s used to construct a SAT for the SBP operator `op` on 
+`grid` associated with the boundary condition `bc`.
+"""
+function sat_tensors end
+
+"""
+    sat(op, grid, bc::BoundaryCondition{ConstantBoundaryData{T}}, params...)
+
+Simultaneous-Approximation-Term for constant boundary data. f = sat(op, grid, bc) returns
+an anonymous function, such that f(u) is a `LazyTensorApplication` weakly imposing the BC.
+"""
+function sat(op, grid, bc::BoundaryCondition{ConstantBoundaryData{T}}, params...) where T
+    closure, penalty = sat_tensors(op, grid, bc, params...)
+    b_sz = size(boundary_grid(grid, bc.id))
+    g = fill(data(bc), b_sz)
+    if iszero(g)
+        return u -> closure(u)
+    else
+        return u -> closure(u) + penalty(g)
+    end
+end
+
+"""
+    sat(op, grid, bc::BoundaryCondition{SpaceDependentBoundaryData{T}}, params...)
+
+Simultaneous-Approximation-Term for space-dependent boundary data. f = sat(op, grid, bc) returns
+an anonymous function, such that f(u) is a `LazyTensorApplication` weakly imposing the BC.
+"""
+function sat(op, grid, bc::BoundaryCondition{SpaceDependentBoundaryData{T}}, params...) where T
+    closure, penalty = sat_tensors(op, grid, bc, params...)
+    g = data(bc)
+    return u -> closure(u) + penalty(g)
+end
+
+"""
+    sat(op, grid, bc::BoundaryCondition{SpaceDependentBoundaryData{T}}, params...)
+
+Simultaneous-Approximation-Term for time-dependent boundary data. f = sat(op, grid, bc) returns
+an anonymous function, such that f(u,t) is a `LazyTensorApplication` weakly imposing the BC at time t.
+"""
+function sat(op, grid, bc::BoundaryCondition{TimeDependentBoundaryData{T}}, params...) where T
+    closure, penalty = sat_tensors(op, grid, bc, params...)
+    b_sz = size(boundary_grid(grid, bc.id))
+    b_vec = ones(eltype(grid), b_sz)
+    g = data(bc)
+    return (u,t) -> closure(u) + g(t)*penalty(b_vec)
+end
+
+"""
+    sat(op, grid, bc::BoundaryCondition{SpaceDependentBoundaryData{T}}, params...)
+
+Simultaneous-Approximation-Term for space-time-dependent boundary data. f = sat(op, grid, bc) returns
+an anonymous function, such that f(u,t) is a `LazyTensorApplication` weakly imposing the BC at time t.
+"""
+function sat(op, grid, bc::BoundaryCondition{SpaceTimeDependentBoundaryData{T}}, params...) where T
+    closure, penalty = sat_tensors(op, grid, bc, params...)
+    g = data(bc)
+    return (u,t) -> closure(u) + penalty(g(t))
+end
\ No newline at end of file

--- a/src/LazyTensors/lazy_tensor_operations.jl	Wed Oct 12 10:44:44 2022 +0200
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Wed Oct 12 10:45:47 2022 +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	Wed Oct 12 10:44:44 2022 +0200
+++ b/src/SbpOperators/SbpOperators.jl	Wed Oct 12 10:45:47 2022 +0200
@@ -23,6 +23,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	Wed Oct 12 10:44:44 2022 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Wed Oct 12 10:45:47 2022 +0200
@@ -53,3 +53,24 @@
     end
     return Δ
 end
+
+"""
+    sat_tensors(Δ::Laplace, g::EquidistantGrid, 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::EquidistantGrid, 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	Wed Oct 12 10:44:44 2022 +0200
+++ b/src/Sbplib.jl	Wed Oct 12 10:45:47 2022 +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