First attempt at linear bc/extrap for B2

Author: Tomas Lycken

src/Interpolations.jl

@@ -12,10 +12,12 @@ export
     ExtrapError,
     ExtrapNaN,
     ExtrapConstant,
+    ExtrapLinear,
     OnCell,
     OnGrid,
     ExtendInner,
-    Flat
+    Flat,
+    LinearBC
 
 abstract Degree{N}
 
@@ -27,6 +29,7 @@ abstract BoundaryCondition
 type None <: BoundaryCondition end
 type ExtendInner <: BoundaryCondition end
 type Flat <: BoundaryCondition end
+type LinearBC <: BoundaryCondition end
 
 abstract InterpolationType{D<:Degree,BC<:BoundaryCondition,GR<:GridRepresentation}
 
@@ -66,9 +69,19 @@ include("quadratic.jl")
 
 
 # This creates getindex methods for all supported combinations
-for IT in (Constant{OnCell},Linear{OnGrid},Quadratic{ExtendInner,OnCell},Quadratic{Flat,OnCell})
-    for EB in (ExtrapError,ExtrapNaN,ExtrapConstant)
-
+for IT in (
+        Constant{OnCell},
+        Linear{OnGrid},
+        Quadratic{ExtendInner,OnCell},
+        Quadratic{Flat,OnCell},
+        Quadratic{LinearBC,OnCell}
+    )
+    for EB in (
+            ExtrapError,
+            ExtrapNaN,
+            ExtrapConstant,
+            ExtrapLinear
+        )
         eval(:(function getindex{T}(itp::Interpolation{T,1,$IT,$EB}, x::Real, d)
             d == 1 || throw(BoundsError())
             itp[x]
@@ -107,7 +120,11 @@ for IT in (Constant{OnCell},Linear{OnGrid},Quadratic{ExtendInner,OnCell},Quadrat
 end
 
 # This creates prefilter specializations for all interpolation types that need them
-for IT in (Quadratic{ExtendInner,OnCell},Quadratic{Flat,OnCell})
+for IT in (
+        Quadratic{ExtendInner,OnCell},
+        Quadratic{Flat,OnCell},
+        Quadratic{LinearBC,OnCell}
+    )
     @ngenerate N promote_type_grid(T, x...) function prefilter{T,N}(A::Array{T,N},it::IT)
         ret = similar(A)
         szs = collect(size(A))

src/extrapolation.jl

@@ -1,4 +1,3 @@
-
 abstract ExtrapolationBehavior
 type ExtrapError <: ExtrapolationBehavior end
 
@@ -25,3 +24,26 @@ function extrap_gen(::OnGrid, ::ExtrapConstant, N)
     end
 end
 extrap_gen(::OnCell, e::ExtrapConstant, N) = extrap_gen(OnGrid(), e, N)
+
+type ExtrapLinear <: ExtrapolationBehavior end
+
+function extrap_gen(::OnGrid, ::ExtrapLinear, N)
+    quote
+        @nexprs $N d->begin
+            if x_d < 1
+                fx_d = x_d - convert(typeof(x_d), 1)
+
+                k = itp[1] - itp[2]
+                return itp[1] - k * fx_d
+            end
+            if x_d > size(itp, d)
+                s_d = size(itp,d)
+                fx_d = x_d - convert(typeof(x_d), s_d)
+
+                k = itp[s_d] - itp[s_d - 1]
+                return itp[s_d] + k * fx_d
+            end
+        end
+    end
+end
+extrap_gen(::OnCell, e::ExtrapLinear, N) = extrap_gen(OnGrid(), e, N)

src/quadratic.jl

@@ -26,16 +26,6 @@ function indices(::Quadratic{ExtendInner,OnCell}, N)
     end
 end
 
-function coefficients(::Quadratic{ExtendInner,OnCell}, N)
-    quote
-        @nexprs $N d->begin
-            cm_d = (fx_d-.5)^2 / 2
-            c_d = 0.75 - fx_d^2
-            cp_d = (fx_d+.5)^2 / 2
-        end
-    end
-end
-
 function bc_gen(::Quadratic{Flat,OnCell}, N)
     quote
         # After extrapolation has been handled, 1 <= x_d <= size(itp,d)
@@ -63,8 +53,42 @@ function indices(::Quadratic{Flat,OnCell}, N)
     end
 end
 
+function bc_gen(::Quadratic{LinearBC,OnCell}, N)
+    quote
+         @nexprs $N d->begin
+            if x_d < 1
+                # extrapolate towards -∞
+                fx_d = x_d - convert(typeof(x_d), 1)
+                k = itp[2] - itp[1]
+
+                return itp[1] + k * fx_d
+            end
+            if x_d > size(itp, d)
+                # extrapolate towards ∞
+                s_d = size(itp, d)
+                fx_d = x_d - convert(typeof(x_d), s_d)
+                k = itp[s_d] - itp[s_d - 1]
+
+                return itp[s_d] + k * fx_d
+            end
+
+            ix_d = iround(x_d)
+         end
+    end
+end
 
-function coefficients(::Quadratic{Flat,OnCell}, N)
+function indices(::Quadratic{LinearBC,OnCell}, N)
+    quote
+        @nexprs $N d->begin
+            ixp_d = ix_d + 1
+            ixm_d = ix_d - 1
+
+            fx_d = x_d - convert(typeof(x_d), ix_d)
+        end
+    end
+end
+
+function coefficients(::Quadratic, N)
     quote
         @nexprs $N d->begin
             cm_d = (fx_d-.5)^2 / 2
@@ -114,3 +138,11 @@ function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnCe
     dl[end] += 1/8
     lufact!(Tridiagonal(dl, d, du))
 end
+
+function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{LinearBC,OnCell})
+    dl,d,du = unmodified_system_matrix(T,n,q)
+
+    d[1] = d[end] = 1
+    du[1] = dl[end] = 0
+    lufact!(Tridiagonal(dl, d, du))
+end