Generalize Neumann handling to higher dimensions

[ChrisRackauckas]

In SciML/DiffEqOperators.jl#349 there is a portion:

        derivars = [[dot(left_weights(j),[depvars[j][central_neighbor_idxs(CartesianIndex(2),1)[1:2][2]],depvars[j][central_neighbor_idxs(CartesianIndex(2),1)[1:2][1]]]),
                    dot(right_weights(j),[depvars[j][central_neighbor_idxs(CartesianIndex(length(space[1])-1),1)[end-1:end][1]],depvars[j][central_neighbor_idxs(CartesianIndex(length(space[1])-1),1)[end-1:end][2]]])]
                    for j in 1:length(pdesys.depvars)]

This is hardcoded to only be the inward facing derivative in a one dimensional sense because CartesianIndex(2) is going to fail against a higher dimensional point. Some starter code is:

edgeindices = [indices[e...] for e in edges]
depvarmaps = reduce(vcat,[substitute.((pdesys.depvars[i],),edgevals) .=> edgevars[i] for i in 1:length(pdesys.depvars)])
left_weights(j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, [space[j][1],space[j][2]])
right_weights(j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, [space[j][end-1],space[j][end]])
central_neighbor_idxs(i,j) = [i+CartesianIndex([ifelse(l==j,-1,0) for l in 1:length(nottime)]...),i,i+CartesianIndex([ifelse(l==j,1,0) for l in 1:length(nottime)]...)]
derivars = [[[[dot(left_weights(j),[depvars[j][central_neighbor_idxs(edgeindices[i][k],1)[1:2][2]],depvars[j][central_neighbor_idxs(edgeindices[i][k],1)[1:2][1]]]),
             dot(right_weights(j),[depvars[j][central_neighbor_idxs(CartesianIndex(length(space[1])-1),1)[end-1:end][1]],depvars[j][central_neighbor_idxs(CartesianIndex(length(space[1])-1),1)[end-1:end][2]]])]
             for j in 1:length(pdesys.depvars)] for k in 1:length(edgeindices[i])] for i in 1:length(edgeindices)]

[ChrisRackauckas]

There are two major difficulties: one is that the boundaries become whole arrays, and the other is that it needs to perturb based off of higher dimensional Cartesians. The Cartesians I show can be handled by grabbing edgeindices like above, so that part isn't hard. So all that's left is munging it into a format so the rules broadcast in the same way as the Dirichlet rules, which is surprisingly difficult but probably just a one line modification in there.

[valentinsulzer]

This is hardcoded to only be the inward facing derivative

FYI this was the thing that made the Neumann and Robin BCs fail. It should just be the derivative in the direction of the dimension

[xtalax]

Fixed by #17, Needs tests