Infrastructure for padding the equation system

The extra equations from the boundary conditions will not always allow indexing as we did before, by e.g. doing cm_d = cp_d at the lower boundary for some BC:s. In order to generalize better, we instead pad the equation system with one line for each extra equation to allow more general BC specs. This had the advantage of significantly simplifying some of the code (and, even better, some of the concepts) required to implement new functionality, so even though some of the underlying infrastructure is now more complicated than before, new contributors will hopefully have a lower threshold to start adding features.
JuliaMath · Dec 26, 2014 · 0c98fef · 0c98fef
1 parent 535e54f
commit 0c98fef
Show file tree

Hide file tree

Showing 2 changed files with 70 additions and 137 deletions.
diff --git a/src/Interpolations.jl b/src/Interpolations.jl
@@ -2,7 +2,7 @@ module Interpolations
 
 using Base.Cartesian
 
-import Base: size, eltype, getindex
+import Base: size, eltype, getindex, ndims
 
 export 
     Interpolation,
@@ -50,8 +50,10 @@ end
 # However, all prefilters copy the array, so do that here as well
 prefilter{T,N,IT<:InterpolationType}(A::AbstractArray{T,N}, ::IT) = copy(A)
 
-size(itp::Interpolation, d::Integer) = size(itp.coefs, d)
-size(itp::Interpolation) = size(itp.coefs)
+size{T,N,IT<:InterpolationType}(itp::Interpolation{T,N,IT}, d::Integer) =
+    size(itp.coefs, d) - 2*padding(IT())
+size(itp::Interpolation) = tuple(collect([size(itp,i) for i in 1:ndims(itp)]))
+ndims(itp::Interpolation) = ndims(itp.coefs)
 eltype(itp::Interpolation) = eltype(itp.coefs)
 
 offsetsym(off, d) = off == -1 ? symbol(string("ixm_", d)) :
@@ -68,6 +70,14 @@ include("constant.jl")
 include("linear.jl")
 include("quadratic.jl")
 
+# If nothing else is specified, don't pad at all
+padding(::InterpolationType) = 0
+padding{T,N,IT<:InterpolationType}(::Interpolation{T,N,IT}) = padding(IT())
+function similar_with_padding(A, it::InterpolationType)
+    pad = padding(it)
+    coefs = Array(eltype(A), [size(A,i)+2pad for i in 1:ndims(A)]...)
+    coefs, pad
+end
 
 # This creates getindex methods for all supported combinations
 for IT in (
@@ -130,22 +140,28 @@ for IT in (
         Quadratic{LinearBC,OnCell}
     )
     @ngenerate N promote_type_grid(T, x...) function prefilter{T,N}(A::Array{T,N},it::IT)
-        ret = similar(A)
+        ret, pad = similar_with_padding(A,it)
         szs = collect(size(A))
         strds = collect(strides(A))
+        strdsR = collect(strides(ret))
 
         for dim in 1:N
             n = szs[dim]
             szs[dim] = 1
 
-            M = prefiltering_system_matrix(eltype(A), n, it)
+            M, b = prefiltering_system(eltype(A), n+2pad, it)
 
             @nloops N i d->1:szs[d] begin
                 cc = @ntuple N i
                 strt = 1 + sum([(cc[i]-1)*strds[i] for i in 1:length(cc)])
+                strtR = 1 + sum([(cc[i]-1)*strdsR[i] for i in 1:length(cc)])
                 rng = range(strt, strds[dim], n)
-                ret[rng] = M \ vec(A[rng])
-            end            
+                rngR = range(strtR, strdsR[dim], size(ret,dim))
+
+                b[1+pad:end-pad] += A[rng]
+                ret[rngR] = M \ b
+                b[1+pad:end-pad], A[rng]
+            end
             szs[dim] = n
         end
         ret

diff --git a/src/quadratic.jl b/src/quadratic.jl
@@ -12,106 +12,12 @@ end
 
 function indices(q::Quadratic, N)
     quote
+        pad = padding($q)
         @nexprs $N d->begin
             ixp_d = ix_d + 1
             ixm_d = ix_d - 1
 
-            # fx_d is the distance from the middle point to the interpolation point
-            fx_d = x_d - convert(typeof(x_d), ix_d)
-        end
-    end
-end
-
-function bc_gen(::Quadratic{Flat,OnGrid}, N)
-    quote
-        # After extrapolation has been handled, 1 <= x_d <= size(itp,d)
-        # The index is simply the closest integer.
-         @nexprs $N d->(ix_d = iround(x_d))
-        # end
-    end
-end
-
-function bc_gen(::Quadratic{Flat,OnCell}, N)
-    quote
-        # After extrapolation has been handled, 0.5 <= x_d <= size(itp,d)+.5
-        @nexprs $N d->begin
-            # The index is the closest integer...
-            ix_d = iround(x_d)
-
-            #...except in the case where x_d is actually at the upper edge;
-            # since size(itp,d)+.5 is rounded toward size(itp,d)+1,
-            # it needs special treatment*.
-            if x_d == size(itp,d)+.5
-                ix_d -= 1
-            end
-        end
-    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 indices(::Quadratic{Flat,OnGrid}, 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)
-
-            if ix_d == size(itp,d)
-                ixp_d = ixm_d
-            end
-            if ix_d == 1
-                ixm_d = ixp_d
-            end
-        end
-    end
-end
-function indices(::Quadratic{Flat,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)
-
-            if ix_d == size(itp,d)
-                ixp_d = ix_d
-            end
-            if ix_d == 1
-                ixm_d = ix_d
-            end
-        end
-    end
-end
-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)
+            fx_d = x_d - convert(typeof(x_d), ix_d - pad)
         end
     end
 end
@@ -126,7 +32,7 @@ function coefficients(::Quadratic, N)
     end
 end
 
-# This assumes integral values ixm_d, ix_d, and ixp_d (typically ixm_d = ix_d-1, ixp_d = ix_d+1, except at boundaries),
+# This assumes integral values ixm_d, ix_d, and ixp_d,
 # coefficients cm_d, c_d, and cp_d, and an array itp.coefs
 function index_gen(degree::QuadraticDegree, N::Integer, offsets...)
     if length(offsets) < N
@@ -140,47 +46,58 @@ function index_gen(degree::QuadraticDegree, N::Integer, offsets...)
     end
 end
 
-function unmodified_system_matrix{T}(::Type{T}, n::Int, ::Quadratic)
+# Quadratic interpolation has 1 extra coefficient at each boundary
+# Therefore, make the coefficient array 1 larger in each direction,
+# in each dimension.
+padding(::Quadratic) = 1
+
+function inner_system_diags{T}(::Type{T}, n::Int, ::Quadratic)
     du = fill(convert(T,1/8), n-1)
     d = fill(convert(T,3/4),n)
     dl = copy(du)
+    d[1] = d[end] = du[1] = dl[end] = convert(T,0)
     (dl,d,du)
 end
 
-function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{ExtendInner})
-    dl,d,du = unmodified_system_matrix(T,n,q)
-    d[1] = d[end] = 9/8
-    du[1] = dl[end] = -1/4
-    MT = lufact!(Tridiagonal(dl, d, du))
-    U = zeros(T,n,2)
-    V = zeros(T,2,n)
-    C = zeros(T,2,2)
-
-    C[1,1] = C[2,2] = 1/8
-    U[1,1] = U[n,2] = 1.
-    V[1,3] = V[2,n-2] = 1.
-
-    Woodbury(MT, U, C, V)
-end
-
-function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnCell})
-    dl,d,du = unmodified_system_matrix(T,n,q)
-    d[1] += 1/8
-    d[end] += 1/8
-    lufact!(Tridiagonal(dl, d, du))
+function prefiltering_system{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnCell})
+    dl,d,du = inner_system_diags(T,n,q)
+    d[1] = d[end] = -1
+    du[1] = dl[end] = 1
+    lufact!(Tridiagonal(dl, d, du)), zeros(T, n)
 end
 
-function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnGrid})
-    dl,d,du = unmodified_system_matrix(T,n,q)
-    du[1] += 1/8
-    dl[end] += 1/8
-    lufact!(Tridiagonal(dl, d, du))
+function prefiltering_system{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnGrid})
+    dl,d,du = inner_system_diags(T,n,q)
+    d[1] = d[end] = convert(T, -1)
+    du[1] = dl[end] = convert(T, 0)
+
+    rowspec = zeros(T,n,2)
+    # first row     last row
+    rowspec[1,1] = rowspec[n,2] = convert(T, 1)
+    colspec = zeros(T,2,n)
+    # third col     third-to-last col
+    colspec[1,3] = colspec[2,n-2] = convert(T, 1)
+    valspec = zeros(T,2,2)
+    # val for [1,3], val for [n,n-2]
+    valspec[1,1] = valspec[2,2] = convert(T, 1)
+
+    Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(T, n)
 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))
+function prefiltering_system{T}(::Type{T}, n::Int, q::Quadratic{LinearBC})
+    dl,d,du = inner_system_diags(T,n,q)
+    d[1] = d[end] = convert(T, 1)
+    du[1] = dl[end] = convert(T, -2)
+
+    rowspec = zeros(T,n,2)
+    # first row     last row
+    rowspec[1,1] = rowspec[n,2] = convert(T, 1)
+    colspec = zeros(T,2,n)
+    # third col     third-to-last col
+    colspec[1,3] = colspec[2,n-2] = convert(T, 1)
+    valspec = zeros(T,2,2)
+    # val for [1,3], val for [n,n-2]
+    valspec[1,1] = valspec[2,2] = convert(T, 1)
+
+    Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(T, n)
 end