Boundscheck for scale_by_mass_matrix! doesn't work

[JoshuaLampert]

There is a boundscheck in scale_by_mass_matrix!, but it does not trigger an error if it fails:

julia> N = 10
10

julia> D = derivative_operator(MattssonNordström2004(), 1, 4, -1, 1, N)
SBP first-derivative operator of order 4 on a grid in [-1.0, 1.0] using 10 nodes 
and coefficients of Mattsson, Nordström (2004) 
  Summation by parts operators for finite difference approximations of second 
    derivatives. 
  Journal of Computational Physics 199, pp. 503-540.

julia> x = ones(N - 1)

julia> scale_by_mass_matrix!(x, D)
9-element Vector{Float64}:
 0.0787037037037037
 0.27314814814814814
 0.19907407407407407
 0.22685185185185183
 0.2222222222222222
 0.22685185185185183
 0.19907407407407407
 0.27314814814814814
 0.0787037037037037

julia> using LinearAlgebra

julia> diag(mass_matrix(D))
10-element Vector{Float64}:
 0.0787037037037037
 0.27314814814814814
 0.19907407407407407
 0.22685185185185183
 0.2222222222222222
 0.2222222222222222
 0.22685185185185183
 0.19907407407407407
 0.27314814814814814
 0.0787037037037037

This means, in this case, one of the central scaled nodes is simply omitted in the updated vector, which can lead to erroneous behavior.