Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/quadrature/diagonal_quadrature.jl @ 634:fb5ac62563aa feature/volume_and_boundary_operators

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Integrate feature/quadrature_as_outer_product into branch, before closing feature/quadrature_as_outer_product. (It is now obsolete apart from tests)
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 01 Jan 2021 16:39:57 +0100
parents 04d7b4eb63ef
children
comparison
equal deleted inserted replaced
632:bf8b66c596f7 634:fb5ac62563aa
1 """
2 diagonal_quadrature(g,quadrature_closure)
3
4 Constructs the diagonal quadrature operator `H` on a grid of `Dim` dimensions as
5 a `TensorMapping`. The one-dimensional operator is a `DiagonalQuadrature`, while
6 the multi-dimensional operator is the outer-product of the
7 one-dimensional operators in each coordinate direction.
8 """
9 function diagonal_quadrature(g::EquidistantGrid{Dim}, quadrature_closure) where Dim
10 H = DiagonalQuadrature(restrict(g,1), quadrature_closure)
11 for i ∈ 2:Dim
12 H = H⊗DiagonalQuadrature(restrict(g,i), quadrature_closure)
13 end
14 return H
15 end
16 export diagonal_quadrature
17
18 """
19 DiagonalQuadrature{T,M} <: TensorMapping{T,1,1}
20
21 Implements the one-dimensional diagonal quadrature operator as a `TensorMapping`
22 The quadrature is defined by the quadrature interval length `h`, the quadrature
23 closure weights `closure` and the number of quadrature intervals `size`. The
24 interior stencil has the weight 1.
25 """
26 struct DiagonalQuadrature{T,M} <: TensorMapping{T,1,1}
27 h::T
28 closure::NTuple{M,T}
29 size::Tuple{Int}
30 end
31 export DiagonalQuadrature
32
33 """
34 DiagonalQuadrature(g, quadrature_closure)
35
36 Constructs the `DiagonalQuadrature` on the `EquidistantGrid` `g` with
37 closure given by `quadrature_closure`.
38 """
39 function DiagonalQuadrature(g::EquidistantGrid{1}, quadrature_closure)
40 return DiagonalQuadrature(spacing(g)[1], quadrature_closure, size(g))
41 end
42
43 """
44 range_size(H::DiagonalQuadrature)
45
46 The size of an object in the range of `H`
47 """
48 LazyTensors.range_size(H::DiagonalQuadrature) = H.size
49
50 """
51 domain_size(H::DiagonalQuadrature)
52
53 The size of an object in the domain of `H`
54 """
55 LazyTensors.domain_size(H::DiagonalQuadrature) = H.size
56
57 """
58 apply(H::DiagonalQuadrature{T}, v::AbstractVector{T}, i) where T
59 Implements the application `(H*v)[i]` an `Index{R}` where `R` is one of the regions
60 `Lower`,`Interior`,`Upper`. If `i` is another type of index (e.g an `Int`) it will first
61 be converted to an `Index{R}`.
62 """
63 function LazyTensors.apply(H::DiagonalQuadrature{T}, v::AbstractVector{T}, i::Index{Lower}) where T
64 return @inbounds H.h*H.closure[Int(i)]*v[Int(i)]
65 end
66
67 function LazyTensors.apply(H::DiagonalQuadrature{T},v::AbstractVector{T}, i::Index{Upper}) where T
68 N = length(v); #TODO: Use dim_size here?
69 return @inbounds H.h*H.closure[N-Int(i)+1]*v[Int(i)]
70 end
71
72 function LazyTensors.apply(H::DiagonalQuadrature{T}, v::AbstractVector{T}, i::Index{Interior}) where T
73 return @inbounds H.h*v[Int(i)]
74 end
75
76 function LazyTensors.apply(H::DiagonalQuadrature{T}, v::AbstractVector{T}, i) where T
77 N = length(v); #TODO: Use dim_size here?
78 r = getregion(i, closure_size(H), N)
79 return LazyTensors.apply(H, v, Index(i, r))
80 end
81
82 """
83 apply(H::DiagonalQuadrature{T}, v::AbstractVector{T}, I::Index) where T
84 Implements the application (H'*v)[I]. The operator is self-adjoint.
85 """
86 LazyTensors.apply_transpose(H::DiagonalQuadrature{T}, v::AbstractVector{T}, i) where T = LazyTensors.apply(H,v,i)
87
88 """
89 closure_size(H)
90 Returns the size of the closure stencil of a DiagonalQuadrature `H`.
91 """
92 closure_size(H::DiagonalQuadrature{T,M}) where {T,M} = M