--- a/src/BoundaryConditions/boundary_condition.jl	Wed Jul 26 21:35:50 2023 +0200
+++ b/src/BoundaryConditions/boundary_condition.jl	Wed Jul 26 23:11:02 2023 +0200
@@ -1,65 +1,25 @@
-# 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
+    BoundaryCondition
 
-`val` is a scalar value of type T
+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.
 """
-struct ConstantBoundaryData{T<:Number} <: BoundaryData
-    val::T
-end
+abstract type BoundaryCondition end
 
 """
-    SpaceDependentBoundaryData
-
-`val` is a function of dimensionality equal to the boundary
-"""
-struct SpaceDependentBoundaryData{T<:Function} <: BoundaryData
-    val::T
-end
+    id(::BoundaryCondition)
 
+The boundary identifier of the BoundaryCondition.
+Must be implemented by subtypes.
 """
-    TimeDependentBoundaryData
-
-`val` is a scalar function val(t)
-"""
-struct TimeDependentBoundaryData{T<:Function} <: BoundaryData
-    val::T
-end
+function id 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
+    data(::BoundaryCondition)
 
-    function SpaceTimeDependentBoundaryData(f::Function)
-        val(t) = (args...) -> f(t,args...)
-        return new{typeof(val)}(val)
-    end
-end
-
+If implemented, the data associated with the BoundaryCondition
 """
-    ZeroBoundaryData
-"""
-struct ZeroBoundaryData <: BoundaryData end
-
+function data end
 
 """
     discretize(::BoundaryData, boundary_grid)
@@ -67,55 +27,20 @@
 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))
+function discretize_data(grid, bc::BoundaryCondition)
+    return eval_on(boundary_grid(grid, id(bc)), data(bc))
 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
+struct NeumannCondition{DT} <: BoundaryCondition{DT}
+    data::DT
     id::BoundaryIdentifier 
 end
+id(bc::NeumannCondition) = bc.id
+data(bc::NeumannCondition) = bc.data
 
-struct DirichletCondition{BD<:BoundaryData} <: BoundaryCondition{BD}
-    data::BD
+struct DirichletCondition{DT} <: BoundaryCondition{DT}
+    data::DT
     id::BoundaryIdentifier
-end
\ No newline at end of file
+end
+id(bc::NeumannCondition) = bc.id
+data(bc::DirichletCondition) = bc.data
\ No newline at end of file

--- a/src/BoundaryConditions/sat.jl	Wed Jul 26 21:35:50 2023 +0200
+++ b/src/BoundaryConditions/sat.jl	Wed Jul 26 23:11:02 2023 +0200
@@ -19,23 +19,6 @@
 `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)
+    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/SbpOperators/volumeops/laplace/laplace.jl	Wed Jul 26 21:35:50 2023 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Wed Jul 26 23:11:02 2023 +0200
@@ -57,8 +57,8 @@
 """
 sat_tensors(Δ::Laplace, g::TensorGrid, bc::NeumannCondition)
 
-Returns anonymous functions for construction the `LazyTensorApplication`s
-recuired in order to impose a Neumann boundary condition.
+Returns the LazyTensors required to impose a Neumann condition
+SAT = sat_op(d*u - g)
 
 See also: [`sat`,`NeumannCondition`](@ref).
 """
@@ -70,7 +70,6 @@
     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
+    sat_tensor = H⁻¹∘e'∘Hᵧ
+    return sat_tensor, d
 end