--- a/Notes.md	Wed Dec 07 21:56:00 2022 +0100
+++ b/Notes.md	Tue Feb 07 21:55:07 2023 +0100
@@ -4,19 +4,19 @@
 
 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.
+ * 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)`. 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).
+* (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.
 
-It is not clear if any of these are preferable as it currently stands.
 
 ## Reading operators
 

--- a/src/BoundaryConditions/sat.jl	Wed Dec 07 21:56:00 2022 +0100
+++ b/src/BoundaryConditions/sat.jl	Tue Feb 07 21:55:07 2023 +0100
@@ -1,65 +1,41 @@
 """
-    sat_tensors(op, grid, bc::BoundaryCondition{T}, params...)
+    sat_tensors(op, grid, bc::BoundaryCondition, params...)
 
-Returns the `LazyTensor`s used to construct a SAT for the SBP operator `op` on 
-`grid` associated with the boundary condition `bc`.
+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
 
-# TODO: Docs must be more specific in what this function does...
+
 """
     sat(op, grid, bc::BoundaryCondition, params...)
 
-Simultaneous-Approximation-Term for general BoundaryCondition bc. f = sat(op, grid, bc) returns
-an anonymous function, such that f(t,u) is a `LazyTensorApplication` weakly imposing bc
-at time t.
+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...)
-    data_array = discretize(data(bc),boundary_grid(grid, bc.id))
-    return (t,u) -> closure(u) + penalty(data_array(t))
-end
-
-function sat(op, grid, bc::BoundaryCondition{ZeroBoundaryData}, params...)
-    closure = sat_tensors(op, grid, bc, params...)
-    return (t,u) -> closure(u)
+    g = discretize(data(bc),boundary_grid(grid, bc.id))
+    return (t,u) -> closure(u) + penalty(g(t))
 end
 
 
-# """
-#     sat(op, grid, bc::BoundaryCondition{SpaceDependentBoundaryData{T}}, params...)
+"""
+    sat(op, grid, bc::BoundaryCondition{ZeroBoundaryData}, 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 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`.
 
-# 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
+`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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/BoundaryConditions/sat_test.jl	Tue Feb 07 21:55:07 2023 +0100
@@ -0,0 +1,40 @@
+using Test
+
+
+using Sbplib.BoundaryConditions
+using Sbplib.Grids
+using Sbplib.RegionIndices
+using Sbplib.LazyTensors
+
+grid = EquidistantGrid(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 = evalOn(grid, x-> -1/2 + x^2)
+    @show f(0.,u)
+end