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

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.

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.

• Thanks, that already shed some light on my question. I updated my question considering your edit. So in my case (where I really want the numeric solution and not the exact solution visualized), would I still have to numerically (or analytically) invert the Vandermonde, or is there a better way out? Also, what I don't understand from your equation: how is it possible for the stated equation (u_eq = V_eq V^-1 u) to hold if u_eq and u are vectors of different size (N != M) and the Vandermondes being square matrices? Commented May 14, 2013 at 18:43
• To answer your question: $\mathbf{V}$ is a $N\times N$ matrix, and $\mathbf{V}_{eq}$ is $M \times N$, so your computed solution comes out $M \times 1$. Also: Do you have an unstructured mesh here and you are trying to get a global equally spaced sampling of your solution? Or are you trying to get equally spaced nodes on each individual element? Also if you are using Lagrange polynomials directly then the first Vandermonde matrix is just the identity. Commented May 14, 2013 at 18:48
• Ah, right. I could have deduced that from your second equation, since $$i = 1,...,M$$ and $$j = 1,...N$$. And yes, I am trying to get the equally spaced nodes on each individual element. Commented May 14, 2013 at 19:11
• In that case the Vandermonde matrix you're using will be block diagonal. You can invert it once on the reference element and just apply that elementwise. It should be a pretty small matrix and easily inverted (even if you are using a massive polynomial order, this is feasible). Since you are using the Lagrange basis directly, it's the identity matrix anyways. But you'll still need the Vandermonde for the equally spaced points. No inversion needed there. Commented May 14, 2013 at 19:27

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.