changeset 515:d55008f5e2f3 feature/boundary_ops

Fix the range of the BoundaryRestriction tensor mapping (the range is zero for the 1D operator). Add some documentation and todos.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 24 Nov 2020 15:28:10 +0100
parents 14e722e8607d
children f932af8c1e56
files src/SbpOperators/boundaryops/boundary_restriction.jl test/testSbpOperators.jl
diffstat 2 files changed, 23 insertions(+), 15 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Mon Nov 23 21:19:08 2020 +0100
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Tue Nov 24 15:28:10 2020 +0100
@@ -1,3 +1,8 @@
+"""
+    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())
@@ -16,15 +21,14 @@
     I = IdentityMapping{T}(size(restrict(grid,1)))
     return I⊗e
 end
-
 export boundary_restriction
 
 """
-    BoundaryRestriction{T,N,R} <: TensorMapping{T,1,1}
+    BoundaryRestriction{T,N,R} <: TensorMapping{T,0,1}
 
 Implements the boundary operator `e` as a TensorMapping
 """
-struct BoundaryRestriction{T,M,R<:Region} <: TensorMapping{T,1,1}
+struct BoundaryRestriction{T,M,R<:Region} <: TensorMapping{T,0,1}
     stencil::Stencil{T,M}
     size::NTuple{1,Int}
 end
@@ -34,9 +38,13 @@
     return BoundaryRestriction{T,M,typeof(region)}(closureStencil,size(grid))
 end
 
-LazyTensors.range_size(e::BoundaryRestriction) = (1,)
+LazyTensors.range_size(e::BoundaryRestriction) = (0,)
 LazyTensors.domain_size(e::BoundaryRestriction) = e.size
 
+# TODO: Currently not working.
+# We need to handle getindex for LazyTensorMappingApplication such that we pass more #indices than the
+# range size of the TensorMapping. Or we need to be able to handle the case where we dont pass any index, for
+# 0-dimensional tensormappings.
 " 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}, i::Index{Lower}) where {T,M}
     @boundscheck if Int(i)!=1
@@ -55,18 +63,18 @@
 
 " 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::AbstractVector{T}, i) where {T,M}
+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[1]
+    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::AbstractVector{T}, i) where {T,M}
+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[1]
+    return e.stencil[e.size[1] - Int(i)]*v[]
 end
--- a/test/testSbpOperators.jl	Mon Nov 23 21:19:08 2020 +0100
+++ b/test/testSbpOperators.jl	Tue Nov 24 15:28:10 2020 +0100
@@ -182,13 +182,13 @@
     e_r = BoundaryRestriction(g,op.eClosure,Upper())
 
     v = evalOn(g,x->1+x^2)
-    u = [3.124] #How to handle scalars having to be arrays? It's kind of ugly.
+    u = fill(3.124)
 
     e_l*v isa LazyTensorMappingApplication
-    @test (e_l*v)[Index{Lower}(1)] == v[1]
-    @test (e_r*v)[Index{Upper}(4)] == v[end]
-    @test e_l'*u == [u[1], 0, 0, 0]
-    @test e_r'*u == [0, 0, 0, u[1]]
+    @test_broken (e_l*v)[Index{Lower}(1)] == v[1]
+    @test_broken (e_r*v)[Index{Upper}(4)] == 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)]
 
@@ -206,8 +206,8 @@
     v[:,2] = [7, 8, 9, 10]
     v[:,1] = [10, 11, 12, 13]
 
-    @test_broken e_w isa TensorMapping{T,1,2} where T
-    @test_broken e_w'  isa TensorMapping{T,2,1} where T
+    @test e_w isa TensorMapping{T,1,2} where T
+    @test e_w'  isa TensorMapping{T,2,1} where T