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 system matrix for the following cases:

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

So IIP is set by default in the compressible flow solver. See issue #306 for details. This MR provides a fix to this issue.

Proposed solution

In the SIP method, we must evaluate the trace and symmetry flux contributions using the same trace quadrature formula. The following cases do not satisfy the condition:

1. p-refined mesh. Then two adjacent elements have different quadrature orders;
2. Even in the uniform order mesh, since Tri/Tet/Pyr/Prism use Gauss-Radau points, the two adjacent traces may still have different quadrature types (GR vs GLL).

In theory, if we can build another trace space, and perform quadrature in the new space for both left and right flux evaluation, then we can fulfill the condition. Such facilities already exist, that is DisContField::m_trace, a shared trace space between adjacent elements. However, to correctly use trace space for flux evaluation, a lot more work needs to be done, which is still on the way.

So instead, we still perform quadrature in the local trace space. If we can set sufficient local quadrature points that are accurate for both sides, then we can get a symmetric system as well. (This approach is unofficially called Point-to-Point interpolation**)**

The original AddSymmFluxIntegralToCoeff perform the tasks in the following steps:

1. First, multiply trace quadrature metric WJ
2. Then interpolate trace data to left/right local trace space
3. Finally do IProduct.

This is wrong! The quadrature formula becomes invalid after the interpolation. The correct way is to first

1. interpolate trace data to local trace space
2. then multiply local trace quadrature metric WJ
3. finally do IProduct

Only fixing that bug is not enough To fully fix GR/GLL mismatch in Tet/Prism/Pyr/Tri, we suggest switching to Gauss-Legendre (GL) points as shown below

The key is that GL points are rotational symmetric to avoid any mismatch shown above.

Implementation

The changes to the codes can be split into two parts:

1. Let GaussLegendre points works in DG;
2. Fix AddSymmFluxIntegralToCoeff by local trace quadrature metric;

In the present DG framework, we extract trace physical data from the field physical space, which requires the quadrature Point to include the endpoints. To fulfil this, we introduce a new combined Point type: Gauss-Legendre with endpoints ( GaussLegendreWithMP and GaussLegendreWithM). The quadrature weights are 0, w1, w2, w3,...wN,0 . It naturally fits into the current framework.

The present implementation is not performance-orientated. We will shortly provide an efficient matrix-free implementation in a later MR.

Tests

A case with a quarter annulus domain and prism mesh is provided to verify the new code can resolve the issue.

rhow now have a much lower magnitude than before (1e-6 vs 1e-2).

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).