How to properly use polynomial projection to get values at visualization nodes?

Question

I am trying to implementing a nodal discontinuous Galerkin spectral element method for linear and non-linear systems of equations. The solution at each time step is given at N nodes, which are located at the integration points for Legendre-Gauss/Legendre-Gauss-Lobatto quadrature.

For visualization purposes I would like to get the solution at M equidistant visualization nodes, where typically M > N. As far as I know, generally L2 projection is preferred to polynomial interpolation (albeit being more expensive), especially when going to higher degrees of freedom, since the projection gives the best approximation w.r.t. the L2 norm.

As far as I understand it, this is a rough outline of the steps I have to take:

1. Assemble the Vandermonde matrix V for the old nodes.
2. Invert V.
3. Use V and the nodal values to get modal coefficients.
4. Project to new nodes.
5. Assemble Vandermonde matrix V' for the new nodes.
6. Use V' and the modal coefficients to get the nodal values.

Finally, my question(s):

• a) Is this approach generally correct?
• b) Are there algorithmic simplifications for my specific problem, e.g.
  • the inversion of V is trivial because I can use some properties of the used polynomials (Legendre) or nodes (L-Gauss/L-Gauss-Lobatto)?
  • the whole algorithm can be simplified because I already start out with a polynomial approximation and not an arbitrary exact function?

Any suggestions or pointers to further reading material are highly appreciated.

Edit:

I am using Lagrangian basis functions for my computations, the GL/LGL nodes only come in for the purpose of numerical integration. So at the point where I want to get the solution at the new nodes, I do not have any information on the solution but the values at the old nodes.

Answer 1

You want to interpolate here, not perform a $L^2$ projection. You want the value of the DG solution at equally spaced nodes. In any case a $L^2$ projection will just feed you back the same polynomial you started with (remember, it minimizes $L^2$ error.. if you're already a polynomial, then you can't do any better).

Here are the matrices you'll need to compute

1. Vandermonde for LGL nodes: $\mathbf{V}$
2. Vandermonde for desired equally spaced noes: $\mathbf{V}_{eq}$

Then given a solution vector $\mathbf{u}$ the interpolation will look like

$$\mathbf{u}_{eq} = \mathbf{V}_{eq} \mathbf{V}^{-1} \mathbf{u} $$

All this assumes that you are not directly using the Lagrange polynomial basis for computations, but instead an orthogonal basis like Legendre polynomials, so your Vandermondes will be

$$\mathbf{V}_{ij} = P_j(x_i)$$

in whatever basis you're using for the computations.

Edit: Reading that over again It isn't clear whether you're trying to visualize a computed DG solution or the analytic solution (or initial condition). In the latter case, still use interpolation. It's much easier, and you won't notice the accuracy difference in a visualization. I'd keep it to the LGL nodes though and then sample the interpolating polynomial for visualization purposes.

Answer 2

In my codes, when I get the FE solution $f(x)$ as a linear combination of the basis functions $\phi_i(x)$: $$ f(x) = \sum_i q_i \phi_i(x) $$ I can easily evaluate it at any point $x_0$ by evaluating the basis functions at this point ($\phi_i(x_0)$) and doing the linear combination: $$ f(x_0) = \sum_i q_i \phi_i(x_0) $$

I also use this approach to transfer the solution to quadrature points (in case they are not the same as the nodes). Obviously one of the advantages of spectral elements is to use the same points as nodes and quadrature points, then one doesn't have to do that --- but you still need this for plotting (as you asked), though again, many times I save myself the trouble and just plot the values at the nodes, that usually gives me a pretty good idea how the solution looks like.