Mercurial > repos > public > sbplib_julia
comparison src/SbpOperators/volumeops/volume_operator.jl @ 651:67639b1c99ea
Merged feature/volume_and_boundary_operators
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Wed, 20 Jan 2021 17:52:55 +0100 |
| parents | 9f27f451d0a0 |
| children | 86bb3606b215 |
comparison
equal
deleted
inserted
replaced
| 615:52749b687a67 | 651:67639b1c99ea |
|---|---|
| 1 """ | |
| 2 volume_operator(grid,inner_stencil,closure_stencils,parity,direction) | |
| 3 Creates a volume operator on a `Dim`-dimensional grid acting along the | |
| 4 specified coordinate `direction`. The action of the operator is determined by the | |
| 5 stencils `inner_stencil` and `closure_stencils`. | |
| 6 When `Dim=1`, the corresponding `VolumeOperator` tensor mapping is returned. | |
| 7 When `Dim>1`, the `VolumeOperator` `op` is inflated by the outer product | |
| 8 of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`, | |
| 9 the boundary restriction operator in the y-direction direction is `Ix⊗op⊗Iz`. | |
| 10 """ | |
| 11 function volume_operator(grid::EquidistantGrid{Dim,T}, inner_stencil::Stencil{T}, closure_stencils::NTuple{M,Stencil{T}}, parity, direction) where {Dim,T,M} | |
| 12 #TODO: Check that direction <= Dim? | |
| 13 | |
| 14 # Create 1D volume operator in along coordinate direction | |
| 15 op = VolumeOperator(restrict(grid, direction), inner_stencil, closure_stencils, parity) | |
| 16 # Create 1D IdentityMappings for each coordinate direction | |
| 17 one_d_grids = restrict.(Ref(grid), Tuple(1:Dim)) | |
| 18 Is = IdentityMapping{T}.(size.(one_d_grids)) | |
| 19 # Formulate the correct outer product sequence of the identity mappings and | |
| 20 # the volume operator | |
| 21 parts = Base.setindex(Is, op, direction) | |
| 22 return foldl(⊗, parts) | |
| 23 end | |
| 24 | |
| 25 """ | |
| 26 VolumeOperator{T,N,M,K} <: TensorOperator{T,1} | |
| 27 Implements a one-dimensional constant coefficients volume operator | |
| 28 """ | |
| 29 struct VolumeOperator{T,N,M,K} <: TensorMapping{T,1,1} | |
| 30 inner_stencil::Stencil{T,N} | |
| 31 closure_stencils::NTuple{M,Stencil{T,K}} | |
| 32 size::NTuple{1,Int} | |
| 33 parity::Parity | |
| 34 end | |
| 35 | |
| 36 function VolumeOperator(grid::EquidistantGrid{1}, inner_stencil, closure_stencils, parity) | |
| 37 return VolumeOperator(inner_stencil, closure_stencils, size(grid), parity) | |
| 38 end | |
| 39 | |
| 40 closure_size(::VolumeOperator{T,N,M}) where {T,N,M} = M | |
| 41 | |
| 42 LazyTensors.range_size(op::VolumeOperator) = op.size | |
| 43 LazyTensors.domain_size(op::VolumeOperator) = op.size | |
| 44 | |
| 45 function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Lower}) where T | |
| 46 return @inbounds apply_stencil(op.closure_stencils[Int(i)], v, Int(i)) | |
| 47 end | |
| 48 | |
| 49 function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Interior}) where T | |
| 50 return apply_stencil(op.inner_stencil, v, Int(i)) | |
| 51 end | |
| 52 | |
| 53 function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Upper}) where T | |
| 54 return @inbounds Int(op.parity)*apply_stencil_backwards(op.closure_stencils[op.size[1]-Int(i)+1], v, Int(i)) | |
| 55 end | |
| 56 | |
| 57 function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i) where T | |
| 58 r = getregion(i, closure_size(op), op.size[1]) | |
| 59 return LazyTensors.apply(op, v, Index(i, r)) | |
| 60 end |