--- a/src/SbpOperators/BoundaryValue.jl	Mon Nov 30 09:13:13 2020 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,127 +0,0 @@
-"""
-    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

--- a/src/SbpOperators/SbpOperators.jl	Mon Nov 30 09:13:13 2020 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -14,5 +14,6 @@
 include("quadrature/quadrature.jl")
 include("quadrature/inverse_diagonal_inner_product.jl")
 include("quadrature/inverse_quadrature.jl")
+include("boundaryops/boundary_restriction.jl")
 
 end # module

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -0,0 +1,74 @@
+"""
+    boundary_restriction(grid,closureStencil,boundary)
+
+Creates a BoundaryRestriction operator for the specified boundary
+"""
+function boundary_restriction(grid::EquidistantGrid{1,T}, closureStencil::Stencil{T,M}, boundary::CartesianBoundary{1}) where {T,M}
+    r = region(boundary)
+    return e = BoundaryRestriction(grid, closureStencil, r())
+end
+
+function boundary_restriction(grid::EquidistantGrid{2,T}, closureStencil::Stencil{T,M}, boundary::CartesianBoundary{1}) where {T,M}
+    r = region(boundary)
+    e = BoundaryRestriction(restrict(grid, 1), closureStencil, r())
+    I = IdentityMapping{T}(size(restrict(grid,2)))
+    return e⊗I
+end
+
+function boundary_restriction(grid::EquidistantGrid{2,T}, closureStencil::Stencil{T,M}, boundary::CartesianBoundary{2}) where {T,M}
+    r = region(boundary)
+    e = BoundaryRestriction(restrict(grid, 2), closureStencil, r())
+    I = IdentityMapping{T}(size(restrict(grid,1)))
+    return I⊗e
+end
+export boundary_restriction
+
+"""
+    BoundaryRestriction{T,N,R} <: TensorMapping{T,0,1}
+
+Implements the boundary operator `e` as a TensorMapping
+"""
+struct BoundaryRestriction{T,M,R<:Region} <: TensorMapping{T,0,1}
+    stencil::Stencil{T,M}
+    size::NTuple{1,Int}
+end
+export BoundaryRestriction
+
+function BoundaryRestriction(grid::EquidistantGrid{1,T}, closureStencil::Stencil{T,M}, region::Region) where {T,M,R}
+    return BoundaryRestriction{T,M,typeof(region)}(closureStencil,size(grid))
+end
+
+LazyTensors.range_size(e::BoundaryRestriction) = ()
+LazyTensors.domain_size(e::BoundaryRestriction) = e.size
+
+# TODO: Should we support indexing into the 0-dimensional lazyarray? This is
+# supported for arrays with linear index style (i.e for e.g
+# u = fill(1), u[] and u[1] are both valid.) This currently not supported by
+# LazyTensorMappingApplication.
+" Restricts a grid function v on a grid of size m to the scalar element v[1]"
+function LazyTensors.apply(e::BoundaryRestriction{T,M,Lower}, v::AbstractVector{T}) where {T,M}
+    apply_stencil(e.stencil,v,1)
+end
+
+" Restricts a grid function v on a grid of size m to the scalar element v[m]"
+function LazyTensors.apply(e::BoundaryRestriction{T,M,Upper}, v::AbstractVector{T}) where {T,M}
+    apply_stencil_backwards(e.stencil,v,e.size[1])
+end
+
+" Transpose of a restriction is an inflation or prolongation.
+  Inflates the scalar (1-element) vector to a vector of size of the grid"
+function LazyTensors.apply_transpose(e::BoundaryRestriction{T,M,Lower}, v::AbstractArray{T,0}, i) where {T,M}
+    @boundscheck if !(0 < Int(i) <= e.size[1])
+        throw(BoundsError())
+    end
+    return e.stencil[Int(i)-1]*v[]
+end
+
+" Transpose of a restriction is an inflation or prolongation.
+  Inflates the scalar (1-element) vector to a vector of size of the grid"
+function LazyTensors.apply_transpose(e::BoundaryRestriction{T,M,Upper}, v::AbstractArray{T,0}, i) where {T,M}
+    @boundscheck if !(0 < Int(i) <= e.size[1])
+        throw(BoundsError())
+    end
+    return e.stencil[e.size[1] - Int(i)]*v[]
+end

--- a/test/testSbpOperators.jl	Mon Nov 30 09:13:13 2020 +0100
+++ b/test/testSbpOperators.jl	Mon Nov 30 18:30:24 2020 +0100
@@ -172,67 +172,86 @@
     @test_broken Qinv*(Q*v) ≈ v
     @test Qinv*v == Qinv'*v
 end
-#
-# @testset "BoundaryValue" begin
-#     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
-#     g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
-#
-#     e_w = BoundaryValue(op, g, CartesianBoundary{1,Lower}())
-#     e_e = BoundaryValue(op, g, CartesianBoundary{1,Upper}())
-#     e_s = BoundaryValue(op, g, CartesianBoundary{2,Lower}())
-#     e_n = BoundaryValue(op, g, CartesianBoundary{2,Upper}())
-#
-#     v = zeros(Float64, 4, 5)
-#     v[:,5] = [1, 2, 3,4]
-#     v[:,4] = [1, 2, 3,4]
-#     v[:,3] = [4, 5, 6, 7]
-#     v[:,2] = [7, 8, 9, 10]
-#     v[:,1] = [10, 11, 12, 13]
-#
-#     @test e_w  isa TensorMapping{T,2,1} where T
-#     @test e_w' isa TensorMapping{T,1,2} where T
-#
-#     @test domain_size(e_w, (3,2)) == (2,)
-#     @test domain_size(e_e, (3,2)) == (2,)
-#     @test domain_size(e_s, (3,2)) == (3,)
-#     @test domain_size(e_n, (3,2)) == (3,)
-#
-#     @test size(e_w'*v) == (5,)
-#     @test size(e_e'*v) == (5,)
-#     @test size(e_s'*v) == (4,)
-#     @test size(e_n'*v) == (4,)
-#
-#     @test e_w'*v == [10,7,4,1.0,1]
-#     @test e_e'*v == [13,10,7,4,4.0]
-#     @test e_s'*v == [10,11,12,13.0]
-#     @test e_n'*v == [1,2,3,4.0]
-#
-#     g_x = [1,2,3,4.0]
-#     g_y = [5,4,3,2,1.0]
-#
-#     G_w = zeros(Float64, (4,5))
-#     G_w[1,:] = g_y
-#
-#     G_e = zeros(Float64, (4,5))
-#     G_e[4,:] = g_y
-#
-#     G_s = zeros(Float64, (4,5))
-#     G_s[:,1] = g_x
-#
-#     G_n = zeros(Float64, (4,5))
-#     G_n[:,5] = g_x
-#
-#     @test size(e_w*g_y) == (UnknownDim,5)
-#     @test size(e_e*g_y) == (UnknownDim,5)
-#     @test size(e_s*g_x) == (4,UnknownDim)
-#     @test size(e_n*g_x) == (4,UnknownDim)
-#
-#     # These tests should be moved to where they are possible (i.e we know what the grid should be)
-#     @test_broken e_w*g_y == G_w
-#     @test_broken e_e*g_y == G_e
-#     @test_broken e_s*g_x == G_s
-#     @test_broken e_n*g_x == G_n
-# end
+
+@testset "BoundaryRestrictrion" begin
+    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+
+    g = EquidistantGrid(4, 0.0, 1.0)
+
+    e_l = BoundaryRestriction(g,op.eClosure,Lower())
+    e_r = BoundaryRestriction(g,op.eClosure,Upper())
+
+    v = evalOn(g,x->1+x^2)
+    u = fill(3.124)
+
+    @test (e_l*v)[] == v[1]
+    @test (e_r*v)[] == v[end]
+    @test e_l'*u == [u[], 0, 0, 0]
+    @test e_r'*u == [0, 0, 0, u[]]
+    @test_throws BoundsError (e_l*v)[Index{Lower}(3)]
+    @test_throws BoundsError (e_r*v)[Index{Upper}(3)]
+
+    g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
+
+    e_w = boundary_restriction(g, op.eClosure, CartesianBoundary{1,Lower}())
+    e_e = boundary_restriction(g, op.eClosure, CartesianBoundary{1,Upper}())
+    e_s = boundary_restriction(g, op.eClosure, CartesianBoundary{2,Lower}())
+    e_n = boundary_restriction(g, op.eClosure, CartesianBoundary{2,Upper}())
+
+    v = zeros(Float64, 4, 5)
+    v[:,5] = [1, 2, 3,4]
+    v[:,4] = [1, 2, 3,4]
+    v[:,3] = [4, 5, 6, 7]
+    v[:,2] = [7, 8, 9, 10]
+    v[:,1] = [10, 11, 12, 13]
+
+    @test e_w isa TensorMapping{T,1,2} where T
+    @test e_w'  isa TensorMapping{T,2,1} where T
+
+    @test domain_size(e_w) == (4,5)
+    @test domain_size(e_e) == (4,5)
+    @test domain_size(e_s) == (4,5)
+    @test domain_size(e_n) == (4,5)
+
+    @test range_size(e_w) == (5,)
+    @test range_size(e_e) == (5,)
+    @test range_size(e_s) == (4,)
+    @test range_size(e_n) == (4,)
+
+    I_w = [(Index{Lower}(1),),
+           (Index{Interior}(2),),
+           (Index{Interior}(3),),
+           (Index{Interior}(4),),
+           (Index{Upper}(5),)]
+    v_w = [10,7,4,1.0,1];
+    for i = 1:length(I_w)
+        @test (e_w*v)[I_w[i]...] == v_w[i];
+    end
+    @test e_w*v == [10,7,4,1.0,1]
+    @test e_e*v == [13,10,7,4,4.0]
+    @test e_s*v == [10,11,12,13.0]
+    @test e_n*v == [1,2,3,4.0]
+
+    g_x = [1,2,3,4.0]
+    g_y = [5,4,3,2,1.0]
+
+    G_w = zeros(Float64, (4,5))
+    G_w[1,:] = g_y
+
+    G_e = zeros(Float64, (4,5))
+    G_e[4,:] = g_y
+
+    G_s = zeros(Float64, (4,5))
+    G_s[:,1] = g_x
+
+    G_n = zeros(Float64, (4,5))
+    G_n[:,5] = g_x
+
+    @test e_w'*g_y == G_w
+    @test e_e'*g_y == G_e
+    @test e_s'*g_x == G_s
+    @test e_n'*g_x == G_n
+end
 #
 # @testset "NormalDerivative" begin
 #     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")