Quadrature order for finite elements and time dependent discontinuous Galerkin

Question

When setting up a finite element system you have to use quadrature to calculate the integrals. I'm having trouble understanding what order rule to use.

I know of some rules of thumb, for example with the 1D Hemholtz you end up with the matrix $M + S$, where $M_{ij} = (\phi_i, \phi_j)$ and $S_{ij} = (\phi_i', \phi_j')$. In this case if $\phi$ is a polynomial of degree $n$, then the integrals in $M$ will of degree $2n$ and the integrals in $S$ will be of degree $2(n-1)$. Using a Gauss rule of order $n$ is recommended here, which integrates $S$ exactly.

When using DG on the advection equation the "stiffness" matrix becomes $S_{ij} = (\phi_i', \phi_j)$. In this case I found that using an order $n$ Gauss quadrature rule (which should still integrate $S$ exactly) gives me a singular matrix, and that I have to use a rule of order $n+1$ to get an invertible matrix. Is there some way to know beforehand what order rule to use, or it just determined experimentally?

EDIT: As per comments below, I'm using DG with Langrangian elements on a uniform grid to solve the time dependent advection equation $u_t + au_x = 0$. I have periodic boundary conditions on the domain $[0, 2\pi]$ and the initial condition $u(x,0) = \sin(x)$.

This is my first attempt coding DG myself, so I'm using central flux and backwards Euler for the time step. Using basis functions of degree $n$ and $\Delta t = h^{n+1}$, the solutions converge to the exact solution in the $L^2$ norm as $O(h^{n+1})$.

Answer 1

Since you have periodic boundary conditions, the system matrix (arising from the volume and flux terms) cannot have full rank (the solution of the stationary equation is unique only up to constants). Hence the well-posedness of the backward Euler step relies on the invertibility of the mass matrix (or at least that the two kernels have only a trivial intersection).

Now you can fairly easily show that the mass matrix will be positive definite if the quadrature is exact. As the product of two polynomials of order $n$ is $2n$, Gauß quadrature of order $n$ is not sufficient and you need at least order $n+1$ (as you have already noted).

The converse is not as obvious, but there's a bit of "finite element folklore" (at least in the engineering community) that each Gauß quadrature point adds $d$ (the space dimension) to the rank of the elemental mass matrix, and full rank is equal to the number of degrees of freedom in each element. In 1D, this is $n+1$, so you need $n+1$ points.

Usually, it is sufficient that the highest-order term (i.e., the one in the bilinear form involving the hightest-order derivatives) is integrated exactly (ignoring variable coefficients) and surface terms are integrated using one order higher, but if the highest-order term is itself not invertible (e.g., because of pure Neumann or periodic boundary conditions), you need to make sure the mass matrix is integrated exactly as well.