##
Quadrature rules for N-th order consistency of discontinuous Galerkin on rectangles

We want to solve

∂_{t}u + Div(f) = s.

Multiply by a test function p and integrate by parts
over a mesh cell:

∂_{t}(∫ u*p) + ∫_{∂}(**n**•f*p) = ∫(f•∇ p) + ∫(s*p)

Restrict p to N-th order polynomials. Then we have an exact
unclosed system of equations for the coefficients of the Taylor
expansion of u out to order N. Replacing f with a numerical flux
function F (and specifying a numerical integration rule, which
can be viewed as part of specifying F if we project F onto the
subspace for which the numerical integration rule is exact)
and likewise replacing s with a numerical source term S
gives a closed system:

∂_{t}(∫ u*p) + ∫_{∂}(**n**•F*p) = ∫(F•∇ p) + ∫(S*p).

We seek conditions to ensure that the error
is order M:=(N+1), i.e. the maximal order of accuracy
that the representation space of the solution is capable of.

A k-point integration rule (Gaussian Quadrature)
is needed to exactly integrate
one-dimensional polynomials of order 2*k-1. That is, to exactly
integrate a one-dimensional polynomial of order κ you need
an integration rule with ceiling((κ+1)/2) many points.
If we use a tensor product integration rule for a rectangular cell
with n dimensions then we need M^n points.

Taking F to be an N-th order polynomial (projection),
to exactly integrate each term in a cartesian mesh cell
in n-dimensional space we need

term order of integrand number of points needed
============ ================== =======================
∫_{∂}(**n**•F*p) 2N (N+1)^{n-1}
∫(F•∇ p) 2N-1 N^{n}
∫(S*p) 2N (N+1)^{n}

This number of points is in general necessary,
since the number of weights plus points equals the
number of coefficients of the polynomial basis
in each case.