Compute mass matrices accurately for tensor products

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 :drum: 1.) That's what we're currently using to form reference mass matrices, by :facepalm: 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](https://doi.org/10.1002/fld.3923)) 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