We are trying to use the high-order Lagrange VTK cells first (?) presented in this blog post. We’ve had some success in getting curved 4th-order VTK_LAGRANGE_QUADRILATERAL cells to be visualized in ParaView, but when comparing it on a rectilinear domain against the result of a classical, piece-wise constant VTK_PIXEL mesh, the results are somewhat underwhelming. We thus figure there might be something wrong with the way we create the Lagrange cells.
Just to recap, we create each VTK_LAGRANGE_QUADRILATERAL with 4 by 4 points, as shown in the blog:
However, the resulting images suffer from weird “spots” and a clear degradation of the quality near the element boundaries (here on a 2-by-2 mesh, each with a 4-by-4 Lagrange quad):
Is there anything that you need to be aware of when creating these cells besides the order of the points? For example, we are using the Legendre-Gauss-Lobatto nodes as the points for visualization - is this correct or would we have to use a different set of nodes (e.g., Chebyshev-Gauss-Lobatto, equidistant etc.)? Any kind of input would be highly appreciated!
@ssss Since you seem to be quite experienced with the use of high-order cells (as evidenced here), would you by any chance also know the answer to this thread?
vtkLangrange* cells use equi-distant nodes (in the reference element). Therefore, it is (usually) necessary to interpolate solution data to such points. You may want to take a look at what we recently developed in PyFR to output high-order vtkLagrange cells using FR-DG simulations.
@ssss How should the nodes look in a curved cell, then?
Once we allow non-Cartesian grids, we could define a mapping that maps the equidistant nodes in the reference element to any nodes in the physical domain.
Suppose we have one element, which maps onto [-1,1]^2 in the physical domain. We run two simulations with different mappings. One with the identity mapping, one with a mapping that transforms the element such that we still get the same physical domain, but our nodes get mapped to different nodes (for example LGL nodes).
We evaluate our initial condition at the physical coordinates of our nodes, so we obviously get different nodal values in both simulations. However, both mappings of the interpolation polynomials should be identical in the physical domain.
From VTK you get the parametric coordinates of the VTK reference element nodes.
In your code you interpolate your metric polynomials to the the parametric coordinates of the VTK reference node (beware, VTK uses [0, 1] as parametric coordinates I think) resulting in the physical coordinates of the VTK cell node.
