comparison src/SbpOperators/boundaryops/boundary_restriction.jl @ 511:f5f3b832f9c4 feature/boundary_ops

Fix boundschecking for BoundaryRestriction
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 23 Nov 2020 13:09:57 +0100
parents db64cfe4d9de
children 5a8cfcc0765d
comparison
equal deleted inserted replaced
510:db64cfe4d9de 511:f5f3b832f9c4
2 BoundaryRestriction{T,N,R} <: TensorMapping{T,1,1} 2 BoundaryRestriction{T,N,R} <: TensorMapping{T,1,1}
3 3
4 Implements the boundary operator `e` as a TensorMapping 4 Implements the boundary operator `e` as a TensorMapping
5 """ 5 """
6 struct BoundaryRestriction{T,M,R<:Region} <: TensorMapping{T,1,1} 6 struct BoundaryRestriction{T,M,R<:Region} <: TensorMapping{T,1,1}
7 closureStencil::Stencil{T,M} 7 stencil::Stencil{T,M}
8 size::NTuple{1,Int} 8 size::NTuple{1,Int}
9 end 9 end
10 export BoundaryRestriction 10 export BoundaryRestriction
11 11
12 12 function BoundaryRestriction(grid::EquidistantGrid{1,T}, closureStencil::Stencil{T,M}, region::Region) where {T,M}
13 function BoundaryRestriction(grid::EquidistantGrid{1,T}, closureStencil::Stencil{T,M},region::Region) where {T,M} 13 return BoundaryRestriction{T,M,typeof(region)}(stencil,size(grid))
14 return BoundaryRestriction{T,M,typeof(region)}(closureStencil,size(grid))
15 end 14 end
16 15
17 LazyTensors.range_size(e::BoundaryRestriction) = (1,) 16 LazyTensors.range_size(e::BoundaryRestriction) = (1,)
18 LazyTensors.domain_size(e::BoundaryRestriction) = e.size 17 LazyTensors.domain_size(e::BoundaryRestriction) = e.size
19 18
20 " Restricts a grid function v on a grid of size (m,) to the scalar element v[1]" 19 " Restricts a grid function v on a grid of size m to the scalar element v[1]"
21 function LazyTensors.apply(e::BoundaryRestriction{T,M,Lower}, v::AbstractVector{T}, i) where {T,M} 20 function LazyTensors.apply(e::BoundaryRestriction{T,M,Lower}, v::AbstractVector{T}, i::Index{Lower}) where {T,M}
22 # TODO: Currently closuresize = 4 for the stencil, 21 @boundscheck if Int(i)!=1
23 # but we should only be able to apply this for v[1] since the result is scalar.
24 @boundscheck if !(0 < Int(i) <= closuresize(e))
25 throw(BoundsError()) 22 throw(BoundsError())
26 end 23 end
27 apply_stencil(e.closureStencil,v,Int(i)) 24 apply_stencil(e.stencil,v,Int(i))
28 end 25 end
29 26
30 " Restricts a grid function v on a grid of size (m,) to the scalar element v[m]" 27 " Restricts a grid function v on a grid of size m to the scalar element v[m]"
31 function LazyTensors.apply(e::BoundaryRestriction{T,M,Upper}, v::AbstractVector{T}, i) where {T,M} 28 function LazyTensors.apply(e::BoundaryRestriction{T,M,Upper}, v::AbstractVector{T}, i::Index{Upper}) where {T,M}
32 # TODO: Currently closuresize = 4 for the stencil, 29 @boundscheck if Int(i) != e.size[1]
33 # but we should only be able to apply this for v[1] since the result is scalar.
34 @boundscheck if !(e.size[1] - closuresize(e) < Int(i) <= e.size[1])
35 throw(BoundsError()) 30 throw(BoundsError())
36 end 31 end
37 apply_stencil_backwards(e.closureStencil,v,Int(i)) 32 apply_stencil_backwards(e.stencil,v,Int(i))
38 end 33 end
39 34
40 " Transpose of a restriction is an inflation or prolongation. 35 " Transpose of a restriction is an inflation or prolongation.
41 Inflates the scalar (1-element) vector to a vector of size of the grid" 36 Inflates the scalar (1-element) vector to a vector of size of the grid"
42 function LazyTensors.apply_transpose(e::BoundaryRestriction{T,M,Lower}, v::AbstractVector{T}, i) where {T,M} 37 function LazyTensors.apply_transpose(e::BoundaryRestriction{T,M,Lower}, v::AbstractVector{T}, i) where {T,M}
43 @boundscheck if !(0 < Int(i) <= e.size[1]) 38 @boundscheck if !(0 < Int(i) <= e.size[1])
44 throw(BoundsError()) 39 throw(BoundsError())
45 end 40 end
46 return e.closureStencil[Int(i)-1]*v[1] 41 return e.stencil[Int(i)-1]*v[1]
47 end 42 end
48
49 closuresize(e::BoundaryRestriction{T,M}) where {T,M} = M