Mercurial > repos > public > sbplib_julia
diff SbpOperators/src/laplace/secondderivative.jl @ 293:f63232aeb1c6
Move laplace.jl from DiffOps to SbpOperators. Rename constandlaplace to secondderivative
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Mon, 22 Jun 2020 22:08:56 +0200 |
parents | SbpOperators/src/constantlaplace.jl@3747e5636eef |
children | 27a0bca5e1f2 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/SbpOperators/src/laplace/secondderivative.jl Mon Jun 22 22:08:56 2020 +0200 @@ -0,0 +1,48 @@ +""" + SecondDerivative{T<:Real,N,M,K} <: TensorOperator{T,1} +Implements the Laplace tensor operator `L` with constant grid spacing and coefficients +in 1D dimension +""" +struct SecondDerivative{T<:Real,N,M,K} <: TensorOperator{T,1} + h_inv::T # The grid spacing could be included in the stencil already. Preferable? + innerStencil::Stencil{T,N} + closureStencils::NTuple{M,Stencil{T,K}} + parity::Parity + #TODO: Write a nice constructor +end + +@enum Parity begin + odd = -1 + even = 1 +end + +LazyTensors.domain_size(D2::SecondDerivative, range_size::NTuple{1,Integer}) = range_size + +function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T + return apply(D2, v, I[1]) +end + +# Apply for different regions Lower/Interior/Upper or Unknown region +@inline function LazyTensors.apply(D2::SecondDerivative, v::AbstractVector, i::Index{Lower}) + return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.closureStencils[Int(i)], v, Int(i)) +end + +@inline function LazyTensors.apply(D2::SecondDerivative, v::AbstractVector, i::Index{Interior}) + return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.innerStencil, v, Int(i)) +end + +@inline function LazyTensors.apply(D2::SecondDerivative, v::AbstractVector, i::Index{Upper}) + N = length(v) # TODO: Use domain_size here instead? + return @inbounds D2.h_inv*D2.h_inv*Int(D2.parity)*apply_stencil_backwards(D2.closureStencils[N-Int(i)+1], v, Int(i)) +end + +@inline function LazyTensors.apply(D2::SecondDerivative, v::AbstractVector, index::Index{Unknown}) + N = length(v) # TODO: Use domain_size here instead? + r = getregion(Int(index), closuresize(L), N) + i = Index(Int(index), r) + return apply(D2, v, i) +end + +function closuresize(D2::SecondDerivative{T<:Real,N,M,K}) where T,N,M,K + return M +end