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