diff --git a/src/HRSC/bernstein.jl b/src/HRSC/bernstein.jl
index c5873cb8d94980f6078db64bbfff29d6c8fa5f25..22c4cfca3d665ef2d988de281acd65dcd2da155e 100644
--- a/src/HRSC/bernstein.jl
+++ b/src/HRSC/bernstein.jl
@@ -55,7 +55,7 @@ end
 @doc """
     reconstruct_cell!(u!, w, P, V; kwargs...)
 
-Compute inplace a convex combination of `u!` and its Bernstein reconstruction weighted by `w`.  
+Compute inplace a convex combination of `u!` and its Bernstein reconstruction weighted by `w`.
 `P, V` are the project and transfer matrices that transfer data from the sample
 grid to the bernstein grid and back, respectively.
 `m` and `M` keywords allow to enforce lower and upper bound in reconstruction.
diff --git a/src/glgl.jl b/src/glgl.jl
index 980b8e927c1335fa4312358a000efed04441618b..47d47cd5e93664ea5c664d66877502772089ecc7 100644
--- a/src/glgl.jl
+++ b/src/glgl.jl
@@ -14,7 +14,7 @@ x ... quadrature nodes, vector
 w ... quadrature weights, vector
 v ... derivative weight, number
 """
-function rule(n_pts)
+function rule(n_pts::Integer)
 
   N = n_pts - 2
   gb_x, gb_w = [], []
@@ -43,14 +43,13 @@ end # function
 
 
 """
-Integrate with Legendre weight w(x) = 1.
-x ... Vector of GLGL nodes.
-w ... Vector of GLGL weights.
-v ... Number, GLGL boundary weight.
-D ... Derivative matrix based on x.
-f1, f2 ... Vector of function values evaluated on grid x.
+  integrate(f, x, w, v, D)
+
+Integrate `f` sampled on GLGL quadrature nodes `x` with quadrature weights `w`,
+boundary weight `v` and derivative matrix `D` based on `x`.
+Use `rule` to obtain `x, w, v` and use `first_derivative_matrix` to obtain `D`.
 """
-function integrate(w, v, D, f)
+function integrate(f::AbstractVector, x::AbstractVector, w::AbstractVector, v::Real, D::AbstractVector)
 
   N = length(w)
   @toggled_assert (N,N) == size(D)
@@ -62,7 +61,43 @@ function integrate(w, v, D, f)
 
   return sum
 end
-function integrate(w, v, D, fn::Function, idxs::AbstractVector{<:Integer})
+
+
+"""
+  integrate(f1, f2, w, v, D)
+
+Similar to `integrate(f, w)`, but computes the inner product of `f1*f2` where both
+`f1, f2` must be sampled on the GLGL quadrature nodes `x` associated with the quadrature
+weights `w`, boundary weight `v` and derivative matrix `D` based on `x`.
+Use `rule` to obtain `x, w, v` and use `first_derivative_matrix` to obtain `D`.
+"""
+function integrate(f1::AbstractVector, f2::AbstractVector, x::AbstractVector, w::AbstractVector,
+    v::Real, D::AbstractMatrix)
+
+  N = length(w)
+  @toggled_assert (N,N) == size(D)
+  @toggled_assert N == length(f1)
+  @toggled_assert N == length(f2)
+
+  sum = dot(f1.*f2, w)
+  vD1, vDN = view(D,1,:), view(D,N,:)
+  sum += v * (f1[1] * dot(vD1, f2) + f2[1] * dot(vD1, f1))
+  sum -= v * (f1[end] * dot(vDN, f2) + f2[end] * dot(vDN, f1))
+
+  return sum
+end
+
+
+"""
+  integrate(f::Function, idxs::AbstractVector{<:Integer}, x, w)
+
+Integrate `f` sampled on GLGL quadrature nodes `x` with quadrature weights `w`,
+boundary weight `v` and derivative matrix `D` based on `x`.
+`f(idxs[i])` shall return the value of the integrand evaluated at `x[i]`.
+Use `rule` to obtain `x, w, v` and use `first_derivative_matrix` to obtain `D`.
+"""
+function integrate(fn::Function, idxs::AbstractVector{<:Integer}, x::AbstractVector, w::AbstractVector,
+    v::Real, D::AbstractMatrix)
 
   N = length(w)
   @toggled_assert (N,N) == size(D)
@@ -78,18 +113,25 @@ function integrate(w, v, D, fn::Function, idxs::AbstractVector{<:Integer})
 end
 
 
-function integrate(w, v, D, f1, f2)
+"""
+  integrate(f::Function, x, w, v, D)
 
-  N = length(w)
-  @toggled_assert (N,N) == size(D)
-  @toggled_assert N == length(f1)
-  @toggled_assert N == length(f2)
-
-  sum = dot(f1.*f2, w)
-  vD1, vDN = view(D,1,:), view(D,N,:)
-  sum += v * (f1[1] * dot(vD1, f2) + f2[1] * dot(vD1, f1))
-  sum -= v * (f1[end] * dot(vDN, f2) + f2[end] * dot(vDN, f1))
+Integrate `f` sampled on GLGL quadrature nodes `x` with quadrature weights `w`,
+boundary weight `v` and derivative matrix `D` based on `x`.
+`f(x[i])` shall return the value of the integrand evaluated at `x[i]`.
+Use `rule` to obtain `x, w, v` and use `first_derivative_matrix` to obtain `D`.
 
+This is similar to `integrate(f, idxs, x, w, v, D)`, but is better suited if the integrand is
+given in closed form.
+"""
+function integrate(f::Function, x::AbstractVector, w::AbstractVector, v::Real, D::AbstractMatrix)
+  @toggled_assert length(w) == length(x) == size(D, 1) == size(D, 2)
+  sum = 0.0
+  N = length(w)
+  for i in 1:N
+    fi = f(x[i])
+    sum += (w[i] + v * (D[1,i] - D[N,i])) * fi
+  end
   return sum
 end
 
@@ -101,12 +143,12 @@ w ... GLGL weights
 v ... GLGL derivative weight
 D ... Derivative matrix based on x.
 """
-function mass_matrix(x, w, v, D)
+function mass_matrix(x::AbstractVector, w::AbstractVector, v::Real, D::AbstractMatrix)
 
   N = length(x)
   @assert N == length(w)
   @assert (N,N) == size(D)
-  
+
   M = diagm(w)
   M[1,:] += v * D[1,:]
   M[:,1] += v * D[1,:]
@@ -128,13 +170,14 @@ v ... GLGL derivative weight
 D ... First derivative matrix based on x.
 D2 ... Second derivative matrix based on x.
 """
-function stiffness_matrix(x, w, v, D, D2)
+function stiffness_matrix(x::AbstractVector, w::AbstractVector, v::Real,
+    D::AbstractMatrix, D2::AbstractMatrix)
 
   N = length(x)
   @assert N == length(w)
   @assert (N,N) == size(D)
   @assert (N,N) == size(D2)
-  
+
   S = w .* D
   S[1,:]    += v * D2[1,:]
   S[end,:]  -= v * D2[end,:]
@@ -145,7 +188,7 @@ function stiffness_matrix(x, w, v, D, D2)
 end
 
 
-# TODO Add versions of mass_matrix, stiffness_matrix that take 
+# TODO Add versions of mass_matrix, stiffness_matrix that take
 # a single argument n_pts.
 
 
@@ -159,14 +202,14 @@ w ... GLGL weights
 v ... GLGL derivative weight
 D ... First derivative matrix based on x.
 """
-function linear_dependency_matrix(x, w, v, D)
+function linear_dependency_matrix(x::AbstractVector, w::AbstractVector, v::Real, D::AbstractMatrix)
   N = length(x)
   L = zeros(N,N)
   for i in 1:N
     Pi = map(x->legendre(x, i-1), x)
     for j in i:N
       Pj = map(x->legendre(x, j-1), x)
-      lij = integrate(w, v, D, Pi, Pj)
+      lij = integrate(Pi, Pj, x, w, v, D)
       L[i,j] = lij
       if i != j
         L[j,i] = L[i,j]
@@ -188,14 +231,14 @@ w ... GLGL weights
 v ... GLGL derivative weight
 D ... First derivative matrix based on x.
 """
-function expansion_matrix(x, w, v, D)
+function expansion_matrix(x::AbstractVector, w::AbstractVector, v::Real, D::AbstractMatrix)
   N = length(x)
   T = zeros(N,N)
   for i in 1:N
     li = map(xi->lagrange_poly(i, xi, x), x)
     for j in 1:N
       Pj = map(x->legendre(x, j-1), x)
-      tij = integrate(w, v, D, li, Pj)
+      tij = integrate(li, Pj, x, w, v, D)
       T[i,j] = tij
     end
   end
diff --git a/src/lgl.jl b/src/lgl.jl
index 493d7f4b1f562d4a60b63f475d6479adacc0a0dc..31309a4173066a78e23a5f26f5b05577a94ab6ff 100644
--- a/src/lgl.jl
+++ b/src/lgl.jl
@@ -11,7 +11,7 @@ Quadrature nodes and weights for Gauss-Legendre rule.
 n_pts ... number of quadrature points
 Exact up to order 2*n_pts - 1.
 """
-function rule(n_pts)
+function rule(n_pts::Integer)
   α = 0.0; β = 0.0
   x = zglj(n_pts, α, β)
   w = wglj(x, α, β)
@@ -21,16 +21,44 @@ end
 
 
 """
-Integrate with Legendre weight w(x) = 1.
-x ... Vector of LGL nodes
-w ... Vector of LGL weights
-f ... Vector with function values evaluated at nodes x.
+  integrate(f, x, w)
+
+Integrate `f` sampled on LGL quadrature nodes `x` with quadrature weights `w`.
+Use `rule` to obtain `x, w`.
 """
-function integrate(w, f)
+function integrate(f::AbstractVector, x::AbstractVector, w::AbstractVector)
   @toggled_assert length(w) == length(f)
   return dot(w, f)
 end
-function integrate(w, f::Function, idxs::AbstractVector{<:Integer})
+
+
+"""
+  integrate(f1, f2, x, w)
+
+Similar to `integrate(f, w)`, but computes the inner product of `f1*f2` where both
+`f1, f2` must be sampled on the LGL quadrature nodes `x` associated with the quadrature
+weights `w`.
+Use `rule` to obtain `x, w`.
+"""
+function integrate(f1::AbstractVector, f2::AbstractVector, x::AbstractVector, w::AbstractVector)
+  @toggled_assert length(w) == length(f1)
+  @toggled_assert length(w) == length(f2)
+  I = 0.0
+  for i = 1:length(w)
+    I += w[i] * f1[i] * f2[i]
+  end
+  return I
+end
+
+
+"""
+  integrate(f::Function, idxs::AbstractVector{<:Integer}, x, w)
+
+Integrate `f` sampled on LGL 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)
@@ -40,14 +68,23 @@ function integrate(w, f::Function, idxs::AbstractVector{<:Integer})
 end
 
 
-function integrate(w, f, g)
-  @toggled_assert length(w) == length(f)
-  @toggled_assert length(w) == length(g)
-  I = 0.0
-  for i = 1:length(w)
-    I += w[i] * f[i] * g[i]
+"""
+  integrate(f::Function, x, w)
+
+Integrate `f` sampled on LGL 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 I
+  return sum
 end
 
 
@@ -57,7 +94,7 @@ 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)
+function inverse_mass_matrix_exact(n_pts::Integer)
   z, w = rule(n_pts)
   N = n_pts - 1
   invM = zeros(n_pts, n_pts)
@@ -72,7 +109,7 @@ end
 """
 n_pts ... number of qudrature nodes = number of lagrange polynomials
 """
-function inverse_mass_matrix_lumped(n_pts)
+function inverse_mass_matrix_lumped(n_pts::Integer)
   _, w = rule(n_pts)
   invM = diagm(1.0 ./ w)
   return invM
@@ -84,7 +121,7 @@ end
 M_ij = (li, lj)
 n_pts ... number of qudrature nodes = number of lagrange polynomials
 """
-function mass_matrix_exact(n_pts)
+function mass_matrix_exact(n_pts::Integer)
   z, w = rule(n_pts)
   M = zeros(n_pts, n_pts)
   aux_z, aux_w = rule(n_pts+2)
@@ -94,7 +131,7 @@ function mass_matrix_exact(n_pts)
     li = L.(i, aux_z)
     for j = i:n_pts
       lj = L.(j, aux_z)
-      M[i,j] = integrate(aux_w, li, lj)
+      M[i,j] = integrate(li, lj, aux_w)
       if (i != j)
         M[j,i] = M[i,j]
       end
@@ -111,7 +148,7 @@ M_ij = (li, lj)
 Approximated mass matrix due to under integration of LGL quadrature.
 n_pts ... number of qudrature nodes = number of lagrange polynomials
 """
-function mass_matrix_lumped(n_pts)
+function mass_matrix_lumped(n_pts::Integer)
   z, w = rule(n_pts)
   M = diagm(w)
   return M
@@ -125,7 +162,7 @@ Exactly representable, because product is of order 2N - 1 which can be integrate
 exactly using LGL quadrature.
 n_pts ... number of quadrature nodes = number of lagrange polynomials
 """
-function stiffness_matrix(n_pts)
+function stiffness_matrix(n_pts::Integer)
   z, w = rule(n_pts)
   D = first_derivative_matrix(z)
   S = w .* D
@@ -141,14 +178,14 @@ of the numerical representation of the Legendre basis.
 x ... GLGL nodes
 w ... GLGL weights
 """
-function linear_dependency_matrix_lumped(x::tV, w::tV) where {tV}
+function linear_dependency_matrix_lumped(x::AbstractVector, w::AbstractVector)
   N = length(x)
   L = zeros(N,N)
   for i in 1:N
     Pi = map(x->legendre(x, i-1), x)
     for j in i:N
       Pj = map(x->legendre(x, j-1), x)
-      lij = integrate(w, Pi, Pj)
+      lij = integrate(Pi, Pj, x, w)
       L[i,j] = lij
       if i != j
         L[j,i] = L[i,j]
@@ -168,14 +205,14 @@ reverse.
 x ... GLGL nodes
 w ... GLGL weights
 """
-function expansion_matrix_lumped(x, w)
+function expansion_matrix_lumped(x::AbstractVector, w::AbstractVector)
   N = length(x)
   T = zeros(N,N)
   for i in 1:N
     li = map(xi->lagrange_poly(i, xi, x), x)
     for j in 1:N
       Pj = map(x->legendre(x, j-1), x)
-      tij = integrate(w, li, Pj)
+      tij = integrate(li, Pj, x, w)
       T[i,j] = tij
     end
   end
@@ -188,11 +225,11 @@ end
 """
 L_ij = (P_i, P_j)
 Analogon to mass matrix of Lagrange polynomials, expresses the linear dependency
-of the numerical representation of the Legendre basis. 
+of the numerical representation of the Legendre basis.
 x ... GLGL nodes
 w ... GLGL weights
 """
-function linear_dependency_matrix_exact(x, w)
+function linear_dependency_matrix_exact(x::AbstractVector, w::AbstractVector)
   N = length(x)
   N2 = N+2
   x2, w2 = rule(N2)
@@ -201,7 +238,7 @@ function linear_dependency_matrix_exact(x, w)
     Pi = map(x->legendre(x, i-1), x2)
     for j in i:N
       Pj = map(x->legendre(x, j-1), x2)
-      lij = integrate(w2, Pi, Pj)
+      lij = integrate(Pi, Pj, x2, w2)
       L[i,j] = lij
       if i != j
         L[j,i] = L[i,j]
@@ -221,7 +258,7 @@ reverse.
 x ... GLGL nodes
 w ... GLGL weights
 """
-function expansion_matrix_exact(x, w)
+function expansion_matrix_exact(x::AbstractVector, w::AbstractVector)
   N = length(x)
   N2 = N + 2
   x2, w2 = rule(N2)
@@ -230,7 +267,7 @@ function expansion_matrix_exact(x, w)
     li = map(xi->lagrange_poly(i, xi, x), x2)
     for j in 1:N
       Pj = map(x->legendre(x, j-1), x2)
-      tij = integrate(w2, li, Pj)
+      tij = integrate(li, Pj, x2, w2)
       T[i,j] = tij
     end
   end
diff --git a/src/spectralelement.jl b/src/spectralelement.jl
index 1fbbf23f6082e72dc20b1dd95e8386147b934756..dcf0b4186767f86b92bce7675df695bf51f1a907 100644
--- a/src/spectralelement.jl
+++ b/src/spectralelement.jl
@@ -93,9 +93,9 @@ end
   @unpack quadr = el
   @toggled_assert quadr in (:LGL, :LGL_lumped, :GLGL)
   if quadr === :LGL || quadr === :LGL_lumped
-    return LGL.integrate(el.w, u1, u2)
+    return LGL.integrate(u1, u2, el.z, el.w)
   else # quadr === :GLGL
-    return GLGL.integrate(el.w, el.v, el.D, u1, u2)
+    return GLGL.integrate(u1, u2, el.z, el.w, el.v, el.D)
   end
 end
 
@@ -104,9 +104,9 @@ end
   @unpack quadr = el
   @toggled_assert quadr in (:LGL, :LGL_lumped, :GLGL)
   if quadr === :LGL || quadr === :LGL_lumped
-    return LGL.integrate(el.w, u)
+    return LGL.integrate(u, el.z, el.w)
   else # quadr === :GLGL
-    return GLGL.integrate(el.w, el.v, el.D, u)
+    return GLGL.integrate(u, el.z, el.w, el.v, el.D)
   end
 end
 
@@ -115,9 +115,9 @@ end
   @unpack quadr = el
   @toggled_assert quadr in (:LGL, :LGL_lumped, :GLGL)
   if quadr === :LGL || quadr === :LGL_lumped
-    return LGL.integrate(el.w, fn, idxs)
+    return LGL.integrate(fn, idxs, el.z, el.w)
   else # quadr === :GLGL
-    return GLGL.integrate(el.w, el.v, el.D, fn, idxs)
+    return GLGL.integrate(fn, idxs, el.z, el.w, el.v, el.D)
   end
 end
 
diff --git a/test/UnitTests/src/test_spectral.jl b/test/UnitTests/src/test_spectral.jl
index 298c3822905db49d3a3da6fc0aa86a9aecab7fbb..3a2a38601b4d6e46d8767170b3bf5bf527f0c86d 100644
--- a/test/UnitTests/src/test_spectral.jl
+++ b/test/UnitTests/src/test_spectral.jl
@@ -162,4 +162,21 @@ end
 end
 
 
+@testset "More integrator tests" begin
+
+  n_pts = 6
+  p = x -> (x-4)^2 * (x+1)^3
+
+  x, w = dg1d.LGL.rule(n_pts)
+  int = dg1d.LGL.integrate(p, x, w)
+  @test int ≈ 140/3
+
+  x, w, v = dg1d.GLGL.rule(n_pts)
+  D = dg1d.first_derivative_matrix(x)
+  int = dg1d.GLGL.integrate(p, x, w, v, D)
+  @test int ≈ 140/3
+
+end
+
+
 end