Issues with Pyramid elements

This issue is for discussing several problems with pyramid elements. While these problems result from the interaction of dune-uggrid, dune-geometry, dune-localfunctions they become most visible here. Below only the case of Lagrange elements is discussed, but other types may suffer from the same issues. The current situation is as follows:

  • Dune-uggrid implements the pyramid elements by a picewise Q1 geometry transformation by splitting into two tetrahedra.
  • Dune-geometry provides quadrature rules for pyramid elements that are obtained by collapsing tensorial cube rules.
  • Dune-localfunctions provides two types of Lagrange elements, static and dynamic order ones:
    • LagrangePyramidLocalFiniteElement supporting compile time order k=1,2 is a composite element that splits the pyramid into two tetrahedra and uses continuous piecewise Qk shape functions.
    • LagrangeLocalFiniteElement with pyramid geometries implements shape functions of any order by a Duffy transformation. For higher order this suffers from severe accuracy issues unless higher precision types are used that make the code significantly slower. It also is neither thread-safe nor cheaply copyable.

As a consequence, when using the out of the box quadrature rules, one will in general not be able to integrate exactly. While the quadrature rules and dynamic order Lagrange shape functions seem to be compatible, at least the geometry Jacobian may be discontinuous along the diagonal. Even worse than being inaccurate, the tensorial rules contain points on the diagonal, i.e. we evaluate shape function or geometry Jacobian exactly on its jump.

There are multiple possibilities to improve the situation, e.g. one could directly implement the Duffy-based Lagrange basis to solve the accuracy problems. But then the quadrature issue wrt to the geometry jacobian persists. Conversely, if one decides to adjust the default quadrature rules to use a tetrahedral split, they will not be able to integrate the Duffy-transformed rational functions exactly. So whatever we may try to improve the situation should keep in mind geometry, quadrature, and shape functions at once.

Edited by Carsten Gräser