Mercurial > repos > public > sbplib_julia
diff src/SbpOperators/quadrature/inverse_diagonal_quadrature.jl @ 504:21fba50cb5b0 feature/quadrature_as_outer_product
Use LazyOuterProduct to construct multi-dimensional quadratures. This change allwed to:
- Replace the types Quadrature and InverseQuadrature by functions returning outer products of the 1D operators.
- Avoid convoluted naming of types. DiagonalInnerProduct is now renamed to DiagonalQuadrature, similarly for InverseDiagonalInnerProduct.
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Sat, 07 Nov 2020 13:07:31 +0100 |
parents | src/SbpOperators/quadrature/inverse_diagonal_inner_product.jl@0d93d406c222 |
children | c2f991b819fc |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/quadrature/inverse_diagonal_quadrature.jl Sat Nov 07 13:07:31 2020 +0100 @@ -0,0 +1,64 @@ +""" +inverse_diagonal_quadrature(g,quadrature_closure) + +Constructs the diagonal quadrature inverse operator `Qi` on a grid of `Dim` dimensions as +a `TensorMapping`. The one-dimensional operator is a InverseDiagonalQuadrature, while +the multi-dimensional operator is the outer-product of the one-dimensional operators +in each coordinate direction. +""" +function inverse_diagonal_quadrature(g::EquidistantGrid{Dim}, quadrature_closure) where Dim + Hi = InverseDiagonalQuadrature(restrict(g,1), quadrature_closure) + for i ∈ 2:Dim + Hi = Hi⊗InverseDiagonalQuadrature(restrict(g,i), quadrature_closure) + end + return Hi +end +export inverse_diagonal_quadrature + + +""" + InverseDiagonalQuadrature{Dim,T<:Real,M} <: TensorMapping{T,1,1} + +Implements the inverse diagonal inner product operator `Hi` of as a 1D TensorOperator +""" +struct InverseDiagonalQuadrature{T<:Real,M} <: TensorMapping{T,1,1} + h_inv::T + closure::NTuple{M,T} + size::Tuple{Int} +end +export InverseDiagonalQuadrature + +function InverseDiagonalQuadrature(g::EquidistantGrid{1}, quadrature_closure) + return InverseDiagonalQuadrature(inverse_spacing(g)[1], 1 ./ quadrature_closure, size(g)) +end + + +LazyTensors.range_size(Hi::InverseDiagonalQuadrature) = Hi.size +LazyTensors.domain_size(Hi::InverseDiagonalQuadrature) = Hi.size + + +function LazyTensors.apply(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, I::Index{Lower}) where T + return @inbounds Hi.h_inv*Hi.closure[Int(I)]*v[Int(I)] +end + +function LazyTensors.apply(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, I::Index{Upper}) where T + N = length(v); + return @inbounds Hi.h_inv*Hi.closure[N-Int(I)+1]*v[Int(I)] +end + +function LazyTensors.apply(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, I::Index{Interior}) where T + return @inbounds Hi.h_inv*v[Int(I)] +end + +function LazyTensors.apply(Hi::InverseDiagonalQuadrature, v::AbstractVector{T}, index::Index{Unknown}) where T + N = length(v); + r = getregion(Int(index), closuresize(Hi), N) + i = Index(Int(index), r) + return LazyTensors.apply(Hi, v, i) +end + +LazyTensors.apply_transpose(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, I::Index) where T = LazyTensors.apply(Hi,v,I) + + +closuresize(Hi::InverseDiagonalQuadrature{T,M}) where {T,M} = M +export closuresize