[WIP] Correct concretization BCs and ghost derivative operators

[jagot]

Fixes #182.

So far, the different concretizations of boundary conditions have been implemented, ghost derivative operators are up next.

Please check u for PeriodicBCs, I was unsure if they should be padded with zero(T) or one(T); they were of incorrect length N, need to be N+2.

[xtalax]

@jagot

Just a heads up, The GhostDerivativeOperators are going to change quite a lot once #170 is merged, as they currently assume that the domain is one dimensional - this PR sets them up to work in multiple dimensions. The concretizations for these are already implemented there - just needing a concretization for the multiple dimensional boundary conditions to banded.

[jagot]

@xtalax Thanks for the heads-up! As you see in the latest commit to this PR, I did not do anything big, just cleaned the concretization code.

However, I did not yet figure out how to shrink the bandwidths properly. You can remove bstl from this line, which will give the output below, but will break some of the tests, so some more thinking on my part is to be done.

julia> L = 10.0
10.0

julia> N = 9
9

julia> δx = L/(N+1)
1.0

julia> Δ = CenteredDifference(2, 2, δx, N)
DerivativeOperator{Float64,1,false,Float64,StaticArrays.SArray{Tuple{3},Float64,1,3},StaticArrays.SArray{Tuple{0},StaticArrays.SArray{Tuple{4},Float64,1,4},1,0},Nothing,Nothing}(2, 2, 1.0, 9, 3, [1.0, -2.0, 1.0], 4, 0, StaticArrays.SArray{Tuple{4},Float64,1,4}[], StaticArrays.SArray{Tuple{4},Float64,1,4}[], nothing, nothing)

julia> Q = Dirichlet0BC(typeof(δx))
RobinBC{Float64,Array{Float64,1}}([-0.0, 0.0], 0.0, [-0.0, 0.0], 0.0)

julia> A = Δ*Q;

julia> first(sparse(A))
9×9 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
 -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅
  1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅
   ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅
   ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅
   ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅
   ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅
   ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0

However, Neumann BCs (et al.) do not generate optimal bandwidths, since banded–banded multiplication has no way of knowing that most of the sub-/superpdiagonals are actually zero:

julia> Q = Neumann0BC(δx)
RobinBC{Float64,Array{Float64,1}}([1.0], -0.0, [1.0], 0.0)

julia> A = Δ*Q;

julia> first(sparse(A))
9×9 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
 -1.0   1.0   0.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅
  1.0  -2.0   1.0   0.0    ⋅     ⋅     ⋅     ⋅     ⋅
  0.0   1.0  -2.0   1.0   0.0    ⋅     ⋅     ⋅     ⋅
   ⋅    0.0   1.0  -2.0   1.0   0.0    ⋅     ⋅     ⋅
   ⋅     ⋅    0.0   1.0  -2.0   1.0   0.0    ⋅     ⋅
   ⋅     ⋅     ⋅    0.0   1.0  -2.0   1.0   0.0    ⋅
   ⋅     ⋅     ⋅     ⋅    0.0   1.0  -2.0   1.0   0.0
   ⋅     ⋅     ⋅     ⋅     ⋅    0.0   1.0  -2.0   1.0
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    0.0   1.0  -1.0

No tests for this yet.