Mercurial > repos > public > sbplib_julia

comparison SbpOperators/src/InverseQuadrature.jl @ 302:6fa2ba769ae3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create 1D tensor mapping for inverse diagonal norm, and make the multi-dimensional inverse quadrature use those. Move InverseQudrature from laplace.jl into InverseQuadrature.jl
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 23 Jun 2020 18:56:59 +0200
parents
children bd09d67ebb22
comparison
equal deleted inserted replaced
301:417b767c847f 302:6fa2ba769ae3
1 """
2 Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
3
4 Implements the inverse quadrature operator `Qi` of Dim dimension as a TensorOperator
5 The multi-dimensional tensor operator consists of a tuple of 1D InverseDiagonalNorm
6 tensor operators.
7 """
8 struct Quadrature{Dim,T<:Real,N,M} <: TensorOperator{T,Dim}
9 Hi::NTuple{Dim,InverseDiagonalNorm{T,N,M}}
10 end
11 export Quadrature
12
13 LazyTensors.domain_size(Qi::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
14
15 function LazyTensors.apply(Qi::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim}
16 error("not implemented")
17 end
18
19 LazyTensors.apply_transpose(Qi::Quadrature{Dim,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where {Dim,T} = LazyTensors.apply(Q,v,I)
20
21 @inline function LazyTensors.apply(Qi::Quadrature{1,T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
22 @inbounds q = apply(Qi.Hi[1], v , I[1])
23 return q
24 end
25
26 @inline function LazyTensors.apply(Qi::Quadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
27 # Quadrature in x direction
28 @inbounds vx = view(v, :, Int(I[2]))
29 @inbounds qx_inv = apply(Qi.Hi[1], vx , I[1])
30 # Quadrature in y-direction
31 @inbounds vy = view(v, Int(I[1]), :)
32 @inbounds qy_inv = apply(Qi.Hi[2], vy, I[2])
33 return qx_inv*qy_inv
34 end
35
36 """
37 Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
38
39 Implements the quadrature operator `Hi` of Dim dimension as a TensorMapping
40 """
41 struct InverseDiagonalNorm{T<:Real,N,M} <: TensorOperator{T,1}
42 h_inv::T # The reciprocl grid spacing could be included in the stencil already. Preferable?
43 closure::NTuple{M,T}
44 #TODO: Write a nice constructor
45 end
46
47 @inline function LazyTensors.apply(Hi::InverseDiagonalNorm{T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
48 return @inbounds apply(Hi, v, I[1])
49 end
50
51 LazyTensors.apply_transpose(Hi::Quadrature{Dim,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T = LazyTensors.apply(Hi,v,I)
52
53 @inline LazyTensors.apply(Hi::InverseDiagonalNorm, v::AbstractVector{T}, i::Index{Lower}) where T
54 return @inbounds Hi.h_inv*Hi.closure[Int(i)]*v[Int(i)]
55 end
56 @inline LazyTensors.apply(Hi::InverseDiagonalNorm,v::AbstractVector{T}, i::Index{Upper}) where T
57 N = length(v);
58 return @inbounds Hi.h_inv*Hi.closure[N-Int(i)+1]v[Int(i)]
59 end
60
61 @inline LazyTensors.apply(Hi::InverseDiagonalNorm, v::AbstractVector{T}, i::Index{Interior}) where T
62 return @inbounds Hi.h_inv*v[Int(i)]
63 end
64
65 function LazyTensors.apply(Hi::InverseDiagonalNorm, v::AbstractVector{T}, index::Index{Unknown}) where T
66 N = length(v);
67 r = getregion(Int(index), closuresize(Hi), N)
68 i = Index(Int(index), r)
69 return LazyTensors.apply(Hi, v, i)
70 end
71 export LazyTensors.apply
72
73 function closuresize(Hi::InverseDiagonalNorm{T<:Real,N,M}) where {T,N,M}
74 return M
75 end