--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/BoundaryValue.jl	Wed Sep 09 21:06:27 2020 +0200
@@ -0,0 +1,127 @@
+"""
+    BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
+
+Implements the boundary operator `e` as a TensorMapping
+"""
+struct BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
+    eClosure::Stencil{T,M}
+    bId::CartesianBoundary
+end
+export BoundaryValue
+
+# TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
+# Can we give special treatment to TensorMappings that go to a higher dim?
+function LazyTensors.range_size(e::BoundaryValue{T}, domain_size::NTuple{1,Integer}) where T
+    if dim(e.bId) == 1
+        return (UnknownDim, domain_size[1])
+    elseif dim(e.bId) == 2
+        return (domain_size[1], UnknownDim)
+    end
+end
+LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
+# TODO: Make a nicer solution for 3-dim(e.bId)
+
+# TODO: Make this independent of dimension
+function LazyTensors.apply(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T
+    i = I[dim(e.bId)]
+    j = I[3-dim(e.bId)]
+    N_i = size(e.grid)[dim(e.bId)]
+    return apply_boundary_value(e.op, v[j], i, N_i, region(e.bId))
+end
+
+function LazyTensors.apply_transpose(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T
+    u = selectdim(v,3-dim(e.bId),Int(I[1]))
+    return apply_boundary_value_transpose(e.op, u, region(e.bId))
+end
+
+function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Lower})
+    @boundscheck if length(v) < closuresize(op)
+        throw(BoundsError())
+    end
+    apply_stencil(op.eClosure,v,1)
+end
+
+function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Upper})
+    @boundscheck if length(v) < closuresize(op)
+        throw(BoundsError())
+    end
+    apply_stencil_backwards(op.eClosure,v,length(v))
+end
+export apply_boundary_value_transpose
+
+function apply_boundary_value(op::ConstantStencilOperator, v::Number, i::Index, N::Integer, ::Type{Lower})
+    @boundscheck if !(0<length(Int(i)) <= N)
+        throw(BoundsError())
+    end
+    op.eClosure[Int(i)-1]*v
+end
+
+function apply_boundary_value(op::ConstantStencilOperator, v::Number, i::Index,  N::Integer, ::Type{Upper})
+    @boundscheck if !(0<length(Int(i)) <= N)
+        throw(BoundsError())
+    end
+    op.eClosure[N-Int(i)]*v
+end
+export apply_boundary_value
+
+
+"""
+    BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
+
+Implements the boundary operator `e` as a TensorMapping
+"""
+struct BoundaryValue{D,T,M,R} <: TensorMapping{T,D,1}
+    e:BoundaryOperator{T,M,R}
+    bId::CartesianBoundary
+end
+
+function LazyTensors.apply_transpose(bv::BoundaryValue{T,M,Lower}, v::AbstractVector{T}, i::Index) where T
+    u = selectdim(v,3-dim(bv.bId),Int(I[1]))
+    return apply_transpose(bv.e, u, I)
+end
+
+
+"""
+    BoundaryOperator{T,N,R} <: TensorMapping{T,1,1}
+
+Implements the boundary operator `e` as a TensorMapping
+"""
+export BoundaryOperator
+struct BoundaryOperator{T,M,R<:Region} <: TensorMapping{T,1,1}
+    closure::Stencil{T,M}
+end
+
+function LazyTensors.range_size(e::BoundaryOperator, domain_size::NTuple{1,Integer})
+    return UnknownDim
+end
+
+LazyTensors.domain_size(e::BoundaryOperator{T}, range_size::NTuple{1,Integer}) where T = range_size
+
+function LazyTensors.apply_transpose(e::BoundaryOperator{T,M,Lower}, v::AbstractVector{T}, i::Index{Lower}) where T
+    @boundscheck if length(v) < closuresize(e) #TODO: Use domain_size here?
+        throw(BoundsError())
+    end
+    apply_stencil(e.closure,v,Int(i))
+end
+
+function LazyTensors.apply_transpose(e::BoundaryOperator{T,M,Upper}}, v::AbstractVector{T}, i::Index{Upper}) where T
+    @boundscheck if length(v) < closuresize(e) #TODO: Use domain_size here?
+        throw(BoundsError())
+    end
+    apply_stencil_backwards(e.closure,v,Int(i))
+end
+
+function LazyTensors.apply_transpose(e::BoundaryOperator{T}, v::AbstractVector{T}, i::Index) where T
+    @boundscheck if length(v) < closuresize(e) #TODO: Use domain_size here?
+        throw(BoundsError())
+    end
+    return eltype(v)(0)
+end
+
+#TODO: Implement apply in a meaningful way. Should it return a vector or a single value (perferable?) Should fit into  the
+function LazyTensors.apply(e::BoundaryOperator, v::AbstractVector, i::Index)
+    @boundscheck if !(0<length(Int(i)) <= length(v))
+        throw(BoundsError())
+    end
+    return e.closure[Int(i)].*v
+end