Mercurial > repos > public > sbplib_julia

comparison SbpOperators/src/quadrature/quadrature.jl @ 328:9cc5d1498b2d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactor 1D diagonal inner product in quadrature.jl to separate file. Write tests for quadratures. Clean up laplace and secondderivative
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 24 Sep 2020 22:31:48 +0200
parents 802edc9f252e
children
comparison
equal deleted inserted replaced
327:802edc9f252e 328:9cc5d1498b2d
1 # At the moment the grid property is used all over. It could possibly be removed if we implement all the 1D operators as TensorMappings 1 export Quadrature
2 """ 2 """
3 Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} 3 Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
4 4
5 Implements the quadrature operator `Q` of Dim dimension as a TensorMapping 5 Implements the quadrature operator `Q` of Dim dimension as a TensorMapping
6 The multi-dimensional tensor operator consists of a tuple of 1D DiagonalNorm H 6 The multi-dimensional tensor operator consists of a tuple of 1D DiagonalInnerProduct H
7 tensor operators. 7 tensor operators.
8 """ 8 """
9 export Quadrature 9 struct Quadrature{Dim,T<:Real,M} <: TensorOperator{T,Dim}
10 struct Quadrature{Dim,T<:Real,N,M} <: TensorOperator{T,Dim} 10 H::NTuple{Dim,DiagonalInnerProduct{T,M}}
11 H::NTuple{Dim,DiagonalNorm{T,N,M}}
12 end 11 end
13 12
14 LazyTensors.domain_size(Q::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size 13 LazyTensors.domain_size(Q::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where {Dim} = range_size
15 14
16 function LazyTensors.apply(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim} 15 function LazyTensors.apply(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim}
17 error("not implemented") 16 error("not implemented")
18 end 17 end
19 18
20 LazyTensors.apply_transpose(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {Dim,T} = LazyTensors.apply(Q,v,I) 19 function LazyTensors.apply(Q::Quadrature{1,T}, v::AbstractVector{T}, I::Index) where T
21
22 @inline function LazyTensors.apply(Q::Quadrature{1,T}, v::AbstractVector{T}, I::Index) where T
23 @inbounds q = apply(Q.H[1], v , I) 20 @inbounds q = apply(Q.H[1], v , I)
24 return q 21 return q
25 end 22 end
26 23
27 @inline function LazyTensors.apply(Q::Quadrature{2,T}, v::AbstractArray{T,2}, I::Index, J::Index) where T 24 function LazyTensors.apply(Q::Quadrature{2,T}, v::AbstractArray{T,2}, I::Index, J::Index) where T
28 # Quadrature in x direction 25 # Quadrature in x direction
29 @inbounds vx = view(v, :, Int(J)) 26 @inbounds vx = view(v, :, Int(J))
30 @inbounds qx = apply(Q.H[1], vx , I) 27 @inbounds qx = apply(Q.H[1], vx , I)
31 # Quadrature in y-direction 28 # Quadrature in y-direction
32 @inbounds vy = view(v, Int(I), :) 29 @inbounds vy = view(v, Int(I), :)
33 @inbounds qy = apply(Q.H[2], vy, J) 30 @inbounds qy = apply(Q.H[2], vy, J)
34 return qx*qy 31 return qx*qy
35 end 32 end
36 33
37 """ 34 LazyTensors.apply_transpose(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {Dim,T} = LazyTensors.apply(Q,v,I...)
38 DiagonalNorm{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
39
40 Implements the diagnoal norm operator `H` of Dim dimension as a TensorMapping
41 """
42 export DiagonalNorm, closuresize, LazyTensors.apply
43 struct DiagonalNorm{T<:Real,N,M} <: TensorOperator{T,1}
44 h::T # The grid spacing could be included in the stencil already. Preferable?
45 closure::NTuple{M,T}
46 #TODO: Write a nice constructor
47 end
48
49 @inline function LazyTensors.apply(H::DiagonalNorm{T}, v::AbstractVector{T}, I::Index) where T
50 return @inbounds apply(H, v, I)
51 end
52
53 LazyTensors.apply_transpose(H::Quadrature{Dim,T}, v::AbstractArray{T,2}, I::Index) where T = LazyTensors.apply(H,v,I)
54
55 @inline LazyTensors.apply(H::DiagonalNorm, v::AbstractVector{T}, I::Index{Lower}) where T
56 return @inbounds H.h*H.closure[Int(I)]*v[Int(I)]
57 end
58 @inline LazyTensors.apply(H::DiagonalNorm,v::AbstractVector{T}, I::Index{Upper}) where T
59 N = length(v);
60 return @inbounds H.h*H.closure[N-Int(I)+1]v[Int(I)]
61 end
62
63 @inline LazyTensors.apply(H::DiagonalNorm, v::AbstractVector{T}, I::Index{Interior}) where T
64 return @inbounds H.h*v[Int(I)]
65 end
66
67 function LazyTensors.apply(H::DiagonalNorm, v::AbstractVector{T}, index::Index{Unknown}) where T
68 N = length(v);
69 r = getregion(Int(index), closuresize(H), N)
70 i = Index(Int(index), r)
71 return LazyTensors.apply(H, v, i)
72 end
73
74 function closuresize(H::DiagonalNorm{T<:Real,N,M}) where {T,N,M}
75 return M
76 end