--- a/src/BoundaryConditions/BoundaryConditions.jl	Wed Oct 12 10:45:47 2022 +0200
+++ b/src/BoundaryConditions/BoundaryConditions.jl	Wed Dec 07 21:39:07 2022 +0100
@@ -1,10 +1,12 @@
 module BoundaryConditions
 
-export BoundaryDataType
+export BoundaryData
 export ConstantBoundaryData
 export SpaceDependentBoundaryData
 export TimeDependentBoundaryData
-export SpaceDependentBoundaryData
+export SpaceTimeDependentBoundaryData
+export ZeroBoundaryData
+export discretize
 
 export BoundaryCondition
 export NeumannCondition
@@ -12,6 +14,7 @@
 export RobinCondition
 
 export data
+
 export sat
 export sat_tensors
 

--- a/src/BoundaryConditions/boundary_condition.jl	Wed Oct 12 10:45:47 2022 +0200
+++ b/src/BoundaryConditions/boundary_condition.jl	Wed Dec 07 21:39:07 2022 +0100
@@ -1,26 +1,97 @@
+# 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.
 """
-    BoundaryDataType
+    BoundaryData
 
 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`.
+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
+
 """
-abstract type BoundaryDataType end
+    SpaceDependentBoundaryData
 
-struct ConstantBoundaryData{T} <: BoundaryDataType
+`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
 
-struct SpaceDependentBoundaryData{T} <: BoundaryDataType
+"""
+    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
 
-struct TimeDependentBoundaryData{T} <: BoundaryDataType
-    val::T
+"""
+    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 evalOn (and extend evalOn for scalars as well)?
+# I.e. if evalOn 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
 
-struct SpaceTimeDependentBoundaryData{T} <: BoundaryDataType
-    val::T
+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 -> evalOn(boundary_grid, bd.val)
+end
+
+function discretize(bd::SpaceTimeDependentBoundaryData, boundary_grid)
+    return t -> evalOn(boundary_grid, bd.val(t))
+end
+
+function discretize(::ZeroBoundaryData, boundary_grid)
+    return t -> LazyTensors.LazyConstantArray(zero(eltype(boundary_grid)), size(boundary_grid))
 end
 
 """
@@ -29,24 +100,22 @@
 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
+abstract type BoundaryCondition{T<:BoundaryData} end
+
+"""
+    data(::BoundaryCondition)
 
-data(bc::BoundaryCondition) = bc.data.val
+Returns the data stored by the `BoundaryCondition`.
+"""
+data(bc::BoundaryCondition) = bc.data
+ 
 
-struct NeumannCondition{BDT<:BoundaryDataType} <: BoundaryCondition{BDT}
-    id::BoundaryIdentifier
-    data::BDT
+struct NeumannCondition{BD<:BoundaryData} <: BoundaryCondition{BD}
+    data::BD
+    id::BoundaryIdentifier 
 end
 
-struct DirichletCondition{BDT<:BoundaryDataType} <: BoundaryCondition{BDT}
+struct DirichletCondition{BD<:BoundaryData} <: BoundaryCondition{BD}
+    data::BD
     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/test/BoundaryConditions/boundary_condition_test.jl	Wed Dec 07 21:39:07 2022 +0100
@@ -0,0 +1,99 @@
+using Test
+
+using Sbplib.BoundaryConditions
+using Sbplib.Grids
+
+grid_1D = EquidistantGrid(11, 0.0, 1.0)
+grid_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
+grid_3D = EquidistantGrid((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() = 2
+        f1(x) = x.^2
+        f2(x,y) = 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_broken g_1D(1.) == fill(f0()) # Does not work since evalOn for f0 returns ().
+        @test g_2D(2.) ≈ f1.(range(0., 1., 11)) rtol=1e-14
+        @test g_3D(0.) ≈ evalOn(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)*evalOn(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) = 0
+        @test g_3D(3.14) ≈ 0.0*evalOn(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