Mercurial > repos > public > sbplib_julia

comparison SbpOperators/src/laplace/secondderivative.jl @ 296:c5edb432954e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 22 Jun 2020 22:19:10 +0200
parents f63232aeb1c6
children 27a0bca5e1f2
comparison
equal deleted inserted replaced
295:4735abcf5d42 296:c5edb432954e
1 """
2 SecondDerivative{T<:Real,N,M,K} <: TensorOperator{T,1}
3 Implements the Laplace tensor operator `L` with constant grid spacing and coefficients
4 in 1D dimension
5 """
6 struct SecondDerivative{T<:Real,N,M,K} <: TensorOperator{T,1}
7 h_inv::T # The grid spacing could be included in the stencil already. Preferable?
8 innerStencil::Stencil{T,N}
9 closureStencils::NTuple{M,Stencil{T,K}}
10 parity::Parity
11 #TODO: Write a nice constructor
12 end
13
14 @enum Parity begin
15 odd = -1
16 even = 1
17 end
18
19 LazyTensors.domain_size(D2::SecondDerivative, range_size::NTuple{1,Integer}) = range_size
20
21 function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
22 return apply(D2, v, I[1])
23 end
24
25 # Apply for different regions Lower/Interior/Upper or Unknown region
26 @inline function LazyTensors.apply(D2::SecondDerivative, v::AbstractVector, i::Index{Lower})
27 return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.closureStencils[Int(i)], v, Int(i))
28 end
29
30 @inline function LazyTensors.apply(D2::SecondDerivative, v::AbstractVector, i::Index{Interior})
31 return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.innerStencil, v, Int(i))
32 end
33
34 @inline function LazyTensors.apply(D2::SecondDerivative, v::AbstractVector, i::Index{Upper})
35 N = length(v) # TODO: Use domain_size here instead?
36 return @inbounds D2.h_inv*D2.h_inv*Int(D2.parity)*apply_stencil_backwards(D2.closureStencils[N-Int(i)+1], v, Int(i))
37 end
38
39 @inline function LazyTensors.apply(D2::SecondDerivative, v::AbstractVector, index::Index{Unknown})
40 N = length(v) # TODO: Use domain_size here instead?
41 r = getregion(Int(index), closuresize(L), N)
42 i = Index(Int(index), r)
43 return apply(D2, v, i)
44 end
45
46 function closuresize(D2::SecondDerivative{T<:Real,N,M,K}) where T,N,M,K
47 return M
48 end