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