Compute mass matrices accurately for tensor products

[inducer]

Legendre-Gauss-Lobatto with $n$ points integrates polynomials up to degree $2n-3$ exactly. (I.e. for $p=1$ aka $n=2$, it's the trapezoidal rule, accurate to degree 🥁 1.) That's what we're currently using to form reference mass matrices, by 🤦 automatic generalization from simplices. (I.e. I don't think anybody has looked at that code with tensor product goggles on.)

To do even just the reference bilinear form right, we need a degree of exactness of $2(n-1)=2n-2$ (because both polynomials have degree $n-1$). Legendre-Gauss is exact to degree $2n-1$, so that would work.

For Legendre-Gauss points, the 1D nodal mass matrix is just the diagonal matrix with the Gauss weights (because Gaussian quadrature is exact). That's a useful computational shortcut which we should exploit, but out of scope for this issue. I'll file a separate one.

It is common to use diagonal mass matrices also for LGL even though they are not exact in that case, this is called "lumping". (e.g. ref) That's mainly helpful for the inverse, cf. #343.

This affects time-stepping solvers via the mass inverse and the surface mass matrices.

This needs to be fixed in two spots:

• the actual mass matrix, and
• the face mass matrix.

cc @a-alveyblanc @MTCam @majosm @lukeolson