diff src/SbpOperators/volumeops/quadratures/inverse_quadrature.jl @ 642:f4a16b403487 feature/volume_and_boundary_operators

Implement the inverse quadrature operator as a volume operator and update tests.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 04 Jan 2021 17:17:40 +0100
parents src/SbpOperators/volumeops/quadratures/inverse_diagonal_quadrature.jl@a1dfaf305f41
children e14627e79a54
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/quadratures/inverse_quadrature.jl	Mon Jan 04 17:17:40 2021 +0100
@@ -0,0 +1,41 @@
+
+"""
+    InverseQuadrature(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)
+
+Creates the inverse `H⁻¹` of the quadrature operator as a `TensorMapping`
+
+The inverse quadrature approximates the integral operator on the grid using
+`inv_inner_stencil` in the interior and a set of stencils `inv_closure_stencils`
+for the points in the closure regions.
+
+On a one-dimensional `grid`, `H⁻¹` is a `VolumeOperator`. On a multi-dimensional
+`grid`, `H` is the outer product of the 1-dimensional inverse quadrature operators in
+each coordinate direction. Also see the documentation of
+`SbpOperators.volume_operator(...)` for more details.
+"""
+function InverseQuadrature(grid::EquidistantGrid{Dim}, inv_inner_stencil, inv_closure_stencils) where Dim
+    h⁻¹ = inverse_spacing(grid)
+    H⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[1]),scale.(inv_closure_stencils,h⁻¹[1]),even,1)
+    for i ∈ 2:Dim
+        Hᵢ⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[i]),scale.(inv_closure_stencils,h⁻¹[i]),even,i)
+        H⁻¹ = H⁻¹∘Hᵢ⁻¹
+    end
+    return H⁻¹
+end
+export InverseQuadrature
+
+"""
+    InverseDiagonalQuadrature(grid::EquidistantGrid, closure_stencils)
+
+Creates the inverse of the diagonal quadrature operator defined by the inner stencil
+1/h and a set of 1-element closure stencils in `closure_stencils`. Note that
+the closure stencils are those of the quadrature operator (and not the inverse).
+"""
+function InverseDiagonalQuadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {T,M}
+    inv_inner_stencil = Stencil(Tuple{T}(1),center=1)
+    inv_closure_stencils = reciprocal_stencil.(closure_stencils)
+    return InverseQuadrature(grid, inv_inner_stencil, inv_closure_stencils)
+end
+export InverseDiagonalQuadrature
+
+reciprocal_stencil(s::Stencil{T}) where T = Stencil(s.range,one(T)./s.weights)