# Why is FEM theory only ever written down for simplices?

I noticed that in most theoretic literature (Braess, Ciarlet, ...) for Finite Elements, the method is only described in a detailed way for simplex triangulations, while other forms of hexahedral methods are often omitted or only briefly mentioned.

Is there any theoretical differerence between FEM theory of simplices and hexahedrals? Do theoretical results (convergence, error-estimates) not carry over to hexahedrals, or is it just that the generality is so obvious that it is left to the reader.

In my opinion the basics stay the same but there are some differences worth commenting upon.

The basic principles of defining a local coordinate system, integrating local dense contribution matrices using quadrature rules, and assembling them into global sparse matrices are identical. Often, hexahedral basis functions arise from taking tensor products of some underlying one-dimensional basis, and are integrated using tensor products of one-dimensional quadrature rules, which exposes some opportunities for optimizing code.

One notable difference is polynomial completeness. High order tetrahedral basis tend to arise from some sort of pascal's triangle / pascal's tetrahedron arrangement. For instance, a 2D gradient-conforming basis needs to span {$$1$$, $$x$$, $$y$$, $$x^2$$, $$xy$$, $$y^2$$} to be complete to second order, and these 6 functions can be arranged very tidily on a 6-node triangle (3 vertices, 3 edge midpoints). In contrast, the most natural way to write a second order basis for a 2D quadrilateral is via a tensor product of two 1D lines, leading to a span of {$$1$$, $$x$$, $$y$$, $$x^2$$, $$xy$$, $$y^2$$, $$x^2y$$, $$xy^2$$, $$x^2y^2$$}. This contains more functions / requires more work, but is not any more accurate in the asymptotic sense. This over-completeness becomes more exaggerated as you raise the order or migrate from 2D to 3D.

Despite this overcomplete-ness at the element level, hexahedral elements are generally more efficient in total cost, simply because they are "larger" and fill space more efficiently. Model accuracy is typically a function of mesh size $$h$$. A hexahedron with sides of length $$h$$ takes up $$h^3$$ volume, while a tetrahedron with sides of length $$h$$ only takes up a volume of $$h^3/6$$. So given some computational domain with volume $$V$$ and an desired accuracy/mesh size $$h$$, I expect to need 6x as many tetrahedra than hexahedra. Even though this difference is just a constant factor, we're talking about the size of the input here, so it will be further magnified by any downstream algorithmic complexities (for instance, if you're using solver/preconditioning techniques that scale superlinearly).

I do want to comment briefly on accuracy. A surprising aspect of some basis functions is that element shape / distortion can effect their order of accuracy, because of how they depend on the pullback/jacobian. As a concrete example, consider 2D Nedelec (curl-conforming) functions. The curl of these functions (ie what's going into the "stiffness" matrix, to borrow a term) is proportional to the jacobian. For a triangle, the jacobian is constant, so you get curl complete to a constant (this is what you want, a lowest order "$$\star$$-conforming" element should have "$$\star$$" complete to a constant). In contrast, the jacobian for a quadrilateral is only constant if the shape mapping is affine (ie the quadrilateral is a parallelogram), which means this element is non-convergent on general meshes (ie you can refine $$h$$ but the error will not improve). Raviart-Thomas elements are similarly affected, as they have a similar dependence on jacobian. Note that a similar thing happens (non-constant jacobian / non-convergence) if you use curvilinear shape mappings (for either triangles or quadrilaterials), I only call out this case specifically because it's surprising that a linear/straight-line quadrilateral can exhibit this problem too (when it's not a parallelogram).

All this said, I find the biggest obstacle to using hexahedra for FEM are not formulation or convergence issues, but mesh generation and computational geometry issues. There are robust general-purpose algorithms for automatic tetrahedral mesh generation (delaunay refinement in particular), but hexahedral meshing is more brittle and often requires user interaction to decompose complicated geometry into simple parts that are eligible for structured strategies or pave/sweep strategies. Note that you can fuse tetrahedral grids with hexahedral grids using pyramidal elements, but the basis functions on pyramids can be pretty technical (pyramids are often implemented in terms of degenerate/distorted hexadra, and the distortion function ends up contributing non-polynomial terms to the jacobian/basis functions, which can make it difficult to obtain the correct order of accuracy).

In light of these geometrical complexities, I find it reasonable for papers to focus more on simplicial elements, there's just a larger market for them.

• Wow, thank you very much for the elaborate answer! – MPIchael Apr 22 at 14:54
• The Finite Element Method is blessed with regards to meshes, since it performs well with tetrahedral meshes, combined with the fact tet meshing is quite easy and there are plenty of tet meshers out there. If you are doing CFD, i.e. fluid mechanics, hex meshes behave better than tets, although this is field/application dependent. If you are searching for hex meshers, preferably open source, your search is a short one. – Dohn Joe Apr 22 at 15:06

Most FM methods require that we can identify pairs of adjacent elements. This is easy with elements derived from a quadrilateral/hexahedral grid that can be mapped onto an array structure. To deal with the complex domains typical of aircraft, for example, the "multiblock" structure can be adopted where the domain is first mapped onto "blocks" which are very coarse quad/hex grids with an arbitrary connectivity between them. Generating these coarse grids is hard to automate if there must be many blocks, and is commonly done by hand, which is time-consuming. Because of this the multiblock strategy is now little-used in practice, and a single "unstructured" quad/hex grid is used for the entire domain. Individual elements should still be as "regular" as possible in aspect ratio and skewness.

• I fail to see how this answer the OP question. – nicoguaro Apr 25 at 20:03