diff --git a/src/mesh.jl b/src/mesh.jl
index 75ed954155d9ccde2257617119d4c157466e19d9..e527633a1ba650243225281f817cc6b5503d5d79 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,),
+                  basis="lgl_lumped", dg_kind="nodal_legendre", scheme="DG",
                   filter_parameters=parameters(:Mesh)["filter_parameters"])
 
   @toggled_assert N > 0
@@ -45,7 +46,7 @@ function Mesh1d(; N=5, K=10, range=[-1.0,1.0], basis="lgl_lumped", periodic=(tru
     # TODO Remove this after updating parameter names in SpectralElement
     quadrature_method = Dict("lgl" => :LGL, "glgl" => :GLGL, "lgl_lumped" => :LGL_lumped)[basis]
     ηc, sh, α = filter_parameters
-    SpectralElement(N, Symbol(quadrature_method); ηc, sh, α)
+    SpectralElement(N, Symbol(quadrature_method), Symbol(dg_kind); ηc, sh, α)
   elseif scheme == "FV"
     FVElement()
   else
@@ -98,9 +99,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),
+                  basis="lgl_lumped", dg_kind="nodal_legendre", scheme="DG",
                   filter_parameters=parameters(:Mesh)["filter_parameters"])
 
   @toggled_assert length(N) == 2
@@ -124,8 +124,8 @@ function Mesh2d(; N=[5,5], K=[4,4], range=[-1.0,1.0, -1.0,1.0],
     # TODO Remove this after updating parameter names in SpectralElement
     quadrature_method = Dict("lgl" => :LGL, "glgl" => :GLGL, "lgl_lumped" => :LGL_lumped)[basis]
     ηc, sh, α = filter_parameters
-    element_x = SpectralElement(Nx, Symbol(quadrature_method); ηc, sh, α)
-    element_y = SpectralElement(Ny, Symbol(quadrature_method); ηc, sh, α)
+    element_x = SpectralElement(Nx, Symbol(quadrature_method), Symbol(dg_kind); ηc, sh, α)
+    element_y = SpectralElement(Ny, Symbol(quadrature_method), Symbol(dg_kind); ηc, sh, α)
   elseif scheme == "FV"
     TODO()
     element_x = FVElement()
@@ -235,12 +235,13 @@ function Mesh(; N=parameters(:Mesh)["n"],
                 K=parameters(:Mesh)["k"],
                 range=parameters(:Mesh)["range"],
                 basis=parameters(:Mesh)["basis"],
+                dg_kind=parameters(:Mesh)["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
 
 
@@ -269,6 +270,7 @@ extends(m::Mesh,cell::Cell) = m.boxes[cell.index].extends
 
 
 function min_grid_spacing(m::Mesh{SpectralElement})
+  # TODO Use nodal_points(m)
   @unpack z = m.element
   ws = widths(m.boxes[1])
   minw = minimum(ws)
diff --git a/src/spectralelement.jl b/src/spectralelement.jl
index 30000da3d97e5a8e9f8056f2fdb74192a4cbf83b..3a2598566046029c4d5f8c68bc88a0b378eaa8f6 100644
--- a/src/spectralelement.jl
+++ b/src/spectralelement.jl
@@ -36,35 +36,33 @@ struct SpectralElement
       if quadr === :LGL
         qaudr_z, qaudr_w = LGL.rule(Npts); qaudr_v = -1.0 # unused here
         purge_zeros!(qaudr_z)
-        z    = quadr_z
-        D    = first_derivative_matrix(z)
         S    = LGL.stiffness_matrix(Npts)
         invM = LGL.inverse_mass_matrix_exact(Npts)
       elseif quadr === :LGL_lumped
         quadr_z, quadr_w = LGL.rule(Npts); quadr_v = -1.0 # unused here
         purge_zeros!(quadr_z)
-        z    = quadr_z
-        D    = first_derivative_matrix(z)
         S    = LGL.stiffness_matrix(Npts)
         invM = LGL.inverse_mass_matrix_lumped(Npts)
       elseif quadr === :GLGL
         quadr_z, quadr_w, quadr_v = GLGL.rule(Npts)
         purge_zeros!(quadr_z)
-        z    = quadr_z
-        D    = first_derivative_matrix(z)
-        D2   = second_derivative_matrix(z)
-        M    = GLGL.mass_matrix(z, w, v, D)
-        S    = GLGL.stiffness_matrix(z, w, v, D, D2)
+        quadr_D  = first_derivative_matrix(quadr_z)
+        quadr_D2 = second_derivative_matrix(quadr_z)
+        M    = GLGL.mass_matrix(quadr_z, quadr_w, quadr_v, quadr_D)
+        S    = GLGL.stiffness_matrix(quadr_z, quadr_w, quadr_v, quadr_D, quadr_D2)
         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)
+      quadr_D = D
       l_lhs = zeros(Float64, N+1)
       l_rhs = zeros(Float64, N+1)
       l_lhs[1]   = 1.0
       l_rhs[end] = 1.0
-      quadr_D = D
       V = vandermonde_matrix_legendre(z)
+    # In preparation for Bernstein basis
     # elseif kind === :nodal_bernstein
     #   if quadr === :LGL || quadr === :LGL_lumped
     #     qaudr_z, qaudr_w = LGL.rule(Npts)
@@ -82,9 +80,11 @@ struct SpectralElement
     #   S = Bernstein.stiffness_matrix(Npts)
     #   invM = inv(M)
     #   D = Bernstein.first_derivative_matrix(Npts)
-    #   D2 = Bernstein.second_derivative_matrix(Npts)
+    #   TODO("quadr_D")
+    #   TODO("quadr_w")
     #   V  = Bernstein.vandermonde_matrix(z)
     #   I  = stack(Bernstein.bernstein_polynomial.(quadr_z, i-1, N) for i in 1:Npts)
+    #   quadr_w = transpose(I) * quadr_w
     #   quadr_D = D * I
     #   l_lhs   = zeros(Float64, N+1)
     #   l_rhs   = zeros(Float64, N+1)