diff --git a/plot/bspline2.jl b/plot/bspline.jl
similarity index 56%
rename from plot/bspline2.jl
rename to plot/bspline.jl
index 75baad63fb321461493aeee1aa43e5435f816269..1ca886e6415577f2e641644630cce7d095ccdafd 100644
--- a/plot/bspline2.jl
+++ b/plot/bspline.jl
@@ -4,16 +4,20 @@ using GLMakie
 using LinearAlgebra
 
 
+let
+if true
+
 # bs = dg1d.Bspline.Bspline2(5)
 x, _ = dg1d.LGL.rule(5)
-bs = dg1d.Bspline.Bspline2(x)
+bs = dg1d.Bspline.Bspline1(x)
+Npts = bs.Npts
 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 ]
+Bs = [ dg1d.Bspline.polynomial(xx[i], j, bs) for i in 1:length(xx), j in 1:Npts ]
+∂Bs = [ dg1d.Bspline.derivative_polynomial(xx[i], j, bs) for i in 1:length(xx), j in 1:Npts ]
 
-Vs = [ dg1d.Bspline.polynomial(bs.xs[2+i], j, bs) for i in 1:bs.Npts, 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:Npts ]
 # 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 ]
+∂Vs = [ dg1d.Bspline.derivative_polynomial(bs.xs[2+i], j, bs) for i in 1:bs.Npts, j in 1:Npts ]
 # display(∂Vs)
 
 M = dg1d.Bspline.mass_matrix(bs)
@@ -21,7 +25,7 @@ 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 i in 1:Npts, j in 1:Npts
   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)
@@ -73,42 +77,77 @@ lines!(ax, xx, PoU, label=L"\Sigma_n B_n", linestyle=:dot)
 axislegend(ax)
 display(fig)
 
+end # if
+end
+
 
+let
+if false
 # # fn = x -> x < 0.15 ? -x+1 : 0.0
 # # fn = x -> x < 0.0 ? -x : 0.0
-Npts = 4
+Npts = 5
 xs, _ = dg1d.LGL.rule(Npts)
 xx = collect(range(extrema(xs)...,length=100))
 fn = x -> x < -0.5 ? -x+1 : 0.0
 fx = fn.(xs)
 
+T = zeros(Float64, Npts+1, Npts)
+T[1,1] = T[end,end] = 1.0
+for i in 2:Npts
+  T[i,i-1:i] .= 0.5
+end
+pinvT = pinv(T)
+# display(T)
+# display(pinvT)
+
 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])
+# avg_fx = vcat(fx[1], (fx[1:end-1].+fx[2:end])./2, fx[end])
+avg_fx = T * fx
 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 ]
+intrpmat_Bs = dg1d.Bspline.vandermonde_matrix(xx, bs)
 A_lsqr = transpose(Vs)*Vs
 b_lsqr = transpose(Vs)*avg_fx
+intrpmat_Bs_to_collocs = dg1d.Bspline.vandermonde_matrix(xs, bs)
+pinv_intrpmat_Bs_to_collocs = pinv(intrpmat_Bs_to_collocs)
+# display(intrpmat_Bs_to_collocs)
+# display(pinv_intrpmat_Bs_to_collocs)
+# display(pinv_intrpmat_Bs_to_collocs*intrpmat_Bs_to_collocs)
+# display(intrpmat_Bs_to_collocs*pinv_intrpmat_Bs_to_collocs)
 
 # nodal Bspline
 nodal_Bs = Vs \ avg_fx
-fn_nodal_Bs = intrp_Bs * nodal_Bs
+fn_nodal_Bs = intrpmat_Bs * nodal_Bs
 
 # quasi-nodal Bspline
 quasi_nodal_Bs = avg_fx
-fn_quasi_nodal_Bs = intrp_Bs * quasi_nodal_Bs
+fn_quasi_nodal_Bs = intrpmat_Bs * quasi_nodal_Bs
+# quasi-nodal Bspline data
+quasi_nodal_Bs_samples = Vs * quasi_nodal_Bs
+# corrected interpolateD
+c_quasi_nodal_Bs_samples = copy(quasi_nodal_Bs_samples)
+# c_quasi_nodal_Bs_samples[3] = 0
+c_quasi_nodal_Bs_samples[3] /= 10
+c_quasi_nodal_Bs = Vs \ c_quasi_nodal_Bs_samples
+c_quasi_nodal_Bs = max.(c_quasi_nodal_Bs, 0.0)
+c_quasi_nodal_Bs_samples = Vs * c_quasi_nodal_Bs
+fn_c_quasi_nodal_Bs = intrpmat_Bs * c_quasi_nodal_Bs
+
+cc_quasi_nodal_Bs = copy(quasi_nodal_Bs)
+cc_quasi_nodal_Bs[3] = c_quasi_nodal_Bs[3]
+fn_cc_quasi_nodal_Bs = intrpmat_Bs * cc_quasi_nodal_Bs
 
 # corrected nodal Bspline
 c_nodal_Bs = max.(nodal_Bs, 0.0)
-fn_c_nodal_Bs = intrp_Bs * c_nodal_Bs
+fn_c_nodal_Bs = intrpmat_Bs * c_nodal_Bs
 
 # least-square Bspline
 lsqr_Bs = A_lsqr \ b_lsqr
-fn_lsqr_Bs = intrp_Bs * lsqr_Bs
+fn_lsqr_Bs = intrpmat_Bs * lsqr_Bs
 
 # corrected least-sqaure Bspline
 c_lsqr_Bs = max.(lsqr_Bs, 0.0)
-fn_c_lsqr_Bs = intrp_Bs * c_lsqr_Bs
+fn_c_lsqr_Bs = intrpmat_Bs * c_lsqr_Bs
 
 # nodal Lagrange
 intrp_La = dg1d.lagrange_interp_matrix(xx, xs)
@@ -132,14 +171,27 @@ scatterlines!(ax, xs, fx, color=:orange, markersize=22, 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)
+scatter!(ax, knots, quasi_nodal_Bs_samples, color=:brown, markersize=10,
+              label="quasi-nodal Bspline samples")
+scatter!(ax, knots, c_quasi_nodal_Bs_samples, color=:cyan, markersize=8,
+              label="corrected quasi-nodal Bspline samples")
+lines!(ax, xx, fn_c_quasi_nodal_Bs, color=:magenta, label="corrected quasi-nodal Bspline")
+scatterlines!(ax, knots, c_quasi_nodal_Bs, color=:magenta, markersize=6, 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_cc_quasi_nodal_Bs, color=:limegreen, label="a posteriori corrected quasi-nodal Bspline")
+scatterlines!(ax, knots, cc_quasi_nodal_Bs, color=:limegreen, markersize=6, 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)
+lines!(ax, xx, fn_c_nodal_Bs, color=:red, label="corrected nodal Bspline")
+scatterlines!(ax, knots, c_nodal_Bs, color=:red, markersize=6, 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)
+
+end # if
+end
diff --git a/src/ScalarEq/setup.jl b/src/ScalarEq/setup.jl
index 7e85208e48a900fb8ecbf442ecc736491f85eab9..3f906e69283abeb977c2f9e16b40d41a8a53ccfb 100644
--- a/src/ScalarEq/setup.jl
+++ b/src/ScalarEq/setup.jl
@@ -15,7 +15,7 @@ function Project(env::Environment, mesh::Mesh1d, prms)
   analyze_error_norm = prms["ScalarEq"]["analyze_error_norm"]
   V_bernstein = dg1d.Bernstein.vandermonde_matrix(mesh.element.z)
   if mesh.element isa DGElement && mesh.element.kind === :modal_bspline2
-    V_bspline = dg1d.Bspline.vandermonde_matrix(mesh.element.z, mesh.element.bs)
+    V_bspline = dg1d.Bspline.vandermonde_matrix(mesh.element.z, mesh.element.bs2)
     invV_bspline = pinv(V_bspline)
   else
     V_bspline, invV_bspline = nothing, nothing
diff --git a/src/bspline.jl b/src/bspline.jl
index 281336dc93e15b97bf7622234b1521fd436b2006..96e5dd1c6c20d1d737acd1a5f914a032734a97ed 100644
--- a/src/bspline.jl
+++ b/src/bspline.jl
@@ -5,6 +5,25 @@ using dg1d
 using Jacobi
 
 
+struct Bspline1
+  Npts::Int64 # number of unique knots
+  xs::Vector{Float64} # the knots, start and end point have multiplicity three
+end
+
+
+function Bspline1(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, z[end])
+  return Bspline1(Npts, xs)
+end
+
+
 struct Bspline2
   Npts::Int64 # number of unique knots
   xs::Vector{Float64} # the knots, start and end point have multiplicity three
@@ -24,6 +43,31 @@ function Bspline2(z::Vector{Float64})
 end
 
 
+function polynomial(x::Float64, i::Int64, bs::Bspline1)
+  @unpack Npts, xs = bs
+  1 ≤ i ≤ Npts || return 0.0
+  @inbounds begin
+    xi, xi1, xi2 = xs[i], xs[i+1], xs[i+2]
+  end
+  if x < xi || x > xi2
+    return 0.0
+  end
+  if i == 1 && x == xi1
+    # compute value on lhs bdry from the right
+    @goto sub2
+  end
+  if x ≤ xi1
+    return (x-xi)/(xi1-xi)
+  end
+  if x ≤ xi2
+    @label sub2
+    return (xi2-x)/(xi2-xi1)
+  end
+  # unreachable
+  return NaN
+end
+
+
 function polynomial(x::Float64, i::Int64, bs::Bspline2)
   @unpack Npts, xs = bs
   1 ≤ i ≤ Npts+1 || return 0.0
@@ -34,11 +78,11 @@ function polynomial(x::Float64, i::Int64, bs::Bspline2)
     return 0.0
   end
   if i == 1 && x == xi2
-    # compute derivative on lhs bdry from the right
+    # compute value on lhs bdry from the right
     @goto sub3
   end
   if i == 2 && x == xi1
-    # compute derivative on lhs bdry from the right
+    # compute value on lhs bdry from the right
     @goto sub2
   end
   if x ≤ xi1
@@ -57,6 +101,31 @@ function polynomial(x::Float64, i::Int64, bs::Bspline2)
 end
 
 
+function derivative_polynomial(x::Float64, i::Int64, bs::Bspline1)
+  @unpack Npts, xs = bs
+  1 ≤ i ≤ Npts || return 0.0
+  @inbounds begin
+    xi, xi1, xi2 = xs[i], xs[i+1], xs[i+2]
+  end
+  if x < xi || x > xi2
+    return 0.0
+  end
+  if i == 1 && x == xi1
+    # compute derivative on lhs bdry from the right
+    @goto sub2
+  end
+  if x ≤ xi1
+    return 1/(xi1-xi)
+  end
+  if x ≤ xi2
+    @label sub2
+    return -1/(xi2-xi1)
+  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
@@ -90,6 +159,40 @@ function derivative_polynomial(x::Float64, i::Int64, bs::Bspline2)
 end
 
 
+"""
+  interpolate(x::Real, fn::AbstractVector, bs::Bspline1)
+
+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::Bspline1)
+  @unpack Npts, xs = bs
+  if length(coeffs) != Npts
+    error("expected Npts = $Npts coefficients, received $(length(coeffs))")
+  end
+  return unsafe_interpolate(x, coeffs, bs)
+end
+function unsafe_interpolate(x::Real, coeffs::AbstractVector, bs::Bspline1)
+  @unpack Npts, xs = bs
+  @inbounds xl, xr = xs[1], xs[end]
+  xl == x && return coeffs[1]
+  xr == x && return coeffs[end]
+  xl < x < xr || return 0.0
+  o = 1 # offset duplicated bdry knots
+  v_xs = view(xs, o .+ (1:Npts))
+  ir = searchsortedfirst(v_xs, x)
+  il = ir-1
+  is_supp = max(il-1,1):min(ir+1,Npts)
+  B = 0.0
+  @inbounds for i in is_supp
+    B += coeffs[i] * polynomial(x, i, bs)
+  end
+  return B
+end
+
+
 """
   interpolate(x::Real, fn::AbstractVector, bs::Bspline2)
 
@@ -124,6 +227,34 @@ function unsafe_interpolate(x::Real, coeffs::AbstractVector, bs::Bspline2)
 end
 
 
+"""
+  mass_matrix(bs::Bspline1)
+
+Compute `M_ij = (Bi, Bj)` exactly where `Bi, Bj` are `i`th and `j`th Bspline of order 1.
+"""
+function mass_matrix(bs::Bspline1)
+  @unpack Npts, xs = bs
+  M = zeros(Float64, Npts, Npts)
+  zs, ws = dg1d.LG.rule(3) # 4th order accurate
+  o = 1 # 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, j in 1:Npts
+      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
+
+
 """
   mass_matrix(bs::Bspline2)
 
@@ -132,7 +263,7 @@ Compute `M_ij = (Bi, Bj)` exactly where `Bi, Bj` are `i`th and `j`th Bspline of
 function mass_matrix(bs::Bspline2)
   @unpack Npts, xs = bs
   M = zeros(Float64, Npts+1, Npts+1)
-  zs, ws = dg1d.LG.rule(2)
+  zs, ws = dg1d.LG.rule(3) # 4th order accurate
   o = 2 # offset bdry knots
   for k in 1:Npts-1
     il, ir = o+k, o+k+1
@@ -152,6 +283,34 @@ function mass_matrix(bs::Bspline2)
 end
 
 
+"""
+  stiffness_matrix(bs::Bspline1)
+
+Compute `S_ij = (Bi, dBj/dx)` exactly where `Bi, Bj` are `i`th and `j`th Bspline of order 1.
+"""
+function stiffness_matrix(bs::Bspline1)
+  @unpack Npts, xs = bs
+  S = zeros(Float64, Npts, Npts)
+  zs, ws = dg1d.LG.rule(3) # 4th order accurate
+  o = 1 # 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, j in 1:Npts
+      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
+
+
 """
   stiffness_matrix(bs::Bspline2)
 
@@ -160,7 +319,7 @@ Compute `S_ij = (Bi, dBj/dx)` exactly where `Bi, Bj` are `i`th and `j`th Bspline
 function stiffness_matrix(bs::Bspline2)
   @unpack Npts, xs = bs
   S = zeros(Float64, Npts+1, Npts+1)
-  zs, ws = dg1d.LG.rule(3)
+  zs, ws = dg1d.LG.rule(3) # 4th order accurate
   o = 2 # offset bdry knots
   for k in 1:Npts-1
     il, ir = o+k, o+k+1
@@ -180,11 +339,36 @@ function stiffness_matrix(bs::Bspline2)
 end
 
 
+function vandermonde_matrix(z::AbstractVector, bs::Bspline1)
+  return [ Bspline.polynomial(z[i], j, bs) for i in 1:length(z), j in 1:bs.Npts ]
+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::Bspline1)
+  @toggled_assert length(f) == bs.Npts
+  sum = 0.0
+  @inbounds for i in bs.Npts
+    fi, xl, xr = f[i], bs.xs[i], bs.xs[i+2]
+    sum += fi * (xr-xl)/2
+  end
+  return sum
+end
+function integrate(f::Function, idxs::AbstractVector{<:Integer}, bs::Bspline1)
+  @toggled_assert length(idxs) == bs.Npts
+  sum = 0.0
+  @inbounds for i in bs.Npts
+    fi, xl, xr = f[idxs[i]], bs.xs[i], bs.xs[i+2]
+    sum += fi * (xr-xl)/2
+  end
+  return sum
+end
+
+
 # analytical integration formula from
 # Splines and PDEs: From Approximation Theory to Numerical Linear Algebra 2017, p. 13
 function integrate(f::AbstractVector, bs::Bspline2)
diff --git a/src/dgelement.jl b/src/dgelement.jl
index b5698cf43ae34bedfaf5100f0f02089ae22b3058..3cd707001f81dbe0aed61b84da9cffac751cd753 100644
--- a/src/dgelement.jl
+++ b/src/dgelement.jl
@@ -29,7 +29,8 @@ 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
+  bs1::Bspline.Bspline1    # only set for :modal_bspline1
+  bs2::Bspline.Bspline2    # only set for :modal_bspline2
 
   # TODO Replace N with Npts and remove N where ever possible.
   function DGElement(N::Integer, quadr::Symbol, kind::Symbol;
@@ -39,13 +40,13 @@ struct DGElement
     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.")
+    if kind ∉ (:nodal_lagrange, :modal_bernstein, :modal_bspline2, :modal_bspline1)
+      error("Invalid DG formulation specified: $kind. Use one of nodal_lagrange, modal_bernstein, modal_bspline2, modal_bspline1.")
     end
 
     Npts = N + 1
 
-    local bs
+    local bs1, bs2
     if kind === :nodal_lagrange
       if quadr === :lgl
         quadr_z, quadr_w = LGL.rule(Npts)
@@ -91,6 +92,33 @@ struct DGElement
       l_rhs      = zeros(Float64, Npts)
       l_lhs[1]   = 1.0
       l_rhs[end] = 1.0
+    elseif kind === :modal_bspline1
+      if N < 1
+        error("Requested N+1=$(N+1) DOFs for DG formulation of kind $kind, but require at least 2.")
+      end
+      # User asked for Nth order polynomial sampled on Npts=N+1 points,
+      # we supply Npts 1st order Bsplines -- one for each DOF.
+      # bs_z is now the grid across which the Bsplines are C^0 continuous.
+      if quadr === :lgl || quadr === :lgl_lumped
+        z, _ = LGL.rule(Npts)
+      elseif quadr === :glgl
+        z, _ = GLGL.rule(Npts)
+      end
+      purge_zeros!(z)
+      bs1 = Bspline.Bspline1(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(3)
+      D = [ Bspline.derivative_polynomial(z[i], j, bs1) for i in 1:Npts, j in 1:Npts ]
+      M = Bspline.mass_matrix(bs1)
+      S = Bspline.stiffness_matrix(bs1)
+      invM = inv(M)
+      V = Bspline.vandermonde_matrix(z, bs1)
+      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.")
@@ -105,7 +133,7 @@ struct DGElement
         bs_z, _ = GLGL.rule(N)
       end
       purge_zeros!(bs_z)
-      bs = Bspline.Bspline2(bs_z)
+      bs2 = 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])
@@ -113,21 +141,24 @@ struct DGElement
       # 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)
+      quadr_z, quadr_w = LG.rule(3)
+      D = [ Bspline.derivative_polynomial(z[i], j, bs2) for i in 1:Npts, j in 1:Npts ]
+      M = Bspline.mass_matrix(bs2)
+      S = Bspline.stiffness_matrix(bs2)
       invM = inv(M)
-      V = Bspline.vandermonde_matrix(z, bs)
+      V = Bspline.vandermonde_matrix(z, bs2)
       l_lhs      = zeros(Float64, Npts)
       l_rhs      = zeros(Float64, Npts)
       l_lhs[1]   = 1.0
       l_rhs[end] = 1.0
     end
 
+    # just set up something, its unused for the other kinds
+    if kind !== :modal_bspline1
+      bs1 = Bspline.Bspline1(z)
+    end
     if kind !== :modal_bspline2
-      # just set up something, its unused for the other kinds
-      bs = Bspline.Bspline2(z)
+      bs2 = Bspline.Bspline2(z)
     end
 
     # TODO Rename MDM to MST or so, because the identity MDM = M^{-1} S^{T} only applies for nodal DG
@@ -152,7 +183,7 @@ struct DGElement
     return new(N, Npts, kind,
                quadr, quadr_z, quadr_w,
                z, D, S, M, invM, MDM, Ml_lhs, Ml_rhs,
-               V, invV, F, bs)
+               V, invV, F, bs1, bs2)
   end
 end
 
@@ -179,7 +210,7 @@ end
 
 
 @inline function integrate(u::AbstractVector, el::DGElement)
-  if el.kind === :modal_bspline2
+  if el.kind ∈ (:modal_bspline2, :modal_bspline1)
     return Bspline.integrate(u, el.bs)
   else
     if el.quadr === :lgl || el.quadr === :lgl_lumped
@@ -192,7 +223,7 @@ end
 
 
 @inline function integrate(fn::Function, idxs::AbstractVector{<:Integer}, el::DGElement)
-  if el.kind === :modal_bspline2
+  if el.kind ∈ (:modal_bspline2, :modal_bspline1)
     return Bspline.integrate(fn, idxs, el.bs)
   else
     if el.quadr === :lgl || el.quadr === :lgl_lumped
diff --git a/src/parameters.jl b/src/parameters.jl
index 2160f97ce7ee979b46a3ab44f5b0821350e70c0b..da44b1818c9c5989afdd81ebcc149c7cf390361f 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -516,6 +516,15 @@ end
     but the Bernstein polynomials are not interpolating.
     We refer to these coefficients as being quasi-nodal.
     See [2,3,4].
+  - `modal_bspline1`: A modal approximation using 1st order Bspline polynomials.
+    The coefficients correspond to the nodal values sampled on the requested quadrature grid.
+    I.e. a `modal_bspline1` 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.
   - `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`.
@@ -533,7 +542,7 @@ end
   [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" ]
+  @check scheme in [ "nodal_lagrange", "modal_bernstein", "modal_bspline1", "modal_bspline2" ]
 
   """
   Periodicity in cartesian directions `x,y`
diff --git a/test/IntegrationTests/refs/modal_bspline2_advection_sine/output1d.h5 b/test/IntegrationTests/refs/modal_bspline2_advection_sine/output1d.h5
index d0e6758edad6c361ad5196a967dd777fb917170e..d960106c40355565a66ad8120a2c0e35f791f645 100644
Binary files a/test/IntegrationTests/refs/modal_bspline2_advection_sine/output1d.h5 and b/test/IntegrationTests/refs/modal_bspline2_advection_sine/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/modal_bspline2_advection_sine/substep_output.h5 b/test/IntegrationTests/refs/modal_bspline2_advection_sine/substep_output.h5
index 14ee01b05bc913a144ae5617141275c5929794b4..02c0423ec92d2fb84d4e9102a9cb5c5626e1e40d 100644
Binary files a/test/IntegrationTests/refs/modal_bspline2_advection_sine/substep_output.h5 and b/test/IntegrationTests/refs/modal_bspline2_advection_sine/substep_output.h5 differ
diff --git a/test/UnitTests/src/test_dgelement.jl b/test/UnitTests/src/test_dgelement.jl
index 210ac1997eac5f04dbe922fef402bbf98c612151..c14d817bd7db66255b8495c73aca017f3134b983 100644
--- a/test/UnitTests/src/test_dgelement.jl
+++ b/test/UnitTests/src/test_dgelement.jl
@@ -1,4 +1,4 @@
-@testset "Spectral Element" begin
+@testset "DGElement" begin
 
 
 function evaluate_derivatives(el::DGElement, f, df)
@@ -202,7 +202,7 @@ end
   @test all(el.S.≈TST)
 
   el = dg1d.DGElement(N, :lgl, :modal_bspline2)
-  bs = el.bs
+  bs = el.bs2
   M = el.M
   S = el.S
   ST = transpose(S)
@@ -221,6 +221,26 @@ end
   dg1d.purge_zeros!(exact_BB)
   @test all(exact_BB .≈ BB)
 
+  el = dg1d.DGElement(N, :lgl, :modal_bspline1)
+  bs = el.bs1
+  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 = 1
+  for i in 1:bs.Npts, j in 1:bs.Npts
+    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
 
 
@@ -267,6 +287,20 @@ end
   POU = dg1d.Bspline.interpolate.(xx, Ref(coeffs), Ref(bs))
   @test all(POU .≈ ones(Float64, length(POU)))
 
+
+  z, _ = dg1d.LGL.rule(5)
+  bs = dg1d.Bspline.Bspline1(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)
+  POU = dg1d.Bspline.interpolate.(xx, Ref(coeffs), Ref(bs))
+  @test all(POU .≈ ones(Float64, length(POU)))
+
 end