Hp-FEM


hp-FEM is a general version of the finite element method, a numerical method for solving partial differential equations based on piecewise-polynomial approximations that employs elements of variable size
' and polynomial degree '. The origins of hp-FEM date back to the pioneering work of Ivo Babuška et al.
who discovered that the finite element method converges exponentially fast when
the mesh is refined using a suitable combination of h-refinements
and p-refinements. The exponential convergence makes the method
a very attractive choice compared to most other finite element methods
which only converge with an algebraic rate. The exponential convergence
of the hp-FEM was not only predicted theoretically but also observed
by numerous independent researchers.

Differences from standard FEM

The hp-FEM differs from the standard FEM in many aspects.
The Fichera problem is a standard benchmark problem for adaptive FEM codes. One can use it to show the dramatic difference in the performance of standard FEM and the hp-FEM. The problem geometry is a cube with missing corner. The exact solution has a singular gradient at the center. The knowledge of the exact solution makes it possible to calculate the approximation error exactly and thus compare various numerical methods. For illustration, the problem was solved using three different versions of adaptive FEM: with linear elements, quadratic elements, and the hp-FEM.
The convergence graphs show the approximation error as a function of the number of degrees of freedom. By DOF we mean parameters that are needed to define the approximation. The number of DOF equals the size of the stiffness matrix. The reader can see in the graphs that the convergence of the hp-FEM is much faster than the convergence of both other methods. Actually, the performance gap is so huge that the linear FEM might not converge at all in reasonable time and the quadratic FEM would need hundreds of thousands or perhaps millions of DOF to reach the accuracy that the hp-FEM attained with approximately 17,000 DOF. Obtaining very accurate results using relatively few DOF is the main strength of the hp-FEM.

Why is hp-FEM so efficient?

Smooth functions can be approximated much more efficiently using large high-order elements than small piecewise-linear ones. This is illustrated in the figure below, where a 1D Poisson equation with zero Dirichlet boundary conditions is solved on two different meshes. The exact solution is the sine function.
While the number of unknowns is the same in both cases, the errors in the corresponding norm are 0.68 and 0.20, respectively. This means that the quadratic approximation was roughly 3.5-times more efficient than the piecewise-linear one. When we proceed one step further and compare four linear elements to one quartic element, then both discrete problems will have three DOF but the quartic approximation will be approximately 40-times more efficient. When performing few more steps like this, the reader will see that the efficiency gap opens extremely fast.
On the contrary, small low-order elements can capture small-scale features such as singularities much better than large high-order ones. The hp-FEM is based on an optimal combination of these two approaches which leads to exponential convergence. Note, that this exponential convergence is expressed in axis of error vs degrees of freedom. For real-life applications we usually consider computational time needed to reach the same level of accuracy. For this performance indicator h- and hp-refinement can provide similar results, e.g. see the final figure at . As soon as it is harder to program and parallelize hp-FEM compared to h-FEM the convergence excellence of hp-refinement may turn to be unpractical.

What is hp-adaptivity?

Some FEM sites describe hp-adaptivity as a combination of h-adaptivity and p-adaptivity. This is not entirely accurate. The hp-adaptivity is significantly different from both h- and p-adaptivity since the hp-refinement of an element can be done in many different ways. Besides a p-refinement, the element can be subdivided in space, but there are many combinations for the polynomial degrees on the subelements. This is illustrated in the figure on the right. For example, if a triangular or quadrilateral element is subdivided into four subelements where the polynomial degrees are allowed to vary by at most two, then this yields 3^4 = 81 refinement candidates. Analogously, splitting a hexahedron into eight subelements and varying their polynomial degrees by at most two yields 3^8 = 6,561 refinement candidates. Clearly, standard FEM error estimates providing one constant number per element are not enough to guide automatic hp-adaptivity.

Higher-order shape functions

In standard FEM one only works with shape functions associated with grid
vertices. In contrast to that,
in the hp-FEM one moreover regards edge functions, face functions,
and bubble functions. The following images show these functions
:
Note: all these functions are defined in the entire element interior!

Open source hp-FEM codes