I would like to ask how is the evaluation of integrals over inter-element interfaces implemented in a typical DG code. I can think of two basic approaches (I assume 2D mesh here):
1) Perform 1D interpolation of the interface data seen from 'left' and 'right' of each edge and numerically integrate. For this, I don't need to evaluate the element shape functions in the elements adjacent to the edge. This will probably work fine for Lagrange shape functions (element interpolant restricted to edge and 1D edge interpolant should be identical).
2) Interpret the 1D edge quadrature rule as a quadrature with points placed on the boundary of each element and use the element shape functions to compute edge integrals.
Option 1) does not require extra storage for jump/average terms: evaluation on each edge is performed on the fly and I can directly add boundary integrals to the weak form. Option 2) will require a lot of extra memory allocation for interface integrals: jump and average terms are only complete after I visit both left and right neighbour of each edge and temporary data evaluated at quadrature points have to wait somewhere before the interface integrals are assembled. I could also loop over all edges and obtain directly the left-and right interpolant from the element shape functions, but that would mean that I re-evaluate the element shape functions every time I visit one of its edges (computational overhead).
I think that the disadvantage of 1) is that it won't be possible to evaluate shape function gradients on element interfaces for higher-order equations. In addition, it would probably not work with some modal bases (all modes from the element interior are needed in the expansion).
I cannot decide which approach is better (or if they are good at all). Maybe there's a more clever way I did not think of. Can you share your opinion?