Mercurial > repos > public > sbplib_julia
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)