calculate correct boundary stencils for irregular grid

Author: fgerick

src/derivative_operators/derivative_operator.jl

@@ -90,16 +90,18 @@ function CenteredDifference{N}(derivative_order::Int,
 
     interior_x              = boundary_point_count+2:len+1-boundary_point_count
     dummy_x                 = -div(stencil_length,2) : div(stencil_length,2)-1
-    boundary_x              = -boundary_stencil_length+1:0
-
+    low_boundary_x          = [zero(T); cumsum(dx[1:boundary_stencil_length-1])]
+    high_boundary_x         = cumsum(dx[end-boundary_stencil_length+1:end])
     # Because it's a N x (N+2) operator, the last stencil on the sides are the [b,0,x,x,x,x] stencils, not the [0,x,x,x,x,x] stencils, since we're never solving for the derivative at the boundary point.
     deriv_spots             = (-div(stencil_length,2)+1) : -1
-    boundary_deriv_spots    = boundary_x[2:div(stencil_length,2)]
 
     stencil_coefs           = [convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T), generate_coordinates(i, stencil_x, dummy_x, dx))) for i in interior_x]
-    _low_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, calculate_weights(derivative_order, oneunit(T)*x0, boundary_x)) for x0 in boundary_deriv_spots]
+    _low_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T},
+                                                                  calculate_weights(derivative_order, low_boundary_x[i+1], low_boundary_x)) for i in 1:boundary_point_count]
     low_boundary_coefs      = convert(SVector{boundary_point_count},_low_boundary_coefs)
-    high_boundary_coefs     = convert(SVector{boundary_point_count},reverse(SVector{boundary_stencil_length, T}[reverse(low_boundary_coefs[i]) for i in 1:boundary_point_count]))
+    _high_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T},
+                                                                  calculate_weights(derivative_order, high_boundary_x[end-i], high_boundary_x)) for i in boundary_point_count:-1:1]
+    high_boundary_coefs      = convert(SVector{boundary_point_count},_high_boundary_coefs)
 
     coefficients            = coeff_func isa Nothing ? nothing : Vector{T}(undef,len)