diff --git a/src/EulerEq/rhs.jl b/src/EulerEq/rhs.jl
index 6920c55f7715315a45b939f3f608f9b25e67a0fb..2d78894584c1da9cec7a6b89709dca14098cfde0 100644
--- a/src/EulerEq/rhs.jl
+++ b/src/EulerEq/rhs.jl
@@ -21,7 +21,7 @@ function timestep(env, P::Project, hrsc::Maybe{HRSC.AbstractHRSC}, mesh::Mesh1d)
   return dt
 end
 dtfactor(mesh::Mesh1d{FVElement}) = 2
-dtfactor(mesh::Mesh1d{SpectralElement}) = mesh.element.N
+dtfactor(mesh::Mesh1d{DGElement}) = mesh.element.N
 
 
 function timestep(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, ::Mesh1d)
diff --git a/src/GRHD/callbacks.jl b/src/GRHD/callbacks.jl
index 9192c61ca4275d44e2846bfc147fc9966657ca0e..66cf7115fb40a6568308a8ac4f8203c53e894e25 100644
--- a/src/GRHD/callbacks.jl
+++ b/src/GRHD/callbacks.jl
@@ -48,7 +48,7 @@ function callback_equation(state_u, state_t, env, P)
 end
 
 compute_cellmax_v(P, mesh) = nothing
-function compute_cellmax_v(P, mesh::Mesh1d{SpectralElement})
+function compute_cellmax_v(P, mesh::Mesh1d{DGElement})
   @unpack max_v = get_static_variables(mesh)
   @unpack cellmax_v = get_cell_variables(mesh)
   for (k,vi) in enumerate(eachcell(mesh, max_v))
@@ -57,7 +57,7 @@ function compute_cellmax_v(P, mesh::Mesh1d{SpectralElement})
 end
 
 compute_atmosphere_mask(env, P, mesh) = nothing
-function compute_atmosphere_mask(env, P, mesh::Mesh1d{SpectralElement})
+function compute_atmosphere_mask(env, P, mesh::Mesh1d{DGElement})
   enforce        = P.prms.c2p_enforce_atm
   enforce_causal = P.prms.c2p_enforce_causal_atm
   (!enforce && !enforce_causal) && return
diff --git a/src/GRHD/cartoon.jl b/src/GRHD/cartoon.jl
index 37bad364ba27291f3ab231e169110cb754010bc1..26b97b3fbd13c2d1d8d04af7abbb419ba6c92ead 100644
--- a/src/GRHD/cartoon.jl
+++ b/src/GRHD/cartoon.jl
@@ -368,7 +368,7 @@ end
 end
 
 
-function impose_symmetry_bdry_values_cartoon!(mesh::Mesh2d{SpectralElement})
+function impose_symmetry_bdry_values_cartoon!(mesh::Mesh2d{DGElement})
   @unpack cache = mesh
   @unpack D, Sx, Sz, Ï„, x, y, r = get_static_variables(mesh)
   @unpack bdry_D, bdry_Sx, bdry_Sz, bdry_Ï„, nx, ny = get_bdry_variables(mesh)
diff --git a/src/GRHD/doublecartoon.jl b/src/GRHD/doublecartoon.jl
index 31caea98dc4692206adc73d7febc89c802bc828d..f03ca5b817d6b7243d766cdceaa76bb3b06b0628 100644
--- a/src/GRHD/doublecartoon.jl
+++ b/src/GRHD/doublecartoon.jl
@@ -308,7 +308,7 @@ end
 end
 
 
-function impose_symmetry_bdry_values_doublecartoon!(mesh::Mesh1d{SpectralElement})
+function impose_symmetry_bdry_values_doublecartoon!(mesh::Mesh1d{DGElement})
 
   @unpack cache = mesh
   @unpack D, Sx, Ï„ = get_static_variables(mesh)
diff --git a/src/GRHD/rhs.jl b/src/GRHD/rhs.jl
index abe5a67b080438695720fe7fd0a2285a7b52ebb3..09187952e0de799664d0e55680354dab50626cad 100644
--- a/src/GRHD/rhs.jl
+++ b/src/GRHD/rhs.jl
@@ -39,9 +39,9 @@ function timestep(mesh, P::Project, hrsc::Maybe{HRSC.AbstractHRSC})
   return dt
 end
 dtfactor(mesh::Mesh1d{FVElement}) = 2
-dtfactor(mesh::Mesh1d{SpectralElement}) = 1
+dtfactor(mesh::Mesh1d{DGElement}) = 1
 
-function timestep(mesh::Mesh2d{SpectralElement}, P::Project, hrsc::Maybe{HRSC.AbstractHRSC})
+function timestep(mesh::Mesh2d{DGElement}, P::Project, hrsc::Maybe{HRSC.AbstractHRSC})
   @unpack element = mesh
   @unpack N, z = element
   compute_maxspeed(P, P.equation, mesh)
diff --git a/src/GRHD/spherical1d.jl b/src/GRHD/spherical1d.jl
index ce0bfd81c01a6f52e861f86be485f74602d8a8b1..9f262ae9f1751b91df9fb0d1e650b44831de1111 100644
--- a/src/GRHD/spherical1d.jl
+++ b/src/GRHD/spherical1d.jl
@@ -271,7 +271,7 @@ end
 end
 
 
-function impose_symmetry_bdry_values_spherical1d!(mesh::Mesh1d{SpectralElement})
+function impose_symmetry_bdry_values_spherical1d!(mesh::Mesh1d{DGElement})
 
   @unpack cache = mesh
   @unpack D, Sr, Ï„ = get_dynamic_variables(mesh)
diff --git a/src/GRHD/valencia1d.jl b/src/GRHD/valencia1d.jl
index 9a7153b46afaba647955b7acb693c3642175c9ef..4e2ba921382c77ab9c13ced83026c2eb7c95c308 100644
--- a/src/GRHD/valencia1d.jl
+++ b/src/GRHD/valencia1d.jl
@@ -275,7 +275,7 @@ end
 end
 
 
-function impose_symmetry_bdry_values_valencia1d!(mesh::Mesh1d{SpectralElement})
+function impose_symmetry_bdry_values_valencia1d!(mesh::Mesh1d{DGElement})
 
   @unpack cache = mesh
   @unpack D, Sr, Ï„ = get_static_variables(mesh)
diff --git a/src/HRSC/bernstein.jl b/src/HRSC/bernstein.jl
index 22c4cfca3d665ef2d988de281acd65dcd2da155e..d8a6a1d0ffa6b72d31beec96aa4f37b2554c3b06 100644
--- a/src/HRSC/bernstein.jl
+++ b/src/HRSC/bernstein.jl
@@ -28,7 +28,7 @@ function BernsteinReconstruction(mesh::Mesh)
   rng = collect(0:Npts-1)
   bernstein_z = @. -1 + 2 * rng / (Npts-1)
   P = dg1d.lagrange_interp_matrix(bernstein_z, z)
-  V = vandermonde_matrix_bernstein(z)
+  V = dg1d.Bernstein.vandermonde_matrix(z)
   return BernsteinReconstruction(mesh, P, V)
 end
 
@@ -38,20 +38,6 @@ end
 #######################################################################
 
 
-# Reference contains a typo in the normalization of the polynomials,
-# it should be n instead of N.
-function bernstein_polynomial(x, n, N)
-  return binomial(N, n) / 2^N * (1 + x)^n * (1 - x)^(N-n)
-end
-
-
-function vandermonde_matrix_bernstein(z)
-  Npts = length(z)
-  V = [ bernstein_polynomial(z[i], j-1, Npts-1) for i=1:Npts, j=1:Npts ]
-  return V
-end
-
-
 @doc """
     reconstruct_cell!(u!, w, P, V; kwargs...)
 
diff --git a/src/HRSC/weno.jl b/src/HRSC/weno.jl
index cb9d2eb6176b7f7bf54b445f29619df65581d315..cf506b231b11a6d2b40f6fb7abf744e1d6b3de73 100644
--- a/src/HRSC/weno.jl
+++ b/src/HRSC/weno.jl
@@ -122,7 +122,7 @@ function smoothness_estimator_matrix(w)
   @assert w <= 2 "Smoothness estimator only implemented up to w = 2"
 
   # setup a local quadrature rule to evaluate integrals over derivatives of lagrange polynomials
-  element = SpectralElement(2*w, :LGL)
+  element = DGElement(2*w, :LGL)
   @unpack iproduct, z, Npts = element
 
   # build w+1 subcell grids which all overlap with x
diff --git a/src/SRHD/rhs.jl b/src/SRHD/rhs.jl
index 87453de5385cb79f5e8d01757a32766339f457a7..271cf2ff5765e0d6662cd3fb97e5f43de830731a 100644
--- a/src/SRHD/rhs.jl
+++ b/src/SRHD/rhs.jl
@@ -27,7 +27,7 @@ function timestep(env, P::Project, hrsc::Maybe{HRSC.AbstractHRSC}, mesh::Mesh1d)
   return dt
 end
 dtfactor(mesh::Mesh1d{FVElement}) = 2
-dtfactor(mesh::Mesh1d{SpectralElement}) = mesh.element.N
+dtfactor(mesh::Mesh1d{DGElement}) = mesh.element.N
 
 
 function timestep(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, mesh::Mesh1d)
diff --git a/src/ScalarEq/rhs.jl b/src/ScalarEq/rhs.jl
index 756811d0cb7222a5366fca46a87f0bf7419e2eba..75100c42b3a8369bd1470fe804dc7319af3a2d23 100644
--- a/src/ScalarEq/rhs.jl
+++ b/src/ScalarEq/rhs.jl
@@ -25,7 +25,7 @@ function timestep(env, P::Project, mesh::Mesh1d, hrsc::Maybe{HRSC.AbstractHRSC})
   return dt
 end
 dtfactor(mesh::Mesh1d{FVElement}) = 2
-dtfactor(mesh::Mesh1d{SpectralElement}) = mesh.element.N
+dtfactor(mesh::Mesh1d{DGElement}) = mesh.element.N
 
 
 function timestep(env, P::Project, mesh::Mesh1d, hrsc::HRSC.AbstractArtificialViscosity)
@@ -50,7 +50,7 @@ function timestep(env, P::Project, mesh::Mesh1d, hrsc::HRSC.AbstractArtificialVi
 end
 
 
-function timestep(env, P::Project, ::Mesh2d{SpectralElement}, hrsc::Maybe{HRSC.AbstractHRSC})
+function timestep(env, P::Project, ::Mesh2d{DGElement}, hrsc::Maybe{HRSC.AbstractHRSC})
   @unpack cache, mesh = env
   @unpack equation = P
 
diff --git a/src/bernstein.jl b/src/bernstein.jl
new file mode 100644
index 0000000000000000000000000000000000000000..5ccaec4a24f36133fca082833f56308780eb884c
--- /dev/null
+++ b/src/bernstein.jl
@@ -0,0 +1,127 @@
+module Bernstein
+
+
+using dg1d
+using LinearAlgebra
+using Jacobi
+
+
+# Normalized so that bernstein_polynomial(-1,0,N) == bernstein_polynomial(1,N,N) == 1
+function bernstein_polynomial(x::Real, n::Integer, N::Integer)
+  return binomial(N, n) / 2^N * (1 + x)^n * (1 - x)^(N-n)
+end
+
+
+function derivative_bernstein_polynomial(x::Real, n::Integer, N::Integer)
+  B = 0.0
+  if n > 0
+    B += n*(1+x)^(n-1) * (1-x)^(N-n)
+  end
+  if n < N
+    B += (1+x)^n * (N-n)*(1-x)^(N-n-1)*(-1)
+  end
+  return binomial(N, n) / 2^N * B
+end
+
+
+function vandermonde_matrix(z::AbstractVector)
+  Npts = length(z)
+  V = [ bernstein_polynomial(z[i], j-1, Npts-1) for i=1:Npts, j=1:Npts ]
+  return V
+end
+
+
+"""
+  mass_matrix(Npts::Integer)
+
+Compute M_ij = (Bi, Bj) exactly where `Bi, Bj` are (i-1)th and (j-1)th Bernstein polynomials.
+"""
+function mass_matrix(Npts::Integer)
+  N = Npts-1
+  M = zeros(Float64, Npts, Npts)
+  xs, ws = dg1d.LGL.rule(Npts+1)
+  for i in 1:Npts, j in i:Npts
+    mij = dg1d.LGL.integrate(xs, ws) do xx
+      bi = bernstein_polynomial(xx, i-1, N)
+      bj = bernstein_polynomial(xx, j-1, N)
+      return bi*bj
+    end
+    M[i,j] = mij
+    if j > i
+      M[j,i] = mij
+    end
+  end
+  return M
+end
+
+
+"""
+  stiffness_matrix(Npts::Integer)
+
+Compute S_ij = (Bi, dBj/dx) exactly where `Bi, Bj` are (i-1)th and (j-1)th Bernstein polynomials.
+"""
+function stiffness_matrix(Npts::Integer)
+  N = Npts-1
+  S = zeros(Float64, Npts, Npts)
+  xs, ws = dg1d.LGL.rule(Npts+1)
+  for i in 1:Npts, j in 1:Npts
+    sij = dg1d.LGL.integrate(xs, ws) do xx
+      bi = bernstein_polynomial(xx, i-1, N)
+      ∂bj = derivative_bernstein_polynomial(xx, j-1, N)
+      return bi*∂bj
+    end
+    S[i,j] = sij
+  end
+  return S
+end
+
+
+"""
+  bernstein_legendre_transformation(Npts::Integer) -> (T, invT)
+
+Compute a transformation matrix `T` and its inverse `invT` to map
+quasi-nodal Bernstein coefficients `Bn` to modal Legendre coefficients `Ln`
+by `Bn = T*Ln`.
+"""
+function bernstein_legendre_transformation(Npts::Integer)
+  N = Npts-1
+  xs, ws = dg1d.LGL.rule(Npts+1)
+  BP = (i,j) -> dg1d.LGL.integrate(xs, ws) do xi
+    Bim1 = bernstein_polynomial(xi, i-1, N)
+    Pjm1 = Jacobi.legendre(xi, j-1)
+    Bim1*Pjm1
+  end
+  PP = i -> dg1d.LGL.integrate(xs, ws) do xi
+    Pim1 = Jacobi.legendre(xi, i-1)
+    Pim1*Pim1
+  end
+
+  T = zeros(Float64, Npts, Npts)
+  for i in 1:Npts, j in 1:Npts
+    T[i,j] = BP(i,j)/PP(j)
+  end
+  invT = inv(T)
+
+  return T, invT
+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`].
+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, fn::AbstractVector)
+  Npts = length(fn)
+  N = Npts-1
+  return sum(enumerate(fn)) do (i,fi)
+    Bim1 = bernstein_polynomial(x, i-1, N)
+    fi*Bim1
+  end
+end
+
+
+end
diff --git a/src/dg1d.jl b/src/dg1d.jl
index a881e5dc222f63ae3b9c21198493ca5314f9798a..39ac3bfb894a3eb5876c03cd7115ad28f02e1af0 100644
--- a/src/dg1d.jl
+++ b/src/dg1d.jl
@@ -51,13 +51,15 @@ include("output.jl")
 include("lagrange.jl")
 include("glgl.jl")
 include("lgl.jl")
+include("bernstein.jl")
+include("legendre.jl")
 
 #######################################################################
 #                             Components                              #
 #######################################################################
 
-export SpectralElement
-include("spectralelement.jl")
+export DGElement
+include("dgelement.jl")
 
 export FVElement
 include("fvelement.jl")
@@ -112,7 +114,7 @@ export fv_update_step!
 include("fv_rhs.jl")
 
 # TODO Move this to top; also requires gather all type defintions
-# like SpectralElement, AbstractMesh, Tree etc. into types.jl
+# like DGElement, AbstractMesh, Tree etc. into types.jl
 export Maybe, Environment
 include("types.jl")
 
diff --git a/src/dg_rhs.jl b/src/dg_rhs.jl
index d45c126efb4bdaa2f5f5c7dc04ae5ed9ce8b8fb7..73b81a958cfae0d5a39bd0ae0e657ef124fcc2cc 100644
--- a/src/dg_rhs.jl
+++ b/src/dg_rhs.jl
@@ -1,11 +1,11 @@
 """
-    compute_rhs_weak_form!(rhs, f, nf, mesh::Mesh1d{SpectralElement})
-    compute_rhs_weak_form!(rhs, f, s, nf, mesh::Mesh1d{SpectralElement})
+    compute_rhs_weak_form!(rhs, f, nf, mesh::Mesh1d{DGElement})
+    compute_rhs_weak_form!(rhs, f, s, nf, mesh::Mesh1d{DGElement})
 
 Compute the `rhs` of the weak DG formulation using
 the (bulk) flux `f`, the numerical flux `nf` and sources `s` for a `mesh`.
 """
-function compute_rhs_weak_form!(rhs, f, nf, mesh::Mesh1d{SpectralElement})
+function compute_rhs_weak_form!(rhs, f, nf, mesh::Mesh1d{DGElement})
   @unpack element = mesh
   @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
   # @unpack K = mesh
@@ -24,7 +24,7 @@ function compute_rhs_weak_form!(rhs, f, nf, mesh::Mesh1d{SpectralElement})
   end
   return
 end
-function compute_rhs_weak_form!(rhs, f, s, nf, mesh::Mesh1d{SpectralElement})
+function compute_rhs_weak_form!(rhs, f, s, nf, mesh::Mesh1d{DGElement})
   compute_rhs_weak_form!(rhs, f, nf, mesh)
   rhs .+= s
   return
@@ -32,13 +32,13 @@ end
 
 
 """
-    compute_rhs_weak_form!(rhs, fx, fy, nf, mesh::Mesh2d{SpectralElement})
-    compute_rhs_weak_form!(rhs, fx, fy, s, nf, mesh::Mesh2d{SpectralElement})
+    compute_rhs_weak_form!(rhs, fx, fy, nf, mesh::Mesh2d{DGElement})
+    compute_rhs_weak_form!(rhs, fx, fy, s, nf, mesh::Mesh2d{DGElement})
 
 Compute the `rhs` of the weak DG formulation using
 the (bulk) flux `fx, fy`, the numerical flux `nf` and sources `s` for a `mesh`.
 """
-function compute_rhs_weak_form!(rhs, fx, fy, nf, mesh::Mesh2d{SpectralElement})
+function compute_rhs_weak_form!(rhs, fx, fy, nf, mesh::Mesh2d{DGElement})
   @unpack element = mesh
   @unpack dxdX, dydY = get_static_variables(mesh.cache)
   @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
@@ -109,13 +109,13 @@ function compute_rhs_weak_form!(rhs, fx, fy, nf, mesh::Mesh2d{SpectralElement})
 end
 
 """
-    compute_rhs_weak_form!(rhs, fx, _::Real, nf, mesh::Mesh2d{SpectralElement})
-    compute_rhs_weak_form!(rhs, _::Real, fy, nf, mesh::Mesh2d{SpectralElement})
+    compute_rhs_weak_form!(rhs, fx, _::Real, nf, mesh::Mesh2d{DGElement})
+    compute_rhs_weak_form!(rhs, _::Real, fy, nf, mesh::Mesh2d{DGElement})
 
 Specialized version of `compute_rhs_weak_form!` where only the fluxs
 `fx` or `fy` are applied, respectively.
 """
-function compute_rhs_weak_form!(rhs, fx, _::Real, nf, mesh::Mesh2d{SpectralElement})
+function compute_rhs_weak_form!(rhs, fx, _::Real, nf, mesh::Mesh2d{DGElement})
   @unpack element = mesh
   @unpack dxdX, dydY = get_static_variables(mesh.cache)
   @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
@@ -171,7 +171,7 @@ function compute_rhs_weak_form!(rhs, fx, _::Real, nf, mesh::Mesh2d{SpectralEleme
   end
   return
 end
-function compute_rhs_weak_form!(rhs, _::Real, fy, nf, mesh::Mesh2d{SpectralElement})
+function compute_rhs_weak_form!(rhs, _::Real, fy, nf, mesh::Mesh2d{DGElement})
   @unpack element = mesh
   @unpack dxdX, dydY = get_static_variables(mesh.cache)
   @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
@@ -227,7 +227,7 @@ function compute_rhs_weak_form!(rhs, _::Real, fy, nf, mesh::Mesh2d{SpectralEleme
   end
   return
 end
-function compute_rhs_weak_form!(rhs, fx, fy, s, nf, mesh::Mesh2d{SpectralElement})
+function compute_rhs_weak_form!(rhs, fx, fy, s, nf, mesh::Mesh2d{DGElement})
   compute_rhs_weak_form!(rhs, fx, fy, nf, mesh)
   rhs .+= s
   return
diff --git a/src/dgelement.jl b/src/dgelement.jl
new file mode 100644
index 0000000000000000000000000000000000000000..f9de0ac97e7a7b03d29d73398ca7b5c9f9cfcf9a
--- /dev/null
+++ b/src/dgelement.jl
@@ -0,0 +1,158 @@
+struct DGElement
+  N::Int64                # polynomial order
+  Npts::Int64             # number of quadrature points = N + 1
+  kind::Symbol            # kind of nodal approximation
+
+  # quadrature rule
+  # quadr_w are such that data sampled on `z` can be directly passed to the integrate
+  # functions below. In particular, any interpolation from `z` onto `quadr_z`
+  # (need kind === :nodal_bernstein) is already accounted for.
+  quadr::Symbol             # :lgl, :lgl_mass_lumped, :glgl
+  quadr_z::Vector{Float64}  # quadrature points
+  quadr_w::Vector{Float64}  # quadrature weights
+
+  # nodal points and DG operators
+  z::Vector{Float64}      # nodal points
+  D::Matrix{Float64}      # derivative matrix
+  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
+  Ml_lhs::Vector{Float64} # = M^-1 l(-1)
+  Ml_rhs::Vector{Float64} # = M^-1 l(+1)
+
+  # helpers
+  V::Matrix{Float64}      # = P_(j-1)(z_i), Vandermonde matrix
+  invV::Matrix{Float64}   # inv(V), inverse Vandermonde matrix
+  F::Matrix{Float64}      # spectral filtering matrix
+
+  function DGElement(N::Integer, quadr::Symbol, kind::Symbol;
+      # filter parameters
+      ηc::Float64=0.5, sh::Float64=16.0, α::Float64=36.0)
+
+    Npts = N + 1
+
+    if kind === :nodal_lagrange
+      if quadr === :lgl
+        quadr_z, quadr_w = LGL.rule(Npts)
+        purge_zeros!(quadr_z)
+        S    = LGL.stiffness_matrix(Npts)
+        M    = LGL.mass_matrix_exact(Npts)
+        invM = LGL.inverse_mass_matrix_exact(Npts)
+      elseif quadr === :lgl_lumped
+        quadr_z, quadr_w = LGL.rule(Npts)
+        purge_zeros!(quadr_z)
+        S    = LGL.stiffness_matrix(Npts)
+        M    = LGL.mass_matrix_lumped(Npts)
+        invM = LGL.inverse_mass_matrix_lumped(Npts)
+      elseif quadr === :glgl
+        quadr_z, quadr_w = GLGL.rule(Npts)
+        purge_zeros!(quadr_z)
+        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[1]   = 1.0
+      l_rhs[end] = 1.0
+    elseif kind === :nodal_bernstein
+      if quadr === :lgl || quadr === :lgl_lumped
+        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 ]
+      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[1]   = 1.0
+      l_rhs[end] = 1.0
+    else
+      error("Invalid quadrature rule specified")
+    end
+
+    MDM    = invM * transpose(S)
+    Ml_lhs = invM * l_lhs
+    Ml_rhs = invM * l_rhs
+    invV = inv(V)
+    purge_zeros!(MDM)
+    purge_zeros!(M)
+    purge_zeros!(S)
+    purge_zeros!(D)
+    purge_zeros!(invM)
+    purge_zeros!(V)
+    purge_zeros!(invV)
+    purge_zeros!(Ml_lhs)
+    purge_zeros!(Ml_rhs)
+    Λ = 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)
+  end
+end
+
+
+# (5.16) from Hesthaven, Warburton 2007
+function spectral_filter(η::Real, prms)
+  ηc, sh, α = prms
+  return if 0 ≤ η ≤ ηc
+    1
+  else
+    exp(-α*((η-ηc)/(1-ηc))^(2*sh))
+  end
+end
+
+
+@inline function inner_product(u1::AbstractVector, u2::AbstractVector, el::DGElement)
+  @unpack Npts, M = el
+  result = 0.0
+  @turbo for i in 1:Npts, j in 1:Npts
+    result += u1[i] * M[i,j] * u2[j]
+  end
+  return result
+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)
+  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)
+  end
+end
+
+
+@inline function differentiate(u::AbstractArray, el::DGElement)
+  return el.D * u
+end
+
+
+@inline function differentiate!(du::AbstractArray, u::AbstractArray, el::DGElement)
+  mul!(du, el.D, u)
+end
diff --git a/src/glgl.jl b/src/glgl.jl
index 54207a060c001df34c369f22da71452ddc3f9b9e..22ecdddf7d1127ebd97065e4c9b18a2494b44e6e 100644
--- a/src/glgl.jl
+++ b/src/glgl.jl
@@ -114,10 +114,10 @@ end
 
 
 """
-  mass_matrix(Npts::Integer)
+  mass_matrix_exact(Npts::Integer)
 
-Mass matrix for Lagrange polynomials based on GLGL quadrature nodes.
-M_ij = (li, lj)
+Compute M_ij = (li, lj) exactly where `li, lj` are Lagrange polynomials
+defined on the GLGL quadrature nodes.
 """
 function mass_matrix(Npts::Integer)
   x, w, v, D = inner_product_rule(Npts)
@@ -133,8 +133,8 @@ end
 """
   stiffness_matrix(Npts::Integer)
 
-Stiffness matrix for Lagrange polynomials based on GLGL quadrature nodes.
-S_ij = (li, dlj/dx)
+Compute S_ij = (li, dlj/dx) exactly where `li, lj` are Lagrange polynomials
+defined on the GLGL quadrature nodes.
 """
 function stiffness_matrix(Npts::Integer)
   x, w, v, D = inner_product_rule(Npts)
diff --git a/src/lagrange.jl b/src/lagrange.jl
index fb6b8fb1ba739bdd0bae52d14dc00d3a5f18f3dc..12ec2bc382ece81cb3e01d5a7386b633c3844b19 100644
--- a/src/lagrange.jl
+++ b/src/lagrange.jl
@@ -1,11 +1,11 @@
-@doc """
+"""
 Lagrange interpolation l_i(x, {z}) polynomials.
 Kopriva 2009, p. 74, (3.32)
 m ... polynomial idx
 x ... point on which to evaluate Lagrange polynomial
 z ... grid points on which Lagrange polynomials are defined, l_i(z_j) = delta_ij
 """
-function lagrange_poly(m, x, z)
+function lagrange_poly(m, x, z::AbstractVector)
   result = 1.0
   xm = z[m]
   for j in 1:length(z)
@@ -17,15 +17,37 @@ function lagrange_poly(m, x, z)
 end # function
 
 
+function derivative_lagrange_poly(m, x, z::AbstractVector)
+  Npts = length(z)
+  denom = 1.0
+  zm = z[m]
+  for j in 1:Npts
+    j == m && continue
+    denom *= z[j] - zm
+  end
+  result = 0.0
+  for j in 1:Npts
+    j == m && continue
+    r = 1.0
+    for k in 1:Npts
+      ((k == m) || (k == j)) && continue
+      r *= x - z[k]
+    end
+    result += r
+  end
+  result /= denom
+  return result
+end
+
 
-@doc """
-  `lagrange_interp_matrix(x, z)`
+"""
+  lagrange_interp_matrix(x, z)
 
 Interpolate data from reference grid `z` onto sample grid `x` by matrix multiplication.
 
 TODO Rename this to interpolation matrix.
 """
-function lagrange_interp_matrix(x, z)
+function lagrange_interp_matrix(x::AbstractVector, z::AbstractVector)
   N = length(x)
   M = length(z)
   T = Matrix{Float64}(undef, N, M)
@@ -37,12 +59,12 @@ end # function
 
 
 
-@doc """
+"""
 Needed for derivative matrices for Lagrange polynomials.
 Kopriva 2009, p. 74, (3.34).
 x ... Interpolation nodes.
 """
-function barycentric_weights(x)
+function barycentric_weights(x::AbstractVector)
   N = length(x)
 
   w = zeros(N)
@@ -62,10 +84,10 @@ end
 
 
 
-@doc """
+"""
 Kopriva 2009, p. 74, (3.36).
 """
-function barycentric_interp(x, z, w, f)
+function barycentric_interp(x, z::AbstractVector, w::AbstractVector, f::AbstractVector)
   idx = findfirst(zk->abs(zk - x) ≈ 0, z)
   !isnothing(idx) && return f[idx]
   t = @. w / (x - z)
@@ -75,10 +97,10 @@ function barycentric_interp(x, z, w, f)
 end
 
 
-@doc """
+"""
 Kopriva 2009, p. 75, (3.38).
 """
-function barycentric_interp_matrix(x, z)
+function barycentric_interp_matrix(x::AbstractVector, z::AbstractVector)
   w = barycentric_weights(z)
   Nx = length(x)
   Nz = length(z)
@@ -101,12 +123,12 @@ end
 
 
 
-@doc """
+"""
 Enhance accuracy of derivativ matrix D by computing diagonal elements such
 that derivative of a constant function vanishes.
 Kopriva 2009, p. 55, (2.43).
 """
-function negative_sum_trick!(D)
+function negative_sum_trick!(D::AbstractMatrix)
   N,M = size(D)
   @assert N == M
   for i = 1:N
@@ -123,11 +145,11 @@ end
 
 
 
-@doc """
+"""
 Kopriva 2009, p. 74, (3.34).
 x ... Interpolation nodes
 """
-function first_derivative_matrix(x)
+function first_derivative_matrix(x::AbstractVector)
   N = length(x)
   bw = barycentric_weights(x)
   D = zeros(N, N)
@@ -144,12 +166,12 @@ end # function
 
 
 
-@doc """
+"""
 Kopriva 2009, p. 82, (3.50).
 x ... Interpolation nodes
 bw ... Barycentric weights based on x
 """
-function second_derivative_matrix(x)
+function second_derivative_matrix(x::AbstractVector)
   N = length(x)
   bw = barycentric_weights(x)
   D = first_derivative_matrix(x)
@@ -168,7 +190,7 @@ end # function
 
 
 
-function nth_derivative_matrix(x)
+function nth_derivative_matrix(x::AbstractVector)
   N = lenght(x)
   bw = barycentric_weights(x)
   error("Implement nth derivative matrix!")
@@ -176,7 +198,7 @@ end
 
 
 
-function derivative_matrix(k, x)
+function derivative_matrix(k, x::AbstractVector)
   @assert k > 0 "Invalid dervative order `k = $k` requested!"
   if k == 1
     return first_derivative_matrix(x)
@@ -196,7 +218,7 @@ Establishes connection between nodal and modal coefficients for the Legendre bas
 the assumption that `u ≈ u_nodal(z_i) = u_modal(z_i)`.
 `z` is the nodal point set at which the approximations should agree on.
 """
-function vandermonde_matrix_legendre(z)
+function vandermonde_matrix_legendre(z::AbstractVector)
   Npts = length(z)
   V = [ legendre(z[i], j-1) / sqrt(2 / (2*(j-1)+1)) for i=1:Npts, j=1:Npts ]
   return V
@@ -212,7 +234,7 @@ the assumption that `u ≈ u_nodal(z_i) = u_modal(z_i)`.
 evaluated at node `z[i]`, where the `N+1` basis functions are labeled as `j = 0,...,N`.
 `z` is the nodal point set at which the approximations should agree on.
 """
-function vandermonde_matrix(basis_fn::Function, z)
+function vandermonde_matrix(basis_fn::Function, z::AbstractVector)
   Npts = length(z)
   V = [ basis_fn(z[i], j-1) for i=1:Npts, j=1:Npts ]
   return V
@@ -222,13 +244,13 @@ end
 @deprecate vandermonde_matrix(z) vandermonde_matrix_legendre(z)
 
 
-@doc """
+"""
 Spectral filtering matrix, motivated by an artificial viscosity regularization.
 Cf. Hesthaven, Warburton 2007, section 5.3,
 PhD thesis Glaubitz 2020, section 8.5.
 see hesthaven 2008, section 5.3
 """
-function filter_matrix(z, α = 36.0, Nc = 0, s = 16.0)
+function filter_matrix(z::AbstractVector, α = 36.0, Nc = 0, s = 16.0)
   Npts = length(z)
   N = Npts - 1
   σ(n) = exp(- α * ((n - Nc) / (N - Nc))^s)
diff --git a/src/legendre.jl b/src/legendre.jl
new file mode 100644
index 0000000000000000000000000000000000000000..fab3a9d7853611e17cc996fac6e4e615645cc4cf
--- /dev/null
+++ b/src/legendre.jl
@@ -0,0 +1,102 @@
+module Legendre
+
+
+using dg1d
+using LinearAlgebra
+using Jacobi
+
+
+# Normalized so that bernstein_polynomial(-1,0,N) == bernstein_polynomial(1,N,N) == 1
+function legendre_polynomial(x::Real, n::Integer)
+  Jacobi.legendre(x::Real, n)
+end
+
+
+function derivative_legendre_polynomial(x::Real, n::Integer)
+  Jacobi.dlegendre(x, n)
+end
+
+
+function vandermonde_matrix(z::AbstractVector)
+  Npts = length(z)
+  V = [ legendre_polynomial(z[i], j-1, Npts-1) for i=1:Npts, j=1:Npts ]
+  return V
+end
+
+
+"""
+  mass_matrix(Npts::Integer)
+
+Compute M_ij = (Pi, Pj) exactly where `Pi, Pj` are (i-1)th and (j-1)th Legendre polynomials.
+"""
+function mass_matrix(Npts::Integer)
+  N = Npts-1
+  M = [ i==j ? 2/(2*(i-1)+1) : 0.0 for i in 1:Npts, j in 1:Npts ]
+end
+
+
+"""
+  stiffness_matrix(Npts::Integer)
+
+Compute S_ij = (Pi, dPj/dx) exactly where `Pi, Pj` are (i-1)th and (j-1)th Bernstein polynomials.
+"""
+function stiffness_matrix(Npts::Integer)
+  N = Npts-1
+  S = zeros(Float64, Npts, Npts)
+  xs, ws = dg1d.LGL.rule(Npts+1)
+  for i in 1:Npts, j in 1:Npts
+    sij = dg1d.LGL.integrate(xs, ws) do xx
+      Pi = legendre_polynomial(xx, i-1)
+      ∂Pj = derivative_legendre_polynomial(xx, j-1)
+      return Pi*∂Pj
+    end
+    S[i,j] = sij
+  end
+  return S
+end
+
+
+"""
+  interpolate(x, fn::AbstractVector)
+
+Interpolate a Legendre polynomial with modal coefficients `fn` onto `x`.
+`x` is assumed to lie within `[-1,1]`.
+
+Use `legendre_projection` to compute modal coefficients from an analytical function.
+"""
+function interpolate(x::Real, fn::AbstractVector)
+  Npts = length(fn)
+  return sum(enumerate(fn)) do (i,fi)
+    Lim1 = legendre_polynomial(x, i-1)
+    fi*Lim1
+  end
+end
+
+
+"""
+  legendre_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,
+then use `Npts=N+1` to obtain that exact result.
+
+See `LGL.rule` to obtain the corresponding modal coefficients.
+"""
+function projection(f::Function, Npts::Integer)
+  N = Npts-1
+  xs, ws = dg1d.LGL.rule(Npts+1)
+  fP = j -> dg1d.LGL.integrate(xs, ws) do xi
+    fi = f(xi)
+    Pjm1 = legendre_polynomial(xi, j-1)
+    fi*Pjm1
+  end
+  PP = i -> dg1d.LGL.integrate(xs, ws) do xi
+    Pim1 = legendre_polynomial(xi, i-1)
+    Pim1*Pim1
+  end
+  coeffs = [ fP(i)/PP(i) for i in 1:Npts ]
+  return coeffs
+end
+
+
+end
diff --git a/src/lgl.jl b/src/lgl.jl
index 6a7c97e859347df251d63d181ce136bb110301eb..68a1a5560b863bc7287ca246dad5da6a7f0c4dca 100644
--- a/src/lgl.jl
+++ b/src/lgl.jl
@@ -85,24 +85,29 @@ end
 
 
 """
-Inverse of exact full mass matrix, computed using Sherman-Morrison formula,
-cf. PhD thesis Marcus Bugner, Notes spectral elements Bernd Bruegmann,
-Kopriva (2013), Teukolsky (2015).
-n_pts ... number of qudrature nodes = number of lagrange polynomials
-"""
-function inverse_mass_matrix_exact(n_pts::Integer)
-  z, w = rule(n_pts)
-  N = n_pts - 1
-  invM = zeros(n_pts, n_pts)
+  inverse_mass_matrix_exact(Npts::Integer)
+
+Compute the exact inverse of M_ij = (li, lj) where `li, lj` are Lagrange polynomials
+defined on the LGL quadrature nodes.
+Cf. PhD thesis Marcus Bugner, Kopriva (2013), Teukolsky (2015).
+"""
+function inverse_mass_matrix_exact(Npts::Integer)
+  z, w = rule(Npts)
+  N = Npts - 1
+  invM = zeros(Npts, Npts)
   P = legendre.(z, N)
   invM = diagm(1.0 ./ w)
-  invM += 0.5 * n_pts * P * transpose(P)
+  invM += 0.5 * Npts * P * transpose(P)
   return invM
 end
 
 
 """
-n_pts ... number of qudrature nodes = number of lagrange polynomials
+  inverse_mass_matrix_lumped(Npts::Integer)
+
+Compute an approximate inverse of M_ij = (li, lj) where `li, lj` are Lagrange polynomials
+defined on the LGL quadrature nodes.
+Cf. Teukolsky (2015), Kopriva and Gassner (2011).
 """
 function inverse_mass_matrix_lumped(Npts::Integer)
   _, w = rule(Npts)
@@ -112,8 +117,10 @@ end
 
 
 """
-M_ij = (li, lj)
-n_pts ... number of qudrature nodes = number of lagrange polynomials
+  mass_matrix_exact(Npts::Integer)
+
+Compute M_ij = (li, lj) exactly where `li, lj` are Lagrange polynomials
+defined on the LGL quadrature nodes.
 """
 function mass_matrix_exact(Npts::Integer)
   z, w = rule(Npts)
@@ -136,9 +143,10 @@ end
 
 
 """
-M_ij = (li, lj)
-Approximated mass matrix due to under integration of LGL quadrature.
-n_pts ... number of qudrature nodes = number of lagrange polynomials
+  mass_matrix_lumped(Npts::Integer)
+
+Compute M_ij = (li, lj) approximately where `li, lj` are Lagrange polynomials
+defined on the LGL quadrature nodes.
 """
 function mass_matrix_lumped(Npts::Integer)
   z, w = rule(Npts)
@@ -148,10 +156,10 @@ end
 
 
 """
-S_ij = (li, dlj/dx)
-Exactly representable, because product is of order 2N - 1 which can be integrated
-exactly using LGL quadrature.
-n_pts ... number of quadrature nodes = number of lagrange polynomials
+  stiffness_matrix(Npts::Integer)
+
+Compute S_ij = (li, dlj/dx) exactly where `li, lj` are Lagrange polynomials
+defined on the LGL quadrature nodes.
 """
 function stiffness_matrix(Npts::Integer)
   z, w = rule(Npts)
diff --git a/src/main.jl b/src/main.jl
index ace44703121ae03575ab4b7aea2dbfb4cfb92812..8ae0bbfca6d366de284c83a3e147f920c3268e88 100644
--- a/src/main.jl
+++ b/src/main.jl
@@ -203,9 +203,8 @@ end
 
 function setup_env(project_name, prms)
   mprms = prms["Mesh"]
-  mesh = Mesh(; N=mprms["n"], K=mprms["k"], range=mprms["range"],
-              basis=mprms["basis"], periodic=mprms["periodic"],
-              scheme=mprms["scheme"],
+  mesh = Mesh(; N=mprms["n"], K=mprms["k"], range=mprms["range"], periodic=mprms["periodic"],
+              scheme=mprms["scheme"], basis=mprms["basis"], dg_kind=mprms["dg_kind"],
               filter_parameters=Tuple(mprms["filter_parameters"]))
   callbacks = CallbackSet()
   callbacks_substeps = CallbackSet()
diff --git a/src/mesh.jl b/src/mesh.jl
index 75ed954155d9ccde2257617119d4c157466e19d9..0202180bd55a6462e0cdbe9e1a8f0403f21a0ab5 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -30,7 +30,8 @@ const Mesh2d{Element} = Mesh{Element,2,4}
 #######################################################################
 
 
-function Mesh1d(; N=5, K=10, range=[-1.0,1.0], basis="lgl_lumped", periodic=(true,), scheme="DG",
+function Mesh1d(; N=5, K=10, range=[-1.0,1.0], periodic=(true,),
+                  scheme="DG", basis="lgl_lumped", dg_kind="nodal_lagrange",
                   filter_parameters=parameters(:Mesh)["filter_parameters"])
 
   @toggled_assert N > 0
@@ -42,10 +43,9 @@ function Mesh1d(; N=5, K=10, range=[-1.0,1.0], basis="lgl_lumped", periodic=(tru
   tree = Tree1d(K,periodic=tuple(periodic...))
   boxes = place_boxes(tree, Float64.(range))
   element = if scheme == "DG"
-    # TODO Remove this after updating parameter names in SpectralElement
-    quadrature_method = Dict("lgl" => :LGL, "glgl" => :GLGL, "lgl_lumped" => :LGL_lumped)[basis]
+    # TODO Remove this after updating parameter names in DGElement
     ηc, sh, α = filter_parameters
-    SpectralElement(N, Symbol(quadrature_method); ηc, sh, α)
+    DGElement(N, Symbol(basis), Symbol(dg_kind); ηc, sh, α)
   elseif scheme == "FV"
     FVElement()
   else
@@ -98,9 +98,8 @@ function Mesh1d(; N=5, K=10, range=[-1.0,1.0], basis="lgl_lumped", periodic=(tru
 end
 
 
-function Mesh2d(; N=[5,5], K=[4,4], range=[-1.0,1.0, -1.0,1.0],
-                  basis="lgl_lumped", periodic=(true,true),
-                  scheme="DG",
+function Mesh2d(; N=[5,5], K=[4,4], range=[-1.0,1.0, -1.0,1.0], periodic=(true,true),
+                  scheme="DG", basis="lgl_lumped", dg_kind="nodal_lagrange",
                   filter_parameters=parameters(:Mesh)["filter_parameters"])
 
   @toggled_assert length(N) == 2
@@ -121,11 +120,10 @@ function Mesh2d(; N=[5,5], K=[4,4], range=[-1.0,1.0, -1.0,1.0],
   tree = Tree2d(Kx,Ky,periodic=tuple(periodic...))
   boxes = place_boxes(tree, Float64.(xrange), Float64.(yrange))
   if scheme == "DG"
-    # TODO Remove this after updating parameter names in SpectralElement
-    quadrature_method = Dict("lgl" => :LGL, "glgl" => :GLGL, "lgl_lumped" => :LGL_lumped)[basis]
+    # TODO Remove this after updating parameter names in DGElement
     ηc, sh, α = filter_parameters
-    element_x = SpectralElement(Nx, Symbol(quadrature_method); ηc, sh, α)
-    element_y = SpectralElement(Ny, Symbol(quadrature_method); ηc, sh, α)
+    element_x = DGElement(Nx, Symbol(basis), Symbol(dg_kind); ηc, sh, α)
+    element_y = DGElement(Ny, Symbol(basis), Symbol(dg_kind); ηc, sh, α)
   elseif scheme == "FV"
     TODO()
     element_x = FVElement()
@@ -235,12 +233,13 @@ function Mesh(; N=parameters(:Mesh)["n"],
                 K=parameters(:Mesh)["k"],
                 range=parameters(:Mesh)["range"],
                 basis=parameters(:Mesh)["basis"],
+                dg_kind=parameters(:Mesh)["dg_kind"],
                 periodic=parameters(:Mesh)["periodic"],
                 scheme=parameters(:Mesh)["scheme"],
                 filter_parameters=parameters(:Mesh)["filter_parameters"])
   dims = Int(length(range)/2)
   MType = dims == 1 ? Mesh1d : Mesh2d
-  return MType(; N, K, range, basis, periodic, scheme, filter_parameters)
+  return MType(; N, K, range, basis, dg_kind, periodic, scheme, filter_parameters)
 end
 
 
@@ -268,7 +267,7 @@ end
 extends(m::Mesh,cell::Cell) = m.boxes[cell.index].extends
 
 
-function min_grid_spacing(m::Mesh{SpectralElement})
+function min_grid_spacing(m::Mesh{DGElement})
   @unpack z = m.element
   ws = widths(m.boxes[1])
   minw = minimum(ws)
@@ -395,7 +394,7 @@ function integrate_cell(u, k, mesh::Mesh2d)
 end
 
 
-function integrate(fn::Function, mesh::Mesh1d{SpectralElement})
+function integrate(fn::Function, mesh::Mesh1d{DGElement})
   @unpack invdetJ = get_static_variables(mesh)
   Npts, K = layout(mesh)
   mat_invdetJ = vreshape(invdetJ, (Npts,K))
diff --git a/src/numutils.jl b/src/numutils.jl
index c6e17fc7ff951d69dc2cc204947ecb53c39acfc5..7061c14229ab161c56e158f8fd49f6c7675b18dd 100644
--- a/src/numutils.jl
+++ b/src/numutils.jl
@@ -37,10 +37,12 @@ end
 
 """
     purge_zeros!(M, thrs=1e-13)
+    purge_zeros(M, thrs=1e-13)
 
 Scan `M` for elements smaller than `thrs` and set them to 0.
 """
-purge_zeros!(M, thrs=1e-13) = M[abs.(M).<thrs] .= 0.0
+purge_zeros!(M, thrs=1e-13) = (M[abs.(M).<thrs] .= 0.0; M)
+purge_zeros(M, thrs=1e-13) = purge_zeros!(copy(M), thrs)
 
 
 """
@@ -121,3 +123,43 @@ end
 Return absolute maximum of `vs...`.
 """
 @inline absmax(vs...) = maximum(abs, vs)
+
+
+"""
+  interpolate_piecewise_constant(x::Real, xs::AbstractVector, fn::AbstractVector)
+
+Interpolate function values `fn` sampled on nodes `xs` in a piecewise constant manner
+onto `x`. This function assumes `xs` to be ordered increasingly, and
+`x` is assumed to lie within `extrema(xs)`.
+"""
+function interpolate_piecewise_constant(x::Real, xs::AbstractVector, fn::AbstractVector)
+  Npts = length(fn)
+  N = Npts-1
+  for i in 2:Npts
+    xim1, xi = xs[i-1], xs[i]
+    if x <= xim1+(xi-xim1)/2
+      return fn[i-1]
+    end
+  end
+  return fn[end]
+end
+
+
+"""
+  interpolate_piecewise_linear(x::Real, xs::AbstractVector, fn::AbstractVector)
+
+Interpolate function values `fn` sampled on nodes `xs` in a piecewise linear manner
+onto `x`. This function assumes `xs` to be ordered increasingly, and
+`x` is assumed to lie within `extrema(xs)`.
+"""
+function interpolate_piecewise_linear(x, xs::AbstractVector, fn::AbstractVector)
+  Npts = length(fn)
+  N = Npts-1
+  for i in 1:Npts-1
+    xi, xip1 = xs[i], xs[i+1]
+    if x <= xip1
+      return fn[i] + (fn[i+1]-fn[i])/(xip1-xi)*(x-xi)
+    end
+  end
+  return fn[end]
+end
diff --git a/src/parameters.jl b/src/parameters.jl
index 7539c9e8f7e83e91453e7540a34a878f5f340758..c13ec135e015face3dc637511955a7165168dd59 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -486,6 +486,7 @@ end
   @check length(n) in [ 1, 2 ]
   @check all(n .>= 0)
 
+  # TODO Rename basis to quadrature
   """
   Nodal basis type used for the polynomial approximation. Available options
   - `lgl_lumped` ... mass-lumped Legendre-Gauss-Lobatto grid
@@ -503,6 +504,19 @@ end
   scheme = "DG"
   @check scheme in [ "DG", "FV" ]
 
+  # TODO Rename this to basis
+  """
+  DG scheme used for the evolution. Available options
+  - `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.
+  - `nodal_bernstein` ... (Quasi-)Nodal approximation using Bernstein polynomials.
+    The coefficients correspond to the nodal values sampled on an equidistant grid.
+    See Beisiegel, Behrens, Environ Earth Sci (2015) 74:7275–7284.
+  """
+  dg_kind = "nodal_lagrange"
+  @check scheme in [ "nodal_lagrange", "nodal_bernstein" ]
+
   """
   Periodicity in cartesian directions `x,y`
   - 1d: [ `x_dir` ]
diff --git a/src/spectralelement.jl b/src/spectralelement.jl
deleted file mode 100644
index 81b56b75597663f8e3b9e77288cba8baa86d59b3..0000000000000000000000000000000000000000
--- a/src/spectralelement.jl
+++ /dev/null
@@ -1,121 +0,0 @@
-Base.@kwdef struct SpectralElement
-  N::Int64                # polynomial order
-  Npts::Int64             # number of quadrature points = N + 1
-  z::Vector{Float64}      # quadrature points
-  w::Vector{Float64}      # quadrature weights
-  D::Matrix{Float64}      # derivative matrix
-  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
-  Ml_lhs::Vector{Float64} # = M^-1 l(-1)
-  Ml_rhs::Vector{Float64} # = M^-1 l(+1)
-  V::Matrix{Float64}      # = P_(j-1)(z_i), Vandermonde matrix
-  invV::Matrix{Float64}   # inv(V), inverse Vandermonde matrix
-  quadr::Symbol           # :LGL, :LGL_mass_lumped, :GLGL
-  F::Matrix{Float64}      # spectral filtering matrix
-
-  function SpectralElement(N::Integer, quadr::Symbol;
-      # filter parameters
-      ηc::Float64=0.5, sh::Float64=16.0, α::Float64=36.0)
-    Npts = N + 1
-
-    if quadr === :LGL
-      z, w        = LGL.rule(Npts)
-      v           = -1.0 # unused here
-      purge_zeros!(z) # ensures that the centered node for even N basis sits exactly at x = 0.
-      D           = first_derivative_matrix(z)
-      S           = LGL.stiffness_matrix(Npts)
-      M           = LGL.mass_matrix_exact(Npts)
-      invM        = LGL.inverse_mass_matrix_exact(Npts)
-    elseif quadr === :LGL_lumped
-      z, w        = LGL.rule(Npts)
-      v           = -1.0 # unused here
-      purge_zeros!(z) # ensures that the centered node for even N basis sits exactly at x = 0.
-      D           = first_derivative_matrix(z)
-      S           = LGL.stiffness_matrix(Npts)
-      M           = LGL.mass_matrix_lumped(Npts)
-      invM        = LGL.inverse_mass_matrix_lumped(Npts)
-    elseif quadr === :GLGL
-      z, w        = GLGL.rule(Npts)
-      purge_zeros!(z) # ensures that the centered node for even N basis sits exactly at x = 0.
-      D           = first_derivative_matrix(z)
-      D2          = second_derivative_matrix(z)
-      M           = GLGL.mass_matrix(Npts)
-      S           = GLGL.stiffness_matrix(Npts)
-      invM        = inv(M)
-    else
-      error("Invalid quadrature rule specified! Use one of :LGL_lumped, :LGL, :GLGL.")
-    end
-
-    MDM     = invM * transpose(S)
-    purge_zeros!(MDM)
-    V       = vandermonde_matrix(z) # TODO needs to be replaced with vandermonde_matrix(z, basis_fn)
-    invV    = inv(V)
-    Λ = diagm([ spectral_filter((i-1)/N, (ηc,sh,α)) for i = 1:Npts ])
-    F = V * Λ * invV
-
-    l_lhs   = zeros(N+1)
-    l_rhs   = zeros(N+1)
-    l_lhs[1]    = 1.0
-    l_rhs[end]  = 1.0
-    Ml_lhs = invM * l_lhs
-    Ml_rhs = invM * l_rhs
-
-    return new(N, Npts, z, w, D, S, M, invM, MDM, Ml_lhs, Ml_rhs,
-               V, invV, quadr, F)
-  end
-end
-
-
-# (5.16) from Hesthaven, Warburton 2007
-function spectral_filter(η::Real, prms)
-  ηc, sh, α = prms
-  return if 0 ≤ η ≤ ηc
-    1
-  else
-    exp(-α*((η-ηc)/(1-ηc))^(2*sh))
-  end
-end
-
-
-@inline function inner_product(u1::AbstractVector, u2::AbstractVector, el::SpectralElement)
-  @unpack Npts, M = el
-  result = 0.0
-  @turbo for i in 1:Npts, j in 1:Npts
-    result += u1[i] * M[i,j] * u2[j]
-  end
-  return result
-end
-
-
-@inline function integrate(u::AbstractVector, el::SpectralElement)
-  @unpack quadr = el
-  @toggled_assert quadr in (:LGL, :LGL_lumped, :GLGL)
-  if quadr === :LGL || quadr === :LGL_lumped
-    return LGL.integrate(u, el.z, el.w)
-  else # quadr === :GLGL
-    return GLGL.integrate(u, el.z, el.w)
-  end
-end
-
-
-@inline function integrate(fn::Function, idxs::AbstractVector{<:Integer}, el::SpectralElement)
-  @unpack quadr = el
-  @toggled_assert quadr in (:LGL, :LGL_lumped, :GLGL)
-  if quadr === :LGL || quadr === :LGL_lumped
-    return LGL.integrate(fn, idxs, el.z, el.w)
-  else # quadr === :GLGL
-    return GLGL.integrate(fn, idxs, el.z, el.w)
-  end
-end
-
-
-@inline function differentiate(u::AbstractArray, el::SpectralElement)
-  return el.D * u
-end
-
-
-@inline function differentiate!(du::AbstractArray, u::AbstractArray, el::SpectralElement)
-  mul!(du, el.D, u)
-end
diff --git a/src/tensorbasis.jl b/src/tensorbasis.jl
index fabb58a84e2c5cb18e93961ed13fcea6e42d46b8..1710857ce29b9cde5981a8cfde0d3e527bf23871 100644
--- a/src/tensorbasis.jl
+++ b/src/tensorbasis.jl
@@ -7,9 +7,9 @@ abstract type AbstractBasis end
 
 
 struct TensorBasis{N_Dim}
-  basis::NTuple{N_Dim,SpectralElement}
+  basis::NTuple{N_Dim,DGElement}
 
-  function TensorBasis{N}(bases::Vararg{SpectralElement,N}) where N
+  function TensorBasis{N}(bases::Vararg{DGElement,N}) where N
     @assert 1 <= N <= 3 "Only support 1D, 2D, 3D rectangular basis!"
     return new{N}(bases)
   end
diff --git a/test/UnitTests/src/test_spectral.jl b/test/UnitTests/src/test_dgelement.jl
similarity index 68%
rename from test/UnitTests/src/test_spectral.jl
rename to test/UnitTests/src/test_dgelement.jl
index cdb3f6f47d17015d4325fe53ed2179e7b4dcaa3a..cf43f55b56bf7f5969a20a588c5402392d77fe5d 100644
--- a/test/UnitTests/src/test_spectral.jl
+++ b/test/UnitTests/src/test_dgelement.jl
@@ -1,7 +1,7 @@
 @testset "Spectral Element" begin
 
 
-function evaluate_derivatives(el::SpectralElement, f, df)
+function evaluate_derivatives(el::DGElement, f, df)
   f_ana = f.(el.z)
   df_ana = df.(el.z)
   df_num = el.D * f_ana
@@ -16,15 +16,17 @@ end
   df(x) = 3.0 * x^2 - 10.0 * x
   # ddf(x) = 6.0 * x - 10.0
 
-  el = SpectralElement(N, :LGL)
+  dg_kind = :nodal_lagrange
+
+  el = DGElement(N, :lgl, dg_kind)
   df_ana, df_num = evaluate_derivatives(el, f, df)
   @test isapprox(df_ana, df_num)
 
-  el = SpectralElement(N, :LGL_lumped)
+  el = DGElement(N, :lgl_lumped, dg_kind)
   df_ana, df_num = evaluate_derivatives(el, f, df)
   @test isapprox(df_ana, df_num)
 
-  glgl = SpectralElement(N, :LGL)
+  glgl = DGElement(N, :lgl, dg_kind)
   df_ana, df_num = evaluate_derivatives(el, f, df)
   @test isapprox(df_ana, df_num)
 
@@ -40,16 +42,16 @@ end
   # LGL quadrature is only exact up to order 2*3-1=5, so use N=4 instead
   N = 4
   L = 2.0
-  quadr = :LGL
-  el = SpectralElement(N, quadr)
+  quadr = :lgl
+  el = DGElement(N, quadr, :nodal_lagrange)
   fx = f.(el.z)
   f2 = dg1d.inner_product(fx, fx, el)
   @test isapprox(f2, f2_ana)
 
-  # check that :GLGL quadrature is exact up to and including order 2N
+  # check that :glgl quadrature is exact up to and including order 2N
   N = 3
   L = 2.0
-  el = SpectralElement(N, :GLGL)
+  el = DGElement(N, :glgl, :nodal_lagrange)
   fx = f.(el.z)
   f2 = dg1d.inner_product(fx, fx, el)
   @test isapprox(f2, f2_ana)
@@ -195,4 +197,51 @@ end
 end
 
 
+@testset "verify some DG operator identities" begin
+
+  N = 10
+
+  for basis in (:lgl, :glgl)
+    el = dg1d.DGElement(N, basis, :nodal_lagrange)
+    # verify: M D = S
+    MD = dg1d.purge_zeros(el.M*el.D)
+    S = dg1d.purge_zeros(el.S)
+    @test all(MD.≈S)
+  end
+
+  el = dg1d.DGElement(N, :lgl, :nodal_bernstein)
+  T, invT = dg1d.Bernstein.bernstein_legendre_transformation(el.Npts)
+  M_legendre = dg1d.Legendre.mass_matrix(el.Npts)
+  TMT = T*M_legendre*transpose(T) |> dg1d.purge_zeros
+  @test all(el.M.≈TMT)
+  S_legendre = dg1d.Legendre.stiffness_matrix(el.Npts)
+  TST = T*S_legendre*transpose(T) |> dg1d.purge_zeros
+  @test all(el.S.≈TST)
+
+end
+
+
+@testset "Legendre expansions and Bernstein polynomials" begin
+
+  # approximate the Bernstein basis through Legendre polynomials
+  N = 10
+  Npts = N+1
+  Cs = zeros(Float64, Npts, Npts)
+  for i in 1:Npts, j in 1:Npts
+    # ith Bernstein polynomial, jth Legendre modal coefficient
+    Cs[i,:] .= dg1d.Legendre.projection(Npts) do xi
+      dg1d.Bernstein.bernstein_polynomial(xi, i-1, N)
+    end
+  end
+
+  xx = collect(range(-1.0,1.0,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
+
+end
+
+
 end
diff --git a/test/UnitTests/src/test_tensorbasis.jl b/test/UnitTests/src/test_tensorbasis.jl
index eca39a0833578732a7ef94113c26e10e7d5c7c94..2fef6c8795ac13ceb47040b14370cd8593e0ec06 100644
--- a/test/UnitTests/src/test_tensorbasis.jl
+++ b/test/UnitTests/src/test_tensorbasis.jl
@@ -1,7 +1,7 @@
 @testset "TensorBasis" begin
 
   nx, ny, nz = 7,8,9
-  elx, ely, elz = [ SpectralElement(Np-1, :LGL) for Np in (nx,ny,nz) ]
+  elx, ely, elz = [ DGElement(Np-1, :lgl, :nodal_lagrange) for Np in (nx,ny,nz) ]
   tb1 = TensorBasis1d(elx)
   tb2 = TensorBasis2d(elx, ely)
   tb3 = TensorBasis3d(elx, ely, elz)
diff --git a/test/UnitTests/src/tests.jl b/test/UnitTests/src/tests.jl
index fcdfd78f8df6242384f35f529cc2f67de29f35fe..a35df0d037f2e3a52cc4b22708d25ae64cdac588 100644
--- a/test/UnitTests/src/tests.jl
+++ b/test/UnitTests/src/tests.jl
@@ -10,7 +10,7 @@ include("utils.jl")
 include("test_utils.jl")
 include("test_objectcache.jl")
 include("test_with_signature.jl")
-include("test_spectral.jl")
+include("test_dgelement.jl")
 include("test_tensorbasis.jl")
 include("test_box.jl")
 include("test_tree.jl")