Mercurial > repos > public > sbplib_julia
diff SbpOperators/src/d2.jl @ 249:7cb4492ccd60 boundary_conditions
Refactor package SbpOperators
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Thu, 27 Jun 2019 14:18:48 +0200 |
parents | |
children | 396eadb652bd |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/SbpOperators/src/d2.jl Thu Jun 27 14:18:48 2019 +0200 @@ -0,0 +1,76 @@ +export D2, closureSize, readOperator, apply_e, apply_d, apply_e_T, apply_d_T + +@enum Parity begin + odd = -1 + even = 1 +end + +struct D2{T,N,M,K} <: ConstantStencilOperator + quadratureClosure::NTuple{M,T} + innerStencil::Stencil{T,N} + closureStencils::NTuple{M,Stencil{T,K}} + eClosure::Stencil{T,M} + dClosure::Stencil{T,M} + parity::Parity +end + +function closureSize(D::D2)::Int + return length(D.quadratureClosure) +end + +function apply_e_T(op::D2, v::AbstractVector, ::Type{Lower}) + @boundscheck if length(v) < closureSize(op) + throw(BoundsError()) + end + apply(op.eClosure,v,1) +end + +function apply_e_T(op::D2, v::AbstractVector, ::Type{Upper}) + @boundscheck if length(v) < closureSize(op) + throw(BoundsError()) + end + apply(flip(op.eClosure),v,length(v)) +end + + +function apply_e(op::D2, v::Number, N::Integer, i::Integer, ::Type{Lower}) + @boundscheck if !(0<length(i) <= N) + throw(BoundsError()) + end + op.eClosure[i-1]*v +end + +function apply_e(op::D2, v::Number, N::Integer, i::Integer, ::Type{Upper}) + @boundscheck if !(0<length(i) <= N) + throw(BoundsError()) + end + op.eClosure[N-i]*v +end + +function apply_d_T(op::D2, h_inv::Real, v::AbstractVector, ::Type{Lower}) + @boundscheck if length(v) < closureSize(op) + throw(BoundsError()) + end + h_inv*apply(op.dClosure,v,1) +end + +function apply_d_T(op::D2, h_inv::Real, v::AbstractVector, ::Type{Upper}) + @boundscheck if length(v) < closureSize(op) + throw(BoundsError()) + end + -h_inv*apply(flip(op.dClosure),v,length(v)) +end + +function apply_d(op::D2, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Lower}) + @boundscheck if !(0<length(i) <= N) + throw(BoundsError()) + end + h_inv*op.dClosure[i-1]*v +end + +function apply_d(op::D2, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Upper}) + @boundscheck if !(0<length(i) <= N) + throw(BoundsError()) + end + -h_inv*op.dClosure[N-i]*v +end