Multi-level hp-adaptivity and explicit error estimation

Recently, a multi-level hp-version of the finite element method (FEM) was proposed to ease the difficulties of treating hanging nodes, while providing full hp-approximation capabilities. In the original paper, the refinement procedure made use of a-priori knowledge of the solution. However, adaptive procedures can produce discretizations which are more effective than an intuitive choice of element sizes h and polynomial degree distributions p. This is particularly prominent when a-priori knowledge of the solution is only vague or unavailable. The present contribution demonstrates that multi-level hp-adaptive schemes can be efficiently driven by an explicit a-posteriori error estimator. To this end, we adopt the classical residual-based error estimator. The main insight here is that its extension to multi-level hp-FEM is possible by considering the refined-most overlay elements as integration domains. We demonstrate on several two- and three-dimensional examples that exponential convergence rates can be obtained.

The finite element method (FEM) has been shown to produce particularly efficient approximations when both refinement in element size h and polynomial degree p are considered (hp-FEM). In this way, exponential convergence can be attained also in presence of singularities [1,2,3]. However, the implementation of hp-FEM is challenging. This is due to the fact that local refinements produce hanging nodes, edges and faces. The associated degrees of freedoms destroy the required \( C^0 \)-continuity of the approximations [4, 5] if not treated appropriately. To this end, it is a common approach to constrain the concerned degrees of freedom. However, it turns out that constraining becomes extremely complex, especially in cases where n-irregular meshes have to be handled in three dimensions. Thus, many codes only allow for 1-irregular meshes [3,4,5].

To overcome these difficulties, multi-level approaches have been developed, e.g., [6,7,8,9,10,11,12]. For a in-depth review, see, e.g., [13]. Along this line of research, the recently introduced multi-level hp-method allows to remove the hanging nodes by construction [14], while allowing the approximation capabilities comparable to the classical hp-FEM [15]. The underlying idea concerns the refinement procedure in which the refinement is not performed by replacement of elements, but by hierarchically overlaying a finer mesh. This technique is extensively explained in [13, 14] and briefly recaptured in “Multi-level FEM” section.

In previous contributions, the multi-level hp-refinement procedure made use of a-priori knowledge of the solution, e.g. singularities given by re-entrant corners. However, adaptive procedures can automatically produce discretizations that are more effective than intuitive choices of element sizes h and polynomial degree distributions p. This is particularly prominent when a-priori knowledge of the solution is only vague or unavailable.

The present contribution demonstrates that multi-level hp-adaptive schemes can be driven by an explicit a-posteriori error estimator as well. This kind of estimator was introduced, analyzed theoretically and used for adaptive computations in [16,17,18,19,20]. Successively, error estimation and adaptivity have become to be broadly researched and proven to be robust, see, e.g., [21,22,23,24,25,26,27]. Furthermore a first extension to high-order was given in, e.g, [18, 28, 29]. For a comprehensive survey, see, e.g., [30,31,32,33,34].

In the rest of the paper, “Estimated error-based hp-adaptivity” section briefly introduces the multi-level hp-method and the explicit error estimator. It is then demonstrated that the simple extension of viewing the refined-most elements as integration domains suffices for this classical method to drive multi-level adaptivity. Next, “Implementational aspects” section discusses some important implementational aspects. The article then proceeds with various numerical example in two- and three-dimensions in “Numerical examples” section.

Model problem

We consider the Poisson’s Equation on a d-dimensional bounded domain \(\Omega \subset \mathbb {R}^d\). We assume the boundary \(\Gamma \) of  \(\Omega \) to be Lipschitz and consisting of two disjoint parts \(\Gamma _D\) and \(\Gamma _N\), representing the Dirichlet- and Neumann-Boundary, respectively, where \(\Gamma _D\) is assumed to be closed. We denote by \( \pmb {n} \) the vector normal to the boundary \( \Gamma _N \). The strong form of the considered problem reads

Instead, its weak formulation reads

where

Here, \( H^1(\Omega ) \subset L_2(\Omega ) \) denotes the classical Sobolev space of functions in \( L_2(\Omega ) \) with weak-derivatives up to order one.

Let \(n \in \mathbb {N}, \; n>0\), and \({\mathcal {T}}(\Omega )=\left\{ T_e \right\} _{e=1}^n\) be a finite partition of \(\Omega \) into quadrilaterals, if \( d=2 \), or hexahedra, if \( d=3 \). \({\mathcal {T}}\) is assumed to be regular, i.e., without hanging nodes. Considering polynomial high-order approximations, we associate to each element T a polynomial degree \( p_{{{T}}}\) and a diameter \( h_{{{T}}}\). In the sequel, the polynomial degree distribution is denoted as \( p = \left\{ p_{{{T}}}\right\} _{T\in {\mathcal {T}}(\Omega )} \). In the FEM, it is natural to consider each element \( T_e \) to be the image of a standard reference element \( \hat{T} \) under a mapping \( F_{T_e}: \hat{T} \rightarrow T_e \) [34]. For simplicity, assume \( F_{T_e} \) to be invertible. Furthermore, let \( {{S}_{{\mathcal {T}}(\Omega ), p}}\) denote the corresponding finite element solution space [34].

Here, \( \mathbb {P}_{p_{{{T}}}}(\hat{T}) \) is the space of polynomials of degrees at most \( p_{{{T}}}\) on \(\hat{T}\).

Multi-level FEM

In the standard FEM, the discretization is based on a mesh \({\mathcal {T}}(\Omega )=\lbrace T_e \rbrace _{e=1}^n \). The multi-level FEM generalizes this framework by allowing for an overlay of meshes, as depicted in Fig. 1a.

Example of one-dimensional multi-level mesh and basis functions. a Multi-level mesh. b Basis functions for analysis

Full size image

Starting from a base mesh \({\mathcal {T}}^0(\Omega )=\lbrace T_{e^0} \rbrace _{e^0=1}^{n^0} \), additional meshes can be superimposed on domains of interest. Namely, one element \( T_{a} \in {\mathcal {T}}^0(\Omega ) \) can be refined by an overlay mesh \({\mathcal {T}}^1(T_{a})=\lbrace T_{a, e^1} \rbrace _{e^1=1}^{n^1}\). In the following, \( T_{a, e^1} \) will be called the sub-elements of \( T_{a} \), while \( T_{a} \) is the parent element of \( T_{a, e^1} \). Any element \( T_{a, b} \in {\mathcal {T}}^1(T_a) \) can be further refined by superimposing an additional mesh \( {\mathcal {T}}^2(T_{a, b}) \). For example, \( T_1\in {\mathcal {T}}^0 \) in Fig. 1a is overlaid by \( {\mathcal {T}}^1 = \lbrace T_{1,1}, T_{1,2} \rbrace \), while \( T_{1,1} \) is overlaid by \( {\mathcal {T}}^2=\lbrace T_{1,1,1}, T_{1,1,2}\rbrace \). This procedure of mesh-overlaying can be carried out an arbitrary number of times for each element of each overlay mesh. Finally, a polynomial degree has to be assigned to every element of each level to define its local basis. In this manner, a multi-level hierarchical structure of elements in defined.

The refinement-by-overlay procedure requires some precautions in order to ensure linear independence of the global shape functions and \( C^0 \)-continuity of the numerical solution. See [14] for details. For example, in one-dimensional domains \( C^0 \)-continuity can be guaranteed by requiring the local shape functions to vanish at the nodes of the underlying levels (c.f. Fig. 1a). Linear independence can be guaranteed by allowing high-order modes only on the leaf elements, i.e. the elements without any further refinement (c.f. Fig. 1a). Other polynomial degree distributions are possible, see [14] for details. The described multi-level mesh defines a set of global basis functions defined on the same domain \( \Omega \). For example, Fig. 1b shows the basis functions on \( \Omega =[0,1] \) from Fig. 1a collapsed to one level. This global basis can be used for classic analysis and a element-local basis point-of-view will be considered later in “Error estimator for multi-level FEM” and “Implementational aspects” sections.

Explicit error estimator

Many different techniques in error estimation can be found in the literature. A comprehensive survey is given in e.g.  [30,31,32,33,34]. For the paper at hand, we follow the strategy of a-posteriori residual-based error estimators.

Let \( u_{\mathcal {T}}\in {{S}_{{\mathcal {T}}(\Omega ), p}}\) be a numerical approximation to \( u \in H^1_D \) computed by solving Eq. (2). According to [35, 36] estimates \( \eta \) to the analytical error \( e({ u_{{{\mathcal {T}}}} }) = u - u_{\mathcal {T}}\) in the energy norm \(\Vert \cdot \Vert _{{\mathcal {B}}}\) can then be computed as

for some constant C. The functions \( \mathcal {R}_{T}({ u_{{{\mathcal {T}}}} }) \), \( \mathcal {R}_{\partial {T}}({ u_{{{\mathcal {T}}}} }) \) denote the interior- and boundary- residuals of the element \( {T}\) with respect to the strong form (1):

where \( \pmb {n} \) is the normal to \( \partial {T}\) and \( \left[ \nabla { u_{{{\mathcal {T}}}} }\cdot \pmb {n} \right] \) is the inter-element flux jump. Specifically, \( \left[ \nabla { u_{{{\mathcal {T}}}} }\cdot \pmb {n} \right] = \lim _{t \rightarrow 0} \left( \nabla { u_{{{\mathcal {T}}}} }\cdot \pmb {n} \right) \left( \pmb {x} + t \pmb {n} \right) -\left( \nabla { u_{{{\mathcal {T}}}} }\cdot \pmb {n} \right) \left( \pmb {x} - t \pmb {n} \right) \quad \forall \pmb {x} \in \partial {T}\setminus \Gamma \). See for example [33, 37], among others.

In general, the constant C is unknown and this makes it difficult to use Eq. (4) alone for assessing the quality of the numerical approximation. However, the element error indicators \( \eta ^2_{{{{T}}}} \) can identify the elements accounting for the highest error contributions. As illustrated later in “Smoothness indicator for the multi-level hp-adaptivity” section, this will be used to drive adaptivity.

Note that in [35, 36] Neumann boundary conditions were not considered. In the present paper, a direct heuristic generalization is formulated and investigated numerically.

Note also that it was proven for even degrees, that an explicit error estimator can be constructed just in terms of interior residuals [31, 38]. Instead, for odd degrees an explicit error estimator can be constructed using only the inter-element jumps [31, 39].

Error estimator for multi-level FEM

The first step towards using the error estimator (3) together with the multi-level FEM is to define what an element is. The estimator (3) employs a sum over all elements \( {T}\) of a traditional mesh \( {\mathcal {T}}\). However, as described above, multi-level meshes define a hierarchical structure of overlaid elements. Therefore, it is necessary to translate appropriately the multi-level structure to conventional finite elements. In this context, we will refer to the set \( {\mathcal {T}}\) to be used in (3) as the “partition” of the multi-level mesh \( {{\mathcal {T}}_{\textit{ml}}}\).

In the standard derivation of the explicit error estimator for the model problem, Green’s theorem is applied to express the energy norm of the approximation error in terms of interior- and boundary-residuals of each element [35, 37, 40]. This requires at least \( C^2 \)-continuity of the element local shape functions. In case of multi-level meshes, the local shape functions of a partition element \( {T}\) are the functions of all the levels \( {\mathcal {T}}^k_\textit{ml} \) that are non-zero on \( {T}\). Thus, not all choices for the partition of a multi-level mesh are suitable. For example, choosing \( \lbrace [0, 0.5] \), \( [0.5, 1] \rbrace \) as partition of \( {{\mathcal {T}}_{\textit{ml}}}\) in Fig. 1a only gives \( C^0 \)-continuity at \( x=0.25\).

The coarsest partition of \( \Omega \) ensuring \( C^2 \)-continuity of the local shape functions is the set of leaf elements, i.e., the set of elements of any level that are not further refined by a superimposed mesh. This set is denoted by \( {\mathcal {T}}^\textit{leaf}_{\textit{ml}}\). Note that the leaf elements \( {\mathcal {T}}^\textit{leaf}_{\textit{ml}}\) do not overlap each other and their union covers the whole domain. For example, \( {\mathcal {T}}^\textit{leaf}_{\textit{ml}}= \lbrace [0, 0.25], [0.25, 0.5], [0.5, 1] \rbrace \) in Fig. 1a. Considering \( {\mathcal {T}}^\textit{leaf}_{\textit{ml}}\) as partition \( {\mathcal {T}}\), an explicit error estimator for the multi-level FEM can be defined as follows:

Note that the error indicator is associated to the domains defined by each leaf element and not divided into the different levels.

It is noteworthy that this definition of partition as the traditional “elements” of the hierarchical multi-level mesh also complies with the set of sufficient conditions for the convergence of a finite element discretization given in [41]. Namely, for second-order problems the shape function shall represent exactly all polynomials of order up to one (completeness), be \( C^1 \)-continuous within each element (smoothness) and be \( C^0 \)-continuous across element boundaries (continuity). Despite the fact that these are just sufficient conditions, most of the finite element bases satisfy these requirements. Therefore, the above definition of elements of a multi-level mesh reveals to be meaningful in a more general sense, as it satisfies the smoothness requirement.

Note that, while each overlay mesh \( {\mathcal {T}}^k \) is assumed to be regular, \( {\mathcal {T}}^\textit{leaf}_{\textit{ml}}\) allows for an arbitrary level of hanging nodes (n-irregular mesh). These irregularities do not need complicated constraining algorithms, as \( C^0\)-continuity of the numerical solution is ensured by construction [14].

Smoothness indicator for the multi-level hp-adaptivity

In the context of adaptive refinement procedures, the error estimator allows to identify the elements that account for the highest error. Various methods exist to decide whether these elements should be p- or h-refined, see [2] for a comprehensive overview. In the sequel, the smoothness indicator proposed in [42] is employed. Its underlying idea is that for regions with a non-smooth solution, h-refinement is more effective than p-refinement, while areas with a smooth solution can be better approximated by p-refinement.

For one-dimensional problems, the Legendre coefficients of the numerical approximation \(u_{\mathcal {T}|{{T}}}\) local to \( {T}=[-1, 1] \) are used to indicate its smoothness. Denoting the i-th Legendre polynomial by \( P_i \), the Legendre coefficients can be computed as \( \alpha _i = (2 i + 1)/2 \int _{-1}^{1} P_i(x) u_{\mathcal {T}|{{T}}}(x) \; dx \), for \(0\le i \le p_{{{T}}}\). The decay rate \( \tau _{{T}}\) of these coefficients is estimated by a best fit of \( | \alpha _i | = c e^{\tau _{{T}}i} \). The decision between h- or p-refinement is then carried out depending on a parameter \( C_{\textit{decay}} \in \mathbb {R}\) as follows. If \( \tau _{{T}}\le C_{\textit{decay}} \), then \( u_{\mathcal {T}|{{T}}}\) is considered to be smooth and a p-refinement is performed. Otherwise, \( u_{\mathcal {T}|{{T}}}\) is classified as non-smooth and a h-refinement is chosen.

For higher-dimensional problems, we follow the anisotropic approach proposed in [43, 44]. In a two-dimensional problem, for example, the Legendre expansion coefficients are computed by

for \( 0\le i \le p_{{{T}}x}\), \(0\le j \le p_{{{T}}y} \). Here, \( p_{{{T}}x} \) and \( p_{{{T}}y} \) denote the polynomial degree in the x- and y- direction, respectively. From this second-order tensor of coefficients a one-dimensional sequence of coefficients which represents one spatial direction is obtained by accumulating the coefficients of the other directions. For example, for the x-direction \( | \tilde{\alpha }_i | = \sum _{j=0}^{p_{{{T}}y}} | \alpha _{ij} | \sqrt{ \frac{2}{2 j +1} } \). This allows to compute the decay rate \( \tau _{{{T}}x} \) in x-direction by applying the one-dimensional technique laid out above to \( | \tilde{\alpha }_i | \). Analogously this is carried out for the y-direction. Once a smoothness indicator for each spatial direction is obtained, anisotropic refinements can be performed.

In the context of this work, only isotropic refinements are considered. In particular, p-refinements are performed if \( \tau _{{T}}=\max \lbrace \tau _{{{T}}x}, \tau _{{{T}}y} \rbrace \le C_{\textit{decay}} \). Otherwise, a multi-level h-refinement is carried out, as described below.

The error estimator (6) and smoothness indicator are then combined to drive adaptivity as laid out in Algorithm 1. Here, \( C_{\textit{error}} \in \left[ 0,1 \right] \) is a chosen parameter representing the refinement threshold of the normalized element error. Note that \(C_{\textit{error}} \ne C\). Moreover, we allow high-order shape functions just on leaf elements, while non-leaf element are always linear. The p-refinement can be performed by local degree elevation on any leaf. The multi-level h-refinement of an element \( \bar{T}\) of order p is performed by superimposition on \( \bar{T}\) by smaller elements, all of order p. Namely, the children element inherit the order of the parent. Element \( \bar{T}\) has then superimposed children elements and is not a leaf anymore, therefore its order is set to 1. We assume that for any leaf element \( {T}\) in the mesh \( {\mathcal {T}}^\textit{leaf}_{\textit{ml}}\) produced after the i-th step of Algorithm 1, the ratio of its diameter \( h_{T}\) to the diameter \( \rho _{T}\) of the largest ball inscribed into \( {T}\) is bounded independently of \( {T}\) and i.

In the following, we consider as multi-level h-refinement of an element \( \bar{T}\), the superimposition on \( \bar{T}\) by by \( 2^d \) elements obtained by spatial bi-section of \( \bar{T}\) in all the d directions.

The multi-level FEM is implemented by defining a sequence of nested reference-spaces, each defining local shape functions, as depicted in Fig. 2. This approach furnishes a structure suitable for many recursive algorithms. In this section we present an approach on how to compute the derivative of the numerical solution and the inter-element flux jump in such a recursive setting.

Example of nested reference-spaces

Full size image

Differentiation recursive formulas

The definition of nested reference-spaces involves the usage of a sequence of mappings \( u(x(\xi ^0(\xi ^1(\ldots \xi ^n)))) \) that have to be considered when evaluating derivatives. See Fig. 2. Here, \( x_i \) (\(i=1,\ldots , d\)) denotes the physical coordinates, while \( \xi ^k_i \) \((i=1, \ldots , d) \) denotes the coordinates in the reference space of level k. In order to compute the residuals (4) and (5), first- and second-order derivatives have to be computed. Using Einstein’s notation, a direct application of the chain rule yields

These formulas are recursive and based on quantities directly computable on the element of level n. Iterating the computation until level \( x=\xi ^{-1} \) produces the desired result.

However, in general the backward mapping \( \xi ^k( \xi ^{k-1} ) \) is not always available. Therefore, it is indispensable to express the differentiation by only using forward mappings. To this end, the chain-rule can be applied to \( \partial ^2 u / \partial \xi _s^k \partial \xi _t^k \), yielding

Here, at the \( k{\text {th}}\) level, almost all quantities are known with the exception of \( \partial ^2 \, u/\partial \xi ^{k-1}_i\, \partial \xi ^{k-1}_j \) for which the following tiny systems have to be solved:

In vector notation, Eqs. (11) and (12) read:

where

Note that the term \( \pmb \nabla ^{k-1} u \) in (14) can be computed locally to each level by (13). This completes the necessary ingredients to use (14) in the setting of a non-linear mapping \( x(\xi ^0) \).

Usually, the mappings \( \xi ^{k-1}( \xi ^{k} ) \) for \( k=1 \ldots n \) are linear transformations. In this case Eq. (14) can be simplified to

Evaluation of the flux jump

In this section, it is illustrated how the multi-level definition of the basis functions can be exploited in the computation of the inter-element flux jump (5). Consider Fig. 3 and recall the definition of the entities \( {\mathcal {T}}^\textit{leaf}_{\textit{ml}}\) regarded as partition of a multi-level mesh \( {{\mathcal {T}}_{\textit{ml}}}\), as explained in “Error estimator for multi-level FEM” section. First, observe that edges can contain edges of levels above. For example, the edge AC of \( T_{1,2,1} \) contains the edge AB of \( T_{1,2,1,1} \). We refer to AC as a parent-edge of AB, while AB is a sub-edge of AC. The edges without sub-edges are called leaf edges. Moreover, given an element \( {T}\in {\mathcal {T}}^k_\textit{ml} \) of an overlay mesh, edges shared by sub-elements of \( {T}\) are called internal to \( {T}\). These internal edges are identified by dashed lines in Fig. 3, whereas boundary edges are marked by solid lines. Finally, we denote as ancestor of an edge as the edge itself, if it does not posses a parent, or (recursively) the ancestor of its parent-edge.

Example of two-dimensional multi-level mesh

Full size image

Consider an element \( T_1 \in {\mathcal {T}}^0_\textit{ml} \) of level 0 and its multi-level structure of refinements, as in Fig. 3. The key observation to evaluate the flux is that the ancestor \( \gamma _\textit{anc} \) of each boundary edge \( \gamma \) is internal to some element \( \tilde{T} \in {\mathcal {T}}^k_\textit{ml} \) of level k. Therefore, \( \gamma _\textit{anc} \) is contained in the volume of \( \tilde{T} \). As explained in “Error estimator for multi-level FEM” section, basis functions are assumed to be at least \( C^1 \)-continuous in the interior of \( \tilde{T} \). Hence, the shape functions defined locally on \( \tilde{T} \) do not give any contribution to the flux jump. Analogously, this holds for all the parent-elements of \( \tilde{T} \) of level \( l<k \). Therefore, it is in general sufficient to consider the contribution to the inter-element flux given by the sub-elements of \( \tilde{T} \) belonging to the level \( l>k \). For example, to evaluate the flux residual across AB, it is sufficient to evaluate the shape functions defined in \( T_{1,1} \) and \( T_{1,1,2} \) for the left flux. While to evaluate the right flux it is enough to consider the contribution given by \( T_{1,2} \), \( T_{1,2,1} \) and \( T_{1,2,1,1} \). In particular, no shape function needs to be evaluated in \( T_1\).

A possible algorithm to compute the residual on edges internal to one element \( {T}\) is given in Algorithm 2. This procedure is applied recursively to the sub-elements of \( {T}\).

In this section, the adaptive scheme is tested on two- and three-dimensional problems presenting singularities or steep gradients. In the following, the absolute error in the energy norm is computed as \( \Vert e({ u_{{{\mathcal {T}}}} })\Vert ^2_{{\mathcal {B}}} = \Vert u\Vert ^2_{{\mathcal {B}}} - \Vert { u_{{{\mathcal {T}}}} }\Vert ^2_{{\mathcal {B}}} \). Here, \( \Vert { u_{{{\mathcal {T}}}} }\Vert ^2_{{\mathcal {B}}} \) is computed numerically, while \( \Vert u\Vert ^2_{{\mathcal {B}}} \) is computed analytically whenever possible, or numerically otherwise. In case an explicit form of u is not available, \( \Vert u\Vert ^2_{{\mathcal {B}}} \) is estimated by an extrapolation method, as in [45, section 4.2]. Moreover, the quality of the estimator is assessed by means of the effectivity index \( \theta = \eta / \Vert E({ u_{{{\mathcal {T}}}} })\Vert _{{\mathcal {B}}} \). In the following, the multi-level hp-adaptive procedure is shown to achieve exponential convergence rates. In particular, error bounds of the form

are obtained, where N is the number of degrees of freedom and \( \alpha , \beta , \phi \in \mathbb {R}\).

L-shaped domain (2D)

Consider the Poisson’s problem (1) on the L-shaped domain \( \Omega \) depicted in Fig. 4a where \( f=0 \) and g is expressed in terms of the polar coordinates r, \( \theta \) as

This problem has an analytical solution \( u=r^\frac{2}{3} \sin (\frac{2}{3}\theta ) \) with a singular gradient at \( r=0 \). See Fig. 4b, c.

L-shape domain: problem description and solution. a L-shaped domain with boundary conditions. b Solution. c Gradient magnitude

Full size image

The discussed adaptive procedure is used to solve the above problem using \( C_{\textit{decay}} =-1.75 \), \( C_{\textit{error}} =0.8 \). Starting from a coarse mesh consisting of three elements of order one, the mesh obtained after 85 multi-level hp-adaptive steps is illustrated in Fig. 5a, b. Note the geometric h-refinement toward the re-entrant corner and linearly graded polynomial degree distribution in Fig. 5c. This kind of mesh is known to be well suited to resolve vertex singularities. Figure 6a compares the analytical error \( \Vert e({ u_{{{\mathcal {T}}}} })\Vert _{{\mathcal {B}}} \) to the estimated error \( \eta \) demonstrating two important results. First, the multi-level hp-adaptive procedure yields approximations converging exponentially to the analytical solution. In particular, this is shown by choosing \( \phi =1/3 \) in Eq. (16), in accordance with the theory of [46] for singular problems. Second, the error estimator tracks closely the behavior of the analytical error. In particular, the efficiency index presents a converging behavior, as depicted in Fig. 6b. Finally, in Fig. 6c, the adaptive procedure is compared to uniform h-, manual hp- and p-refinements on a fixed geometric multi-level mesh. Such meshes are obtained by iteratively performing multi-level h-refinement of the elements closer to the singularity. The number of iterations is denoted as depth in the legend. The manual hp-discretization is taken from [15].

L-shape domain: discretization. a Ansatz order. b Ansatz order around singular corner. c Polynomial degree distribution on the line x \(\ge \) 0, y = 0

Full size image

L-shape domain: error estimates, effectivity index and convergence. a Error estimates and convergence rate. b Effectivity Index. c Convergence comparison

Full size image

Singular cube (3D)

To test adaptivity on three-dimensional singular solutions, the problem given by Eq. (1) is considered on \( \Omega =[0,1]^3 \) with \( \Gamma _N = \partial \Omega \), \( \Gamma _D = \left\{ \pmb {0} \right\} \). Moreover, \( f = \lambda (\lambda + 1) r^{\lambda -2} \) and \( g = \lambda r^{\lambda -2} \pmb {x} \cdot \pmb {n} \), where \( \lambda \in \mathbb {R}\), \( \lambda > 0 \), \( r(\pmb {x}) = \Vert \pmb {x}\Vert _2 \). The analytical solution \( u = r^\lambda \) has a gradient singular at \( \pmb {0} \) for \( \lambda < 2 \).

The above problem is solved with \( \lambda =2/3 \), \( C_{\textit{decay}} =-1.75 \), \( C_{\textit{error}} =0.8 \) starting from a single linear element. Figure 7a depicts the numerical solution and Fig. 7b, c illustrate the discretization obtained after 67 adaptive steps. Note the geometrical h-refinement toward the singularity and the linear-like distribution of the polynomial degree shown in Fig. 7d. Furthermore, Fig. 8a shows exponential convergence of the analytical error, along with the error estimates. For this problem which contains a vertex-singularity, exponential convergence is attained for \( \phi =1/4 \) of Eq. (16). This matches the theory given in [46]. Finally, the effectivity index presents a convergent behavior, as shown in Fig. 8b. Convergence compared to uniform h-refinement, uniform p-refinement on graded mesh and manual hp-refinement is shown in Fig. 8c. Here, the manual hp-discretization is obtained by constructing for each value \(1\le D\le 9 \) a multi-level mesh geometrically refined toward the singularity with geometric factor 0.5 and depth D and polynomial degree distribution equal to \( D-L+1 \), for each leaf-element of overlay-level \( 0\le L\le D \). This is analogous to the description given in [15].

Singular cube: solution and discretization. a Numerical solution gradient magnitude. b Ansatz order close to singular corner. c Ansatz order on side. d Polynomial degree distribution on the line x = 0, y = 0

Full size image

Singular cube: error estimates, effectivity index and convergence. a Error estimates and convergence rate. b Effectivity index. c Convergence comparison

Full size image

Shock problem: solution and discretization. a Numerical solution. b Numerical gradient magnitude. c Ansatz order

Full size image

Shock cube (3D)

In this section, a problem with steep concentrated gradients is introduced. Let \( r(\pmb {x}) = \Vert \pmb {x}-\pmb {x_0}\Vert _2 \) for \( \pmb {x_0} \in \mathbb {R}^3 \) and \( \alpha ,r_0 \in \mathbb {R} \), and consider the boundary value problem 1 on \( \Omega =[0,1]^3 \) with \( \Gamma _N = \partial \Omega \), \( \Gamma _D = \left\{ \pmb {0} \right\} \). Moreover, let

The analytical solution to the problem reads \( u = \arctan ( \alpha (r-r_0) ) - \arctan ( \alpha (0-r_0) ) \). Note that \( \nabla u \) is concentrated on the sphere of radius \( r_0 \) with a steepness controllable by \( \alpha \). The problem is solved by the adaptive procedure starting from a coarse mesh of \( 2\times 2\times 2 \) linear elements. Figure 9a shows the gradient of the numerical solution obtained after 60 multi-level hp-adaptive steps for \( \alpha = 80 \), \( r_0 = \sqrt{3} \), \( \pmb {x_0}= [-1/4,-1/4,-1/4]^T \), \(C_{\textit{decay}} =-1 \) and \( C_{\textit{error}} =0.8 \). The obtained discretization is strongly refined in both h and p toward the steep gradient. This is visualized in Fig. 9b, c. Similar discretizations were also obtained for the same problem in [47]. Figure 10a, c show exponential convergence, while Fig. 10b plots the converging behavior of the effectivity index. For smooth problems, a p-refinement on a fixed mesh would produce error bounds as in Eq. (16) with \( \phi =1/3 \) [46]. The first part of the plot in Fig. 10a does not depict a straight line, as in this phase the h-refinements are performed. In the asymptotic part, the grid is fixed and just p-refinements are employed. Here, a clear exponential behavior is shown.

Shock problem: error estimates, effectivity index and convergence. a Error estimates and convergence rate. b Effectivity index. c Convergence comparison

Full size image

Geometrically complex example (3D)

Finally, we consider a problem with a more complex geometry taken from an industrial application. The Field Assisted Sintering Technology (FAST) is a manufacturing process to produce dense artifacts from powder by sintering. To this end, copper powder is contained under pressure in a tool made of graphite and simultaneously heated by an electric field, see Fig. 11a and [48] for further details. This process was simulated numerically in [49] using electro-thermo-mechanically coupled fields. In the present paper, we only consider the Poisson’s model problem on the original geometry to test the multi-level adaptive procedure. Due to the symmetry of the tool, it is possible to simulate one eighth of the whole structure only, as depicted in Fig. 11b. Therefore, we consider the problem 1 where \(\Gamma _N\) is the top surface of the structure (contained in the xy-plane at z \(=\) 67.55), \(\Gamma _D\) is the bottom surface (contained in xy-plane at z \(=\) 0) and g and f are constant: \( g=1 \), \( f=0 \). The remaining faces are subject to homogeneous natural boundary conditions.

Geometrically complex example: problem description, solution and discretization. a Schematic sketch of the FAST process from [49]. b Geometry description in [mm]. c Solution gradient magnitude on initial mesh. d Solution gradient magnitude on refined mesh. e Ansatz order after 30 steps. f Ansatz order after 57 steps

Full size image

The initial mesh is obtained by means of the TUM.GeoFrame mesh generator [50]. This coarse mesh is composed of elements with a polynomial degree \( p=2 \) and a second-order polynomial mapping function describing the geometry, as presented in [51]. The magnitude of the solution gradient is visualized in Fig. 11c along with the mesh. Figure 11d shows the mesh and the solution after 30 adaptive steps with \( C_{\textit{error}} =0.75, C_{\textit{decay}} =-1.5 \). Figure 11e, f illustrate the polynomial degree distribution after 30 and after 57 steps, respectively. Interestingly, in this example the smoothness indicator happens to favor h-refinements over p-extensions in regions where the solution is smooth. This suggests that further improvements are likely to be possible by applying different smoothness indicators. See, e.g.,  [2] for a comparison of different approaches found in the literature.

For this problem, no analytical solution is available such that the reference strain energy is approximated by the extrapolation described in [45, section 4.2]. The energies used for the extrapolation are obtained by uniform p-extensions using the trunk space [45] and listed in Table 1. For this type of geometry with edge singularities, an optimal hp-adaptive scheme is expected to achieve exponential convergence of the form (16) with \( \phi =1/4 \)  [46]. This, however, requires anisotropic h-refinement which is still a subject of further research in the context of the multi-level FEM. An exponential decay of the error can, thus, not be expected. Nevertheless, the error estimator tracks closely the behavior of the error and the adaptive procedure drastically improves the convergence rate of the standard p-extension, as clearly shown in Fig. 12a–c.

Geometrically complex example: error estimates, effectivity index and convergence. a Error estimates and convergence rate. b Effectivity index. c Convergence comparison

Full size image

This work demonstrates that the multi-level hp-adaptive scheme can be driven by means of an explicit error estimator and provides reasoning as well as the necessary formulae for its implementation. In this context, the decision between h- or p-refinement is carried out according to the decay rate of the Legendre coefficients of the numerical solution.

The considered examples include singularities or concentrated steep gradients which the error estimator tracks closely; the effectivity index also shows a converging behavior. Moreover, exponential convergence rates of the multi-level hp-adaptive procedure are shown. However, while the smoothness indicator performs excellently in the benchmark with simple geometries, this was not observed for the example presenting with a complex geometry. Nevertheless, efficient discretizations were found automatically also in this example.

All authors have prepared the manuscript. All authors read and approved the final manuscript.

Acknowledgements

With the support of the Technische Universität München-Institute for Advanced Study, funded by the German Excellence Initiative and the European Union Seventh Framework Programme under Grant Agreement Number 291763. The authors also gratefully acknowledge the financial support of the German Research Foundation (DFG) under Grants RA 624/22-1 and RA 624/27-1, and Project SCHR 1244/4-1 of the priority programme SPP 1748.

Competing interests

The authors declare that they have no competing interests.

Authors and Affiliations

1. Chair for Computation in Engineering, Technische Universität München, Arcisstraße 21, 80333, Munich, Germany

   Davide D’Angella, Nils Zander, Stefan Kollmannsberger, Felix Frischmann & Ernst Rank

2. Institute for Advanced Study, Technische Universität München, Lichtenbergstraße 2 a, Garching, Germany

   Davide D’Angella, Ernst Rank & Alessandro Reali

3. Department of Mathematics, Universität Salzburg, Hellbrunnerstraße 34, Salzburg, Austria

   Andreas Schröder

4. Dipartimento di Ingegneria Civile e Architettura, Università degli Studi di Pavia, Via Ferrata 3, Pavia, Italy

   Alessandro Reali

Corresponding author

Correspondence to Davide D’Angella.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions