diff SbpOperators/src/stencil.jl @ 219:69a6049e14d9 package_refactor

Create package SbpOperators
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 26 Jun 2019 13:12:38 +0200
parents stencil.jl@6b6d921e8f05
children ccef055233a2
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/stencil.jl	Wed Jun 26 13:12:38 2019 +0200
@@ -0,0 +1,38 @@
+struct Stencil{T<:Real,N}
+    range::Tuple{Int,Int}
+    weights::NTuple{N,T}
+
+    function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T <: Real, N}
+        @assert range[2]-range[1]+1 == N
+        new{T,N}(range,weights)
+    end
+end
+
+function flip(s::Stencil)
+    range = (-s.range[2], -s.range[1])
+    return Stencil(range, reverse(s.weights))
+end
+
+# Provides index into the Stencil based on offset for the root element
+@inline function Base.getindex(s::Stencil, i::Int)
+    @boundscheck if i < s.range[1] || s.range[2] < i
+        return eltype(s.weights)(0)
+    end
+    return s.weights[1 + i - s.range[1]]
+end
+
+Base.@propagate_inbounds @inline function apply(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
+    w = s.weights[1]*v[i + s.range[1]]
+    @simd for k ∈ 2:N
+        w += s.weights[k]*v[i + s.range[1] + k-1]
+    end
+    return w
+end
+
+Base.@propagate_inbounds @inline function apply_backwards(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
+    w = s.weights[N]*v[i - s.range[2]]
+    @simd for k ∈ N-1:-1:1
+        w += s.weights[k]*v[i - s.range[1] - k + 1]
+    end
+    return w
+end