dg1d: Add Bsplines (!210)

This adds multiple things:
- A `Bspline` module for working with 2nd order Bsplines.
- Extends `DGElement` with the `kind=:modal_bspline2` option.
- Updates docstrings and refactors some tests.
- Adds a `LG` module for Legendre-Gauss quadrature.
- Adds a plotting script `plot/bspline.jl` to tinker with the `Bspline2` approximation.

Also adds a new option to the `Mesh` parameters section:
- `kind = "modal_bspline2"`: Utilize a 2nd order modal Bspline Ansatz for the DG evolution.

plot/bspline2.jl

+using dg1d
+using dg1d.Jacobi
+using GLMakie
+using LinearAlgebra
+
+
+# bs = dg1d.Bspline.Bspline2(5)
+x, _ = dg1d.LGL.rule(5)
+bs = dg1d.Bspline.Bspline2(x)
+xx = collect(range(-2.0,2.0,length=300))
+Bs = [ dg1d.Bspline.polynomial(xx[i], j, bs) for i in 1:length(xx), j in 1:bs.Npts+1 ]
+∂Bs = [ dg1d.Bspline.derivative_polynomial(xx[i], j, bs) for i in 1:length(xx), j in 1:bs.Npts+1 ]
+
+Vs = [ dg1d.Bspline.polynomial(bs.xs[2+i], j, bs) for i in 1:bs.Npts, j in 1:bs.Npts+1 ]
+# display(Vs)
+∂Vs = [ dg1d.Bspline.derivative_polynomial(bs.xs[2+i], j, bs) for i in 1:bs.Npts, j in 1:bs.Npts+1 ]
+# display(∂Vs)
+
+M = dg1d.Bspline.mass_matrix(bs)
+S = dg1d.Bspline.stiffness_matrix(bs)
+ST = transpose(S)
+BB = similar(M); fill!(BB, 0.0)
+o = 2
+for i in 1:bs.Npts+1, j in 1:bs.Npts+1
+  for k in 1:bs.Npts-1
+    xl, xr = bs.xs[o+k], bs.xs[o+k+1]
+    Nil, Nir = dg1d.Bspline.polynomial(xl, i, bs), dg1d.Bspline.polynomial(xr, i, bs)
+    Njl, Njr = dg1d.Bspline.polynomial(xl, j, bs), dg1d.Bspline.polynomial(xr, j, bs)
+    BB[i,j] += Nir*Njr - Nil*Njl
+  end
+end
+RBR = dg1d.purge_zeros(S+ST)
+display(M)
+display(S)
+display(RBR)
+display(BB)
+
+if any(isnan, Bs)
+  error("polynomial produces NaNs")
+end
+if any(isnan, ∂Bs)
+  error("derivative_polynomial produces NaNs")
+end
+if any(isnan, Vs)
+  error("polynomial produces NaNs")
+end
+if any(isnan, ∂Vs)
+  error("derivative_polynomial produces NaNs")
+end
+
+fig = Figure()
+ax = Axis(fig[1,1])
+vlines!(ax, bs.xs, color=(:black,0.5), linestyle=:dash)
+
+for (i,(B,V)) in enumerate(zip(eachcol(Bs),eachcol(Vs)))
+  lines!(ax, xx, B, label="$i", color=Makie.wong_colors()[i])
+  scatter!(ax, bs.xs[3:2+bs.Npts], V, color=Makie.wong_colors()[i])
+end
+PoU = vec(sum(Bs, dims=2))
+lines!(ax, xx, PoU, label=L"\Sigma_n B_n", linestyle=:dot)
+
+# for (i,(∂B,∂V)) in enumerate(zip(eachcol(∂Bs),eachcol(∂Vs)))
+#   lines!(ax, xx, ∂B, label="$i", linestyle=:dash, color=Makie.wong_colors()[i])
+#   scatter!(ax, bs.xs[3:2+bs.Npts], ∂V, color=Makie.wong_colors()[i])
+# end
+
+# FD_∂B = diff(Bs, dims=1) ./ diff(xx)
+# xx_ctr = xx[1:end-1] .+ (xx[2]-xx[1])/2
+# for (i,∂B) in enumerate(eachcol(FD_∂B))
+#   lines!(ax, xx_ctr, ∂B, label="$i", color=Makie.wong_colors()[i])
+# end
+
+axislegend(ax)
+display(fig)
+
+
+# # fn = x -> x < 0.15 ? -x+1 : 0.0
+# # fn = x -> x < 0.0 ? -x : 0.0
+Npts = 4
+xs, _ = dg1d.LGL.rule(Npts)
+xx = collect(range(extrema(xs)...,length=100))
+fn = x -> x < -0.5 ? -x+1 : 0.0
+fx = fn.(xs)
+
+knots = vcat(xs[1], (xs[1:end-1].+xs[2:end])./2, xs[end])
+avg_fx = vcat(fx[1], (fx[1:end-1].+fx[2:end])./2, fx[end])
+bs = dg1d.Bspline.Bspline2(xs)
+Vs = dg1d.Bspline.vandermonde_matrix(knots, bs)
+intrp_Bs = [ dg1d.Bspline.polynomial(xx[i], j, bs) for i in 1:length(xx), j in 1:bs.Npts+1 ]
+A_lsqr = transpose(Vs)*Vs
+b_lsqr = transpose(Vs)*avg_fx
+
+# nodal Bspline
+nodal_Bs = Vs \ avg_fx
+fn_nodal_Bs = intrp_Bs * nodal_Bs
+
+# quasi-nodal Bspline
+quasi_nodal_Bs = avg_fx
+fn_quasi_nodal_Bs = intrp_Bs * quasi_nodal_Bs
+
+# corrected nodal Bspline
+c_nodal_Bs = max.(nodal_Bs, 0.0)
+fn_c_nodal_Bs = intrp_Bs * c_nodal_Bs
+
+# least-square Bspline
+lsqr_Bs = A_lsqr \ b_lsqr
+fn_lsqr_Bs = intrp_Bs * lsqr_Bs
+
+# corrected least-sqaure Bspline
+c_lsqr_Bs = max.(lsqr_Bs, 0.0)
+fn_c_lsqr_Bs = intrp_Bs * c_lsqr_Bs
+
+# nodal Lagrange
+intrp_La = dg1d.lagrange_interp_matrix(xx, xs)
+fn_La = intrp_La * fx
+
+# modal Legendre
+Ln = dg1d.Legendre.projection(fn, Npts)
+fn_Ln = dg1d.Legendre.interpolate.(xx, Ref(Ln))
+
+fig = Figure()
+ax = Axis(fig[1,1])
+
+vlines!(ax, xs, color=(:black,0.4), linestyle=:dash)
+# lines!(ax, xx, fn_Ln, color=:blue, label="Legendre")
+lines!(ax, xx, fn_La, color=:orange, label="Lagrange")
+scatterlines!(ax, xs, fx, color=:orange, markersize=22, linestyle=:dash)
+
+# lines!(ax, xx, fn_nodal_Bs, color=:blue, label="nodal Bspline")
+# scatterlines!(ax, knots, nodal_Bs, color=:blue, markersize=18, linestyle=:dash)
+
+lines!(ax, xx, fn_quasi_nodal_Bs, color=:green, label="quasi-nodal Bspline")
+scatterlines!(ax, knots, quasi_nodal_Bs, color=:green, markersize=14, linestyle=:dash)
+
+# lines!(ax, xx, fn_c_nodal_Bs, color=:red, label="corrected nodal Bspline")
+# scatterlines!(ax, knots, c_nodal_Bs, color=:red, markersize=10, linestyle=:dash)
+
+lines!(ax, xx, fn_lsqr_Bs, color=:violet, label="least-square Bspline")
+scatterlines!(ax, knots, lsqr_Bs, color=:violet, markersize=14, linestyle=:dash)
+
+lines!(ax, xx, fn_c_lsqr_Bs, color=:purple, label="corrected least-square Bspline")
+scatterlines!(ax, knots, c_lsqr_Bs, color=:purple, markersize=14, linestyle=:dash)
+
+axislegend(ax)
+display(fig)

src/bernstein.jl

@@ -110,7 +110,7 @@ end
   interpolate(x::Real, fn::AbstractVector)
 
 Interpolate a Bernstein polynomial with quasi-nodal coefficients `fn` onto `x`.
-`x` is assumed to lie within `[-1,1`].
+`x` is assumed to lie within `[-1,1]`.
 The `fn` are assumed to be sampled on the nodes `xs = [ -1.0 + 2*i/N for i in 0:N ]`
 where `N = length(fn)-1`.
 """

src/bspline.jl

+module Bspline
+
+
+using dg1d
+using Jacobi
+
+
+struct Bspline2
+  Npts::Int64 # number of unique knots
+  xs::Vector{Float64} # the knots, start and end point have multiplicity three
+end
+
+
+function Bspline2(z::Vector{Float64})
+  Npts = length(z)
+  if Npts < 2
+    error("z must contain at least two points, but you provided $Npts")
+  end
+  if any(d -> d < 0.0, diff(z))
+    error("z must be a strictly increasing point set, but you provided $z")
+  end
+  xs = vcat([z[1],z[1]], z, [z[end],z[end]])
+  return Bspline2(Npts, xs)
+end
+
+
+function polynomial(x::Float64, i::Int64, bs::Bspline2)
+  @unpack Npts, xs = bs
+  1 ≤ i ≤ Npts+1 || return 0.0
+  @inbounds begin
+    xi, xi1, xi2, xi3 = xs[i], xs[i+1], xs[i+2], xs[i+3]
+  end
+  if x < xi || x > xi3
+    return 0.0
+  end
+  if i == 1 && x == xi2
+    # compute derivative on lhs bdry from the right
+    @goto sub3
+  end
+  if i == 2 && x == xi1
+    # compute derivative on lhs bdry from the right
+    @goto sub2
+  end
+  if x ≤ xi1
+    return (x-xi)^2/((xi2-xi)*(xi1-xi))
+  end
+  if x ≤ xi2
+    @label sub2
+    return (x-xi)*(xi2-x)/((xi2-xi)*(xi2-xi1)) + (xi3-x)*(x-xi1)/((xi3-xi1)*(xi2-xi1))
+  end
+  if x ≤ xi3
+    @label sub3
+    return (xi3-x)^2/((xi3-xi1)*(xi3-xi2))
+  end
+  # unreachable
+  return NaN
+end
+
+
+function derivative_polynomial(x::Float64, i::Int64, bs::Bspline2)
+  @unpack Npts, xs = bs
+  1 ≤ i ≤ Npts+1 || return 0.0
+  @inbounds begin
+    xi, xi1, xi2, xi3 = xs[i], xs[i+1], xs[i+2], xs[i+3]
+  end
+  if x < xi || x > xi3
+    return 0.0
+  end
+  if i == 1 && x == xi2
+    # compute derivative on lhs bdry from the right
+    @goto sub3
+  end
+  if i == 2 && x == xi1
+    # compute derivative on lhs bdry from the right
+    @goto sub2
+  end
+  if x ≤ xi1
+    return 2*(x-xi)/((xi2-xi)*(xi1-xi))
+  end
+  if x ≤ xi2
+    @label sub2
+    return 2*(xi2-x)/((xi2-xi)*(xi2-xi1)) - 2*(x-xi1)/((xi3-xi1)*(xi2-xi1))
+  end
+  if x ≤ xi3
+    @label sub3
+    return -2*(xi3-x)/((xi3-xi1)*(xi3-xi2))
+  end
+  # unreachable
+  return NaN
+end
+
+
+"""
+  interpolate(x::Real, fn::AbstractVector, bs::Bspline2)
+
+Interpolate a Bspline polynomial with quasi-nodal coefficients `fn` onto `x`.
+`x` is assumed to lie within `[-1,1]`.
+The `fn` are assumed to be sampled on the nodes `xs = [ -1.0 + 2*i/N for i in 0:N ]`
+where `N = length(fn)-1`.
+"""
+function interpolate(x::Real, coeffs::AbstractVector, bs::Bspline2)
+  @unpack Npts, xs = bs
+  if length(coeffs) != Npts+1
+    error("expected Npts+1 = $Npts coefficients, received $(length(coeffs))")
+  end
+  return unsafe_interpolate(x, coeffs, bs)
+end
+function unsafe_interpolate(x::Real, coeffs::AbstractVector, bs::Bspline2)
+  @unpack Npts, xs = bs
+  @inbounds xl, xr = xs[1], xs[end]
+  xl ≤ x ≤ xr || return 0.0
+  o = 2 # offset duplicated bdry knots
+  v_xs = view(xs, o .+ (1:Npts))
+  ir = searchsortedfirst(v_xs, x)
+  il = ir-1
+  il < 1 && return coeffs[1] # x == xl
+  is_supp = max(il-1,1):min(ir+1,Npts+1)
+  B = 0.0
+  @inbounds for i in is_supp
+    B += coeffs[i] * polynomial(x, i, bs)
+  end
+  return B
+end
+
+
+"""
+  mass_matrix(bs::Bspline2)
+
+Compute `M_ij = (Bi, Bj)` exactly where `Bi, Bj` are `i`th and `j`th Bspline of order 2.
+"""
+function mass_matrix(bs::Bspline2)
+  @unpack Npts, xs = bs
+  M = zeros(Float64, Npts+1, Npts+1)
+  zs, ws = dg1d.LG.rule(2)
+  o = 2 # offset bdry knots
+  for k in 1:Npts-1
+    il, ir = o+k, o+k+1
+    xl, xr = xs[il], xs[ir]
+    for i in 1:Npts+1, j in 1:Npts+1
+      mij = dg1d.LG.integrate(zs, ws) do z
+        x = ((xr-xl)*z + xr+xl)/2
+        dxdz = (xr-xl)/2
+        Bi = polynomial(x, i, bs)
+        Bj = polynomial(x, j, bs)
+        Bi*Bj*dxdz
+      end
+      M[i,j] += mij
+    end
+  end
+  return M
+end
+
+
+"""
+  stiffness_matrix(bs::Bspline2)
+
+Compute `S_ij = (Bi, dBj/dx)` exactly where `Bi, Bj` are `i`th and `j`th Bspline of order 2.
+"""
+function stiffness_matrix(bs::Bspline2)
+  @unpack Npts, xs = bs
+  S = zeros(Float64, Npts+1, Npts+1)
+  zs, ws = dg1d.LG.rule(3)
+  o = 2 # offset bdry knots
+  for k in 1:Npts-1
+    il, ir = o+k, o+k+1
+    xl, xr = xs[il], xs[ir]
+    for i in 1:Npts+1, j in 1:Npts+1
+      sij = dg1d.LG.integrate(zs, ws) do z
+        x = ((xr-xl)*z + xr+xl)/2
+        dxdz = (xr-xl)/2
+        Bi = polynomial(x, i, bs)
+        ∂Bj = derivative_polynomial(x, j, bs)
+        Bi*∂Bj*dxdz
+      end
+      S[i,j] += sij
+    end
+  end
+  return S
+end
+
+
+function vandermonde_matrix(z::AbstractVector, bs::Bspline2)
+  return [ Bspline.polynomial(z[i], j, bs) for i in 1:length(z), j in 1:bs.Npts+1 ]
+end
+
+
+# analytical integration formula from
+# Splines and PDEs: From Approximation Theory to Numerical Linear Algebra 2017, p. 13
+function integrate(f::AbstractVector, bs::Bspline2)
+  @toggled_assert length(f) == bs.Npts+1
+  sum = 0.0
+  @inbounds for i in bs.Npts+1
+    fi, xl, xr = f[i], bs.xs[i], bs.xs[i+3]
+    sum += fi * (xr-xl)/3
+  end
+  return sum
+end
+function integrate(f::Function, idxs::AbstractVector{<:Integer}, bs::Bspline2)
+  @toggled_assert length(idxs) == bs.Npts+1
+  sum = 0.0
+  @inbounds for i in bs.Npts+1
+    fi, xl, xr = f[idxs[i]], bs.xs[i], bs.xs[i+3]
+    sum += fi * (xr-xl)/3
+  end
+  return sum
+end
+
+
+end

src/dg1d.jl

@@ -51,7 +51,9 @@ include("output.jl")
 include("lagrange.jl")
 include("glgl.jl")
 include("lgl.jl")
+include("lg.jl")
 include("bernstein.jl")
+include("bspline.jl")
 include("legendre.jl")
 
 #######################################################################

src/dgelement.jl

+# TODO Make DGElement an abstract type and add the following subtypes
+# - NodalDGElement
+# - ModalBernsteinDGElement
+# - ModalBspline2DGElement
 struct DGElement
   N::Int64                # polynomial order
   Npts::Int64             # number of quadrature points = N + 1
@@ -17,7 +21,7 @@ struct DGElement
   S::Matrix{Float64}      # stiffness matrix
   M::Matrix{Float64}      # mass matrix
   invM::Matrix{Float64}   # inverse mass matrix
-  MDM::Matrix{Float64}    # = (M^-1 D^T M) = M^-1 S^T
+  MDM::Matrix{Float64}    # = (M^-1 D^T M) = M^-1 S^T (TODO: this identity only applies to nodal DG!!!)
   Ml_lhs::Vector{Float64} # = M^-1 l(-1)
   Ml_rhs::Vector{Float64} # = M^-1 l(+1)
 
@@ -25,13 +29,23 @@ struct DGElement
   V::Matrix{Float64}      # = P_(j-1)(z_i), Vandermonde matrix
   invV::Matrix{Float64}   # inv(V), inverse Vandermonde matrix
   F::Matrix{Float64}      # spectral filtering matrix
+  bs::Bspline.Bspline2    # only set for :modal_bspline
 
+  # TODO Replace N with Npts and remove N where ever possible.
   function DGElement(N::Integer, quadr::Symbol, kind::Symbol;
       # filter parameters
       ηc::Float64=0.5, sh::Float64=16.0, α::Float64=36.0)
 
+    if quadr ∉ (:lgl, :lgl_lumped, :glgl)
+      error("Invalid quadrature rule specified for kind $kind: $quadr. Use one of lgl_lumped, lgl, glgl.")
+    end
+    if kind ∉ (:nodal_lagrange, :modal_bernstein, :modal_bspline2)
+      error("Invalid DG formulation specified: $kind. Use one of nodal_lagrange, modal_bernstein, modal_bspline2.")
+    end
+
     Npts = N + 1
 
+    local bs
     if kind === :nodal_lagrange
       if quadr === :lgl
         quadr_z, quadr_w = LGL.rule(Npts)
@@ -51,14 +65,12 @@ struct DGElement
         M    = GLGL.mass_matrix(Npts)
         S    = GLGL.stiffness_matrix(Npts)
         invM = inv(M)
-      else
-        error("Invalid quadrature rule specified for kind $kind: $quadr. Use one of lgl_lumped, lgl, glgl.")
       end
       z = quadr_z
       D = first_derivative_matrix(z)
       V = vandermonde_matrix_legendre(z)
-      l_lhs = zeros(Float64, N+1)
-      l_rhs = zeros(Float64, N+1)
+      l_lhs = zeros(Float64, Npts)
+      l_rhs = zeros(Float64, Npts)
       l_lhs[1]   = 1.0
       l_rhs[end] = 1.0
     elseif kind === :modal_bernstein
@@ -66,25 +78,59 @@ struct DGElement
         quadr_z, quadr_w = LGL.rule(Npts)
       elseif quadr === :glgl
         quadr_z, quadr_w = GLGL.rule(Npts)
-      else
-        error("Invalid quadrature rule specified for kind $kind: $quadr. Use one of lgl_lumped, lgl, glgl.")
       end
       purge_zeros!(quadr_z)
       z = collect(range(-1.0,1.0,Npts))
       purge_zeros!(z)
-      D = [ Bernstein.derivative_bernstein_polynomial.(z[i], j-1, N) for i in 1:Npts, j in 1:Npts ]
+      D = [ Bernstein.derivative_bernstein_polynomial(z[i], j-1, N) for i in 1:Npts, j in 1:Npts ]
       M = Bernstein.mass_matrix(Npts)
       S = Bernstein.stiffness_matrix(Npts)
       invM = inv(M)
       V = Bernstein.vandermonde_matrix(z)
-      l_lhs      = zeros(Float64, N+1)
-      l_rhs      = zeros(Float64, N+1)
+      l_lhs      = zeros(Float64, Npts)
+      l_rhs      = zeros(Float64, Npts)
+      l_lhs[1]   = 1.0
+      l_rhs[end] = 1.0
+    elseif kind === :modal_bspline2
+      if N < 2
+        error("Requested N+1=$(N+1) DOFs for DG formulation of kind $kind, but require at least 3.")
+      end
+      # User asked for Nth order polynomial sampled on Npts=N+1 points,
+      # but we can only provide N 2nd order Bsplines, which would give a total of Npts+1 DOFs.
+      # This is more than requested, so we reduce it by one.
+      # bs_z is now the grid across which the Bsplines are C^1 continuous.
+      if quadr === :lgl || quadr === :lgl_lumped
+        bs_z, _ = LGL.rule(N)
+      elseif quadr === :glgl
+        bs_z, _ = GLGL.rule(N)
+      end
+      purge_zeros!(bs_z)
+      bs = Bspline.Bspline2(bs_z)
+      # The dual grid of bs_z is the one on which we define the Npts quasi-nodal DOFs.
+      # These are the nodes on which we evaluate residuals etc.
+      z = vcat(bs_z[1], (bs_z[1:end-1].+bs_z[2:end])./2, bs_z[end])
+      purge_zeros!(z)
+      # Ignore given quadrature rule since we have to evaluate integrals on the subintervals,
+      # for which a 4th order accurate LG rule should suffice.
+      # We use a Gauss rule to avoid evaluating the Bspline on subcell interfaces and bdrys.
+      quadr_z, quadr_w = LG.rule(2)
+      D = [ Bspline.derivative_polynomial(z[i], j, bs) for i in 1:Npts, j in 1:Npts ]
+      M = Bspline.mass_matrix(bs)
+      S = Bspline.stiffness_matrix(bs)
+      invM = inv(M)
+      V = Bspline.vandermonde_matrix(z, bs)
+      l_lhs      = zeros(Float64, Npts)
+      l_rhs      = zeros(Float64, Npts)
       l_lhs[1]   = 1.0
       l_rhs[end] = 1.0
-    else
-      error("Invalid DG formulation specified: $kind. Use on of nodal_lagrange, modal_bernstein.")
     end
 
+    if kind !== :modal_bspline2
+      # just set up something, its unused for the other kinds
+      bs = Bspline.Bspline2(z)
+    end
+
+    # TODO Rename MDM to MST or so, because the identity MDM = M^{-1} S^{T} only applies for nodal DG
     MDM    = invM * transpose(S)
     Ml_lhs = invM * l_lhs
     Ml_rhs = invM * l_rhs
@@ -98,13 +144,15 @@ struct DGElement
     purge_zeros!(invV)
     purge_zeros!(Ml_lhs)
     purge_zeros!(Ml_rhs)
+
+    # TODO Filter operators should likely go into separate struct
     Λ = diagm([ spectral_filter((i-1)/N, (ηc,sh,α)) for i = 1:Npts ])
     F = V * Λ * invV
 
     return new(N, Npts, kind,
                quadr, quadr_z, quadr_w,
                z, D, S, M, invM, MDM, Ml_lhs, Ml_rhs,
-               V, invV, F)
+               V, invV, F, bs)
   end
 end
 
@@ -131,19 +179,27 @@ end
 
 
 @inline function integrate(u::AbstractVector, el::DGElement)
-  if el.quadr === :lgl || el.quadr === :lgl_lumped
-    return LGL.integrate(u, el.quadr_z, el.quadr_w)
-  else # el.quadr === :glgl
-    return GLGL.integrate(u, el.quadr_z, el.quadr_w)
+  if el.kind === :modal_bspline2
+    return Bspline.integrate(u, el.bs)
+  else
+    if el.quadr === :lgl || el.quadr === :lgl_lumped
+      return LGL.integrate(u, el.quadr_z, el.quadr_w)
+    else # el.quadr === :glgl
+      return GLGL.integrate(u, el.quadr_z, el.quadr_w)
+    end
   end
 end
 
 
 @inline function integrate(fn::Function, idxs::AbstractVector{<:Integer}, el::DGElement)
-  if el.quadr === :lgl || el.quadr === :lgl_lumped
-    return LGL.integrate(fn, idxs, el.quadr_z, el.quadr_w)
-  else # el.quadr === :glgl
-    return GLGL.integrate(fn, idxs, el.quadr_z, el.quadr_w)
+  if el.kind === :modal_bspline2
+    return Bspline.integrate(fn, idxs, el.bs)
+  else
+    if el.quadr === :lgl || el.quadr === :lgl_lumped
+      return LGL.integrate(fn, idxs, el.quadr_z, el.quadr_w)
+    else # el.quadr === :glgl
+      return GLGL.integrate(fn, idxs, el.quadr_z, el.quadr_w)
+    end
   end
 end
 

src/legendre.jl

@@ -74,7 +74,7 @@ end
 
 
 """
-  legendre_projection(f::Function, Npts::Integer)
+  projection(f::Function, Npts::Integer)
 
 Compute the modal Legendre coefficients `f` using a LGL quadrature rule.
 If a smooth `f` is exactly representable as an `N`th polynomial,

src/lg.jl

+module LG
+
+
+using Jacobi, LinearAlgebra, ToggleableAsserts
+
+
+"""
+  rule(Npts::Integer)
+
+Computes quadrature nodes and weights for the Legendre-Gauss quadrature rule.
+"""
+function rule(Npts::Integer)
+  α = 0.0; β = 0.0
+  x = zgj(Npts, α, β)
+  w = wgj(x, α, β)
+  return x, w
+end
+
+
+"""
+  integrate(f, x, w)
+
+Integrate `f` sampled on LG quadrature nodes `x` with quadrature weights `w`.
+Use `rule` to obtain `x, w`.
+"""
+function integrate(f::AbstractVector, x::AbstractVector, w::AbstractVector)
+  @toggled_assert length(w) == length(f) == length(x)
+  return dot(w, f)
+end
+
+
+"""
+  integrate(f::Function, idxs::AbstractVector{<:Integer}, x, w)
+
+Integrate `f` sampled on LG quadrature nodes `x` with quadrature weights `w`.
+`f(idxs[i])` shall return the value of the integrand evaluated at `x[i]`.
+Use `rule` to obtain `x, w`.
+"""
+function integrate(f::Function, idxs::AbstractVector{<:Integer}, x::AbstractVector, w::AbstractVector)
+  @toggled_assert length(w) == length(idxs)
+  sum = 0.0
+  for i in 1:length(w)
+    sum += w[i] * f(idxs[i])
+  end
+  return sum
+end
+
+
+"""
+  integrate(f::Function, x, w)
+
+Integrate `f` sampled on LG quadrature nodes `x` with quadrature weights `w`.
+`f(x[i])` shall return the value of the integrand evaluated at `x[i]`.
+Use `rule` to obtain `x, w`.
+
+This is similar to `integrate(f, idxs, x, w)`, but is better suited if the integrand is
+given in closed form.
+"""
+function integrate(f::Function, x::AbstractVector, w::AbstractVector)
+  @toggled_assert length(w) == length(x)
+  sum = 0.0
+  for i in 1:length(w)
+    sum += w[i] * f(x[i])
+  end
+  return sum
+end
+
+
+"""
+  inner_product(f1, f2, x, w)
+
+Computed the inner product between `f1` and `f2` sampled on LG quadrature nodes `x`
+with quadrature weights `w`.
+Use `rule` to obtain `x, w`.
+"""
+function inner_product(f1::AbstractVector, f2::AbstractVector, x::AbstractVector, w::AbstractVector)
+  @toggled_assert length(w) == length(f1) == length(f2) == length(x)
+  sum = 0.0
+  for i in 1:length(w)
+    sum += w[i] * f1[i] * f2[i]
+  end
+  return sum
+end
+
+
+end

src/lgl.jl

@@ -8,7 +8,7 @@ import dg1d: lagrange_poly, first_derivative_matrix, purge_zeros!
 """
   rule(Npts::Integer)
 
-Computes quadrature nodes and weights for the Legendre-Gauss-Lobatto-Legendre quadrature rule.
+Computes quadrature nodes and weights for the Legendre-Gauss-Lobatto quadrature rule.
 """
 function rule(Npts::Integer)
   α = 0.0; β = 0.0

src/parameters.jl

@@ -489,17 +489,17 @@ end
   # TODO Rename basis to quadrature
   """
   Nodal basis type used for the polynomial approximation. Available options
-  - `lgl_lumped` ... mass-lumped Legendre-Gauss-Lobatto grid
-  - `lgl`        ... Legendre-Gauss-Lobatto grid
-  - `glgl`       ... Generalized Legendre-Gauss-Lobatto grid
+  - `lgl_lumped`: mass-lumped Legendre-Gauss-Lobatto grid
+  - `lgl`: Legendre-Gauss-Lobatto grid
+  - `glgl`: Generalized Legendre-Gauss-Lobatto grid
   """
   basis = "lgl_lumped"
   @check basis in [ "lgl", "glgl", "lgl_lumped" ]
 
   """
   Cellwise approximation scheme of solution. Available options
-  - `DG` ... Discontinuous Galerkin
-  - `FV` ... Finite Volume (central scheme)
+  - `DG`: Discontinuous Galerkin
+  - `FV`: Finite Volume (central scheme)
   """
   scheme = "DG"
   @check scheme in [ "DG", "FV" ]
@@ -507,14 +507,30 @@ end
   # TODO Rename this to basis
   """
   DG scheme used for the evolution. Available options
-  - `nodal_lagrange` ... Nodal approximation using Lagrange polynomials collocated with
+  - `nodal_lagrange`: Nodal approximation using Lagrange polynomials collocated with
     the quadrature nodes provided by `basis`. This is the standard method discussed
     in Hesthaven, Warburton, 2007.
-  - `modal_bernstein` ... A modal approximation using Bernstein polynomials.
+    See [1].
+  - `modal_bernstein`: A modal approximation using Bernstein polynomials.
     The coefficients correspond to the nodal values sampled on an equidistant grid,
     but the Bernstein polynomials are not interpolating.
-    Hence, this Ansatz is also refered to as quasi-nodal.
-    See Beisiegel, Behrens, Environ Earth Sci (2015) 74:7275–7284.
+    We refer to these coefficients as being quasi-nodal.
+    See [2,3,4].
+  - `modal_bspline2`: A modal approximation using 2nd order Bspline polynomials.
+    The coefficients correspond to the nodal values sampled on a grid that is dual (including the bdrys)
+    to the grid of a `nodal_lagrange` Ansatz of order `n-1`.
+    I.e. a `modal_bspline2` element has `n+1` DOF.
+    We adapt a quasi-nodal interpretation for the coefficients.
+    The bdry values correspond to the same bdry values as in a nodal Ansatz,
+    the cell interior values correspond to the cell-averages (determined by linear interpolation)
+    of a nodal Ansatz.
+    See [3,4].
+    Note: Spectral filtering with the `filter_*` parameters is not implemented correctly.
+
+  [1] Hesthaven, Warburton, Nodal Discontinuous Galerkin Methods (2008).
+  [2] Beisiegel, Behrens, Environ Earth Sci (2015) 74:7275–7284.
+  [3] Piegl, Tiller, The NURBS book (1997).
+  [4] Kunoth et al., Splines and PDEs: From Approximation Theory to Numerical Linear Algebra (2018).
   """
   dg_kind = "nodal_lagrange"
   @check scheme in [ "nodal_lagrange", "modal_bernstein" ]
@@ -556,13 +572,13 @@ end
 
   """
   Time stepping method. Available options
-  - `midpt`         ... 2nd order explicit Runge Kutta method -- midpoint formula
-  - `rk4`           ... 4th order explicit Runge Kutta method
-  - `rk4_twothirds` ... 4th order explicit Runge Kutta method -- two thirds version
-  - `rk4_ralston`   ... 4th order explicit Runge Kutta method -- Ralston's method, designed for minimal truncation error
-  - `ssprk3`        ... 3rd order strong stability preserving Runge Kutta method
-  - `lserk4`        ... 4th order low storage explicit Runge Kutta method
-  - `rkf10`         ... 10th order Runge Kutta Feagin method
+  - `midpt`: 2nd order explicit Runge Kutta method -- midpoint formula
+  - `rk4`: 4th order explicit Runge Kutta method
+  - `rk4_twothirds`: 4th order explicit Runge Kutta method -- two thirds version
+  - `rk4_ralston`: 4th order explicit Runge Kutta method -- Ralston's method, designed for minimal truncation error
+  - `ssprk3`: 3rd order strong stability preserving Runge Kutta method
+  - `lserk4`: 4th order low storage explicit Runge Kutta method
+  - `rkf10`: 10th order Runge Kutta Feagin method
   """
   method = "lserk4"
   @check method in [ "midpt", "lserk4", "rk4", "rk4_twothirds", "rk_ralston",
@@ -655,25 +671,25 @@ end
 
   """
   Data format for 0d output. Available options:
-  - h5 ... raw hdf5 output, all fields in one file
+  - `h5`: raw hdf5 output, all fields in one file
   """
   format0d = "h5"
   @check format0d in [ "h5" ]
 
   """
   Data format for 1d output. Available options:
-  - h5 ... raw hdf5 output, all fields in one file
+  - `h5`: raw hdf5 output, all fields in one file
   """
   format1d = "h5"
   @check format1d in [ "h5" ]
 
   """
   Data format for 2d output. Available options:
-  - h5 ... raw hdf5 output, all fields in one file
-  - xdmf+bin ... xdmf data description and binary data dump
-  - xdmf+h5 ... xdmf data description and hdf5 data dump
-  - xdmf ... same as xdmf+h5
-  - vtk ... vtk xml format utilizing unstructured grids (but without subcell resolution)
+  - `h5`: raw hdf5 output, all fields in one file
+  - `xdmf+bin`: xdmf data description and binary data dump
+  - `xdmf+h5`: xdmf data description and hdf5 data dump
+  - `xdmf`: same as xdmf+h5
+  - `vtk`: vtk xml format utilizing unstructured grids (but without subcell resolution)
   """
   format2d = "vtk"
   @check format2d in [ "h5", "xdmf", "xdmf+bin", "xdmf+h5", "vtk" ]

test/UnitTests/src/test_dgelement.jl

@@ -129,27 +129,6 @@ end
 end
 
 
-function lgl_linear_dependency_matrix_exact(Npts::Integer)
-  Npts2 = Npts+2
-  x2, w2 = dg1d.LGL.rule(Npts2)
-  L = zeros(N,N)
-  for i in 1:N
-    Pi = map(x->dg1d.legendre(x, i-1), x2)
-    for j in i:N
-      Pj = map(x->dg1d.legendre(x, j-1), x2)
-      lij = dg1d.LGL.integrate(Pi, Pj, x2, w2)
-      L[i,j] = lij
-      if i != j
-        L[j,i] = L[i,j]
-      end
-    end
-  end
-  dg1d.purge_zeros!(L)
-
-  return L
-end
-
-
 @testset "Test orthogonality relation" begin
 
   Npts = 6
@@ -194,10 +173,14 @@ end
   int = dg1d.GLGL.integrate(p, x, w)
   @test int ≈ 140/3
 
+  x, w = dg1d.LG.rule(Npts)
+  int = dg1d.LG.integrate(p, x, w)
+  @test int ≈ 140/3
+
 end
 
 
-@testset "verify some DG operator identities" begin
+@testset "verify summation-by-parts property of DG operators" begin
 
   N = 10
 
@@ -218,12 +201,32 @@ end
   TST = T*S_legendre*transpose(T) |> dg1d.purge_zeros
   @test all(el.S.≈TST)
 
+  el = dg1d.DGElement(N, :lgl, :modal_bspline2)
+  bs = el.bs
+  M = el.M
+  S = el.S
+  ST = transpose(S)
+  # verify the sum of surface terms
+  BB = dg1d.purge_zeros(S+ST)
+  exact_BB = similar(M); fill!(exact_BB, 0.0)
+  o = 2
+  for i in 1:bs.Npts+1, j in 1:bs.Npts+1
+    for k in 1:bs.Npts-1
+      xl, xr = bs.xs[o+k], bs.xs[o+k+1]
+      Nil, Nir = dg1d.Bspline.polynomial(xl, i, bs), dg1d.Bspline.polynomial(xr, i, bs)
+      Njl, Njr = dg1d.Bspline.polynomial(xl, j, bs), dg1d.Bspline.polynomial(xr, j, bs)
+      exact_BB[i,j] += Nir*Njr - Nil*Njl
+    end
+  end
+  dg1d.purge_zeros!(exact_BB)
+  @test all(exact_BB .≈ BB)
+
 end
 
 
-@testset "Legendre expansions and Bernstein polynomials" begin
+@testset "Bernstein consistency tests" begin
 
-  # approximate the Bernstein basis through Legendre polynomials
+  ## approximate the Bernstein basis through Legendre polynomials
   N = 10
   Npts = N+1
   Cs = zeros(Float64, Npts, Npts)
@@ -234,13 +237,36 @@ end
     end
   end
 
-  xx = collect(range(-1.0,1.0,25))
+  xx = collect(range(-1.0,1.0,length=25))
   for xi in xx
     Bs = [ dg1d.Bernstein.bernstein_polynomial(xi, i-1, N) for i in 1:Npts ] |> dg1d.purge_zeros
     approx_Bs = [ dg1d.Legendre.interpolate(xi, Cs[i,:]) for i in 1:Npts ] |> dg1d.purge_zeros
     @test all(isapprox.(Bs, approx_Bs; atol=1e-10, rtol=1e-8))
   end
 
+  ## partition of unity
+  Bs = [ dg1d.Bernstein.bernstein_polynomial(xx[i], j-1, N) for i in 1:length(xx), j in 1:Npts ]
+  POU = sum(Bs, dims=2)
+  @test all(POU .≈ ones(Float64, length(POU)))
+
+end
+
+
+@testset "Bspline consistency tests" begin
+
+  z, _ = dg1d.LGL.rule(5)
+  bs = dg1d.Bspline.Bspline2(z)
+  xx = collect(range(-1.0,1.0,length=25))
+  Bs = [ dg1d.Bspline.polynomial(xx[i], j, bs) for i in 1:length(xx), j in 1:bs.Npts+1 ]
+
+  ## partition of unity
+  POU = sum(Bs, dims=2)
+  @test all(POU .≈ ones(Float64, length(POU)))
+  # alternative test: evaluate a Bspline polynomial with all coefficients equal to one
+  coeffs = ones(Float64, bs.Npts+1)
+  POU = dg1d.Bspline.interpolate.(xx, Ref(coeffs), Ref(bs))
+  @test all(POU .≈ ones(Float64, length(POU)))
+
 end
 
 

test/UnitTests/src/test_utils.jl

 import dg1d: make_outputdir, flatten, dirname_path, query,
-             splitallext, trim_keys, match_keys, prepend_keys,
-             concretetypes, match_properties, indent
+             splitallext, trim_keys, match_keys, prepend_keys, indent
 
 
 @testset "Utilities" begin