Non linear quads have poor triangulation
This issue is related to issue #17961. Apart from the bug in vtkQuad::IntersectWithLine and vtkQuad::Triangulate, non linear quads perform an arbitrary triangulation relying on the one of vtkQuad, which should be unrelated. If we assume that a vtkQuad should be planar, then performing a Delaunay triangulation makes sense. However, the triangulation of a non planar quad should depend on the local convexity of the shape where the quad lies.
There are 2 ways to triangulate a non planar quad. Given a normal orientation, the triangulation results in a convex or a concave cell. Its convexity is only determined by the splitting diagonal choice. Here are visualization of such cases
When subdividing a non linear quad, each underlying refined quad needs to share the same convexity as the local convexity at the corresponding position in the non linear quad. In other words, the splitting diagonal needs to be chosen case by case, and has no causality with the triangulation of a linear quad. I didn't write the equations down, but we need to decide the triangulation based on sign of the estimate of the mean curvature (or the laplacian equivalently). Ideally we should look at the integrated value to smoothen the estimate, but it might not be necessary in practice.
This error in the triangulation is detrimental when visualizing non linear cells, as when the wrong triangulation is performed, you see bumps caused by the wrong local convexity of the shape, as you can see below. You see it by looking at the normals, which oscillates at every subdivision.
I think that this issue should be fixed at the same time as #17961, as currently non linear quads rely on vtkQuad. I haven't looked if there is the same issue for every non linear cells. This should be looked at as well.