Draft: Fix incorrect symmetry flux for interior penalty methods in DG

Issue/feature addressed

The current implementation of AddSymmFluxToField has defects and leads to a non-symmetrical matrix system for

1. generalized tensorial elements (Tri/Tet/Pyr/Prism) and
2. variable quadrature orders (trace left element has different points from right element) and causes SIP (symmetrical interior penalty) method to fail.

We have to set IIP by default in the compressible flow solver. See issue #306 for details. In this MR, I provide a fully functional solution to this issue.

Cause of the problem

In the IP method, to ensure the symmetry of the system, both trace flux integral and symmetry flux integral should do IProduct in the same trace space, in other words, we must guarantee the quadrature formula is identical, for both flux terms and for both side of elements.

For variable quadrature order cases, left and right local traces will have different quadrature points. In the current design, we create another trace Explist, which is identical to either the left or right local trace. We provide InterpLocTraceToTrace and InterpTraceToLocTrace to transfer the data between local traces and the trace space.

If we first interpolate trace data to both sides and then apply their own quadrature formula respectively, we definitely get different quadrature results; If apply the same quadrature formula, e.g. by multiplying the trace data by trace quadrature weights, and then interpolate to left/right side, the quadrature formula will be no longer valid and the results will be different as well. That's why SIP fails in variable polynomial cases.

For Tri/Tet/Pyr/Prism, it's cumbersome to evaluate the derivative at the collapsed vertex in Tri/Tet elements. So we choose Gauss-Radau points to avoid it. When constructing the trace space, we restore 'n' Gauss-Radau points back to 'n+1' Gauss-Lobatto points just by adding the end point. So basically we only need to interpolate the end point when transferring data between local and shared traces.

This also violates the quadrature formula: removing/interpolating the end point also makes the quadrature rule invalid.

Solution and Implementation

The key ideas are summarized as follows:

1. Use Gauss-Legendre points instead of Gauss-Lobatto points or Gauss-Radau points, to solve the symmetry issue in generalized tensorial elements;

In this way, both local trace and trace can use the same Gauss-Legendre points. No interpolation happens between them. Quadrature results are identical.

2. Create so-called trace "Conformal" element expansions for variable polynomial orders, to avoid InterpTraceToLocTraces

"Conformal" means the new element has the same points as the target trace, but still keeps the same coeff space as the original element. We replace InterpTraceToLocTraces with a new function called InterpTraceToConformalTrace which only involves copy operations. In this way, the same quadrature rule is guaranteed for both sides.

Other implementation details:

• DG method requires the evaluation of physical value on traces, but inherently the Gauss points don't include the end points (trace points). To avoid adding another operation which deliberately evaluates the trace phys. We create a new point set called: GaussLegendreWithMP, which means GL points plus M(x=-1) and P(x=+1) end points. The quadrature weights are defined as 0,w0,w1,...wN-2,0. Therefore, when performing the IProduct, it automatically restores the GL quadrature without any special treatment. I believe it is a clever solution to realize my idea in the current Nektar++ DG framework with minimum effort.

• We provide a matrix-type operator for evaluating symmetrical flux integral and adding it to field space, which may also be called Lifting operator. It is very similar to the GJP implementation. The matrix PhysDerivBaseOnTraceMat will multiply D * B to the conformal trace phys space and the results will be field coeff space. Similarly, we also provide the matrix PhysBaseOnTraceMat which can be used to build an operator to replace AddTraceIntegral. These new matrix-type operators, although they may be inefficient, can be applied to any basis type, Modified, GLL_Lagrange or Orthogonal, regardless of boundary/interior decomposition.

It's possible to convert it into matrix-free operations but this will be implemented in a later MR.

Tests

Suggested reviewers

Please suggest any people who would be appropriate to review your code.

Notes

Please add any other information that could be useful for reviewers.

Checklist

• Functions and classes, or changes to them, are documented.
• User guide/documentation is updated.
• Changelog is updated.
• Suitable tests added for new functionality.
• Contributed code is correctly formatted. (See the contributing guidelines).
• License added to any new files.
• No extraneous files have been added (e.g. compiler output or test data files).