How do hexahedral FEM meshes improve approximation quality per degree of freedom, compared to tetrahredal meshes?

Question

From the deal.II FAQ :

...quadrilaterals and hexahedra typically provide a significantly better approximation quality than triangular meshes with the same number of degrees of freedom; you therefore get more accurate solutions for the same amount of work...

How can this point be made more precise?

Answer 1

Having written this entry, let me also answer you question :-)

The issue in essence boils down to the following: given a function $u(x)$ and its interpolation $u_h(x)$ onto either a triangular or quadrilateral mesh with the same number of unknowns, then which of the two is more accurate? In other words, is $$ \| u - I_h^\text{triangles} u \| < \| u - I_h^\text{quads} u \|, $$ or is it the other way around. We know that they have the same asymptotic order ${\cal O}(h^{p+1})$ when using polynomials of order $p$ (in $L_2$). The question is simply which of the two has the smaller interpolation constant. The answer, in practice, is that the inequality above is wrong in general, and that the following is empirically true: $$ \| u - I_h^\text{triangles} u \| \approx 10 \| u - I_h^\text{quads} u \|. $$ In other words, the extra polynomial $xy$ on quadrilaterals allows for a 10 times better approximation. (I'm sure you can find counterexamples; this is just an empirical observation for "typical" solutions of PDEs.)

Answer 2

I fully agree with Wolfgang's answer, but I'd like to add something of my experience. I had recently to choose between tri and quad based FEM and after googling for it, I decided that I would code both from scratch for a fair comparison (I had a previous experience in writing FEM kernel. I used C++ and the Eigen library for matrices and vectors. I considered only cartesian meshes).

The answer by Wolfgang seemed to hold in my cases (simple advection-diffusion problems) as well, so I can confirm this. However, the "same amount of work" is not exactly true. There are other differences that we need take into account for a true time to error comparison.

• The number of quadrature points might be different. In my case, the tri quadrature had 1 point while the quad one had 4 points. That makes twice less quadrature nodes in the assembly procedure for the tri based FEM. However, this depends on the problem at hand.
• The number of non zero entries in the resulting matrix is different for the two approaches. For a cartesian-like mesh with triangles, a internal degree of freedom typically couples with 6 others, while with quad base, it does so with 8 others. That makes the assembly more expensive and might also affect the linear solver depending on which one is used.
• The linear system for the quad-based FEM was (suprisingly?) easier for the iterative solver than the tri-based (less iterations for a given tolerance). I used an ILUT-preconditioned Krylov-based solver (can't remember which one exactly). I don't know if this is general or specific to the problem, but that is what I saw.

In the end, I still found that quad-based FEM was faster for the time-to-error comparison, but the factor was definitly less than 10, around 3 if I remember correctly (This number of course depends on the specific implementations and problems).

Another point for the quad-based FEM is that it is WAY easier to implement if you want to consider high order or "high" dimensional (>2) FEM.

Hope this is useful!