An hp-adaptive strategy for elliptic problems

TL;DR

This paper introduces an hp-adaptive method for elliptic problems that intelligently refines the mesh or polynomial degree based on error and smoothness estimates, achieving exponential convergence.

Contribution

It presents a novel hp-adaptive strategy utilizing refinement history and smoothness estimates for improved convergence in elliptic problems.

Findings

Exponential convergence achieved with the proposed strategy.

Refinement decisions are guided by a posteriori error and smoothness estimates.

Numerical experiments validate the effectiveness of the method.

Abstract

In this paper a new hp-adaptive strategy for elliptic problems based on refinement history is proposed, which chooses h-, p- or hp-refinement on individual elements according to a posteriori error estimate, as well as smoothness estimate of the solution obtained by comparing the actual and expected error reduction rate. Numerical experiments show that exponential convergence can be achieved with this strategy.

Tables

Table 1. Table 1: The final meshes (Examples 4.1)

	# elements	# DOF	Energy error
HP/PHG	3,772	246,046	1.01e-4
HP/MK	35,696	1,171,216	1.67e-4
HAFEM	1,663,068	2,263,137	3.57e-2

Table 2. Table 2: Final meshes (Example 4.2)

	# elements	# DOF	Energy error
HP/PHG	3,429	155,812	1.07e-5
HP/MK	163,204	1,158,279	1.20e-4
HAFEM	1,377,588	1,904,054	4.44e-4

Equations

$$-\Delta u=f\text{ on }\Omega ,u=g\text{ in }\partial \Omega .$$

$$a(u,v)=L(v),\forall v\in H_0^1(\Omega ),$$

$$a(u,v)={\int }_{\Omega }\nabla u\cdot \nabla vdx,L(v)={\int }_{\Omega }fvdx.$$

$$a(u_{hp},v_{hp})=L(v_{hp}),\forall v_{hp}\in V_{hp},$$

$$\Vert e\Vert \le C_1{N}_d^{-\frac{p}{d}},$$

$$\Vert e\Vert =C_1{N}_d^{-\frac{p}{d}}=C_2(\sum _{K\in T}{\eta }_K^2{)}^{\frac{1}{2}},$$

$$C_1(c_K{N}_d{)}^{\frac{-p}{d}}=C_2(\sum _{K^{\prime }\in T_1}{\eta }_{K^{\prime }}^2{)}^{\frac{1}{2}}=C_2(c_K{\lambda }_h^2\sum _{K\in T}{\eta }_K^2{)}^{\frac{1}{2}}.$$

$${\lambda }_h^2=\frac{1}{c_K}(\frac{1}{c_K}{)}^{\frac{2p}{d}}.$$

$${\lambda }_h=(\frac{1}{c_K}{)}^{\frac{p}{d}}.$$

$$\Vert e{\Vert }_{H^1(\Omega )}\le C\frac{h^{\mu }}{p^{(m-1)}}\Vert u{\Vert }_{H^m(\Omega )},$$

$$\Vert e{\Vert }_{H^1(\Omega )}=C\frac{h^{\mu }}{(p+1{)}^{(m-1)}}\Vert u{\Vert }_{H^m(\Omega )}=C(\frac{p}{p+1}{)}^{m-1}\frac{h^{\mu }}{p^{(m-1)}}\Vert u{\Vert }_{H^m(\Omega )}.$$

$${\lambda }_p=(\frac{p}{p+1}{)}^{m-1}.$$

$${\lambda }_p=(\frac{p}{p+1}{)}^{\frac{p}{2}}.$$

$${\lambda }_{hp}=(\frac{p}{p+1}{)}^{\frac{p}{2}}(\frac{1}{c_K}{)}^{\frac{p}{d}}.$$

$${\eta }_K\ge \alpha \underset{K^{\prime }\in T}{\max }{\eta }_{K^{\prime }}\iff K\text{ is marked for refinement},$$

$${\eta }_K^2\le {\lambda }_h^2{\eta }_{K_m}^2.$$

$${\eta }_K^2\le {\lambda }_p^2{\eta }_{K_m}^2.$$

$${\eta }_K^2\le {\lambda }_{hp}^2{\eta }_{K_m}^2.$$

$${\lambda }_h=(\frac{1}{2}{)}^{\frac{1}{3}{p}_K},$$

$${\lambda }_p=(\frac{p_K-1}{p_K}{)}^{\frac{1}{2}(p_K-1)},$$

$${\lambda }_{hp}=(\frac{p_K-1}{p_K}{)}^{\frac{1}{2}(p_K-1)}(\frac{1}{2}{)}^{\frac{1}{3}(p_K-1)}.$$

$${\eta }_K^2=\frac{h_K^2}{p_K^2}\Vert f_{p_K}+\Delta u_{hp}{\Vert }_{L^2(k)}^2+\sum _{f\subset \partial K\cap \Omega }\frac{h_f}{2p_f}\Vert [\frac{\partial u_{hp}}{\partial n_f}]{\Vert }_{L^2(f)}^2,$$

$$\|e\|\leq C\exp\bigl{(}-\gamma(N_{d})^{1/5}\bigr{)},$$

Full text

††journal: Journal of Computational and Applied Mathematics

1 Introduction

The adaptive finite element method (AFEM) is a widely used numerical method for solving partial differential equations. The $h$-version of AFEM modifies the size of the elements ($h$-refinement) while keeping the polynomial degrees fixed [14]. The $p$-version of AFEM adjusts the polynomial degrees in the elements ($p$-refinement) while keeping the size of the elements fixed. The $hp$-version of AFEM is more general, which consists of combining freely $h$-refinement and $p$-refinement. The $hp$-version of AFEM dates back to 1986, thanks to the pioneering work of Ivo Babuška et al. [9, 10, 11, 12, 13]. With $hp$ AFEM exponential convergence could be achieved if $h$-refinement and $p$-refinement are integrated properly [5, 8, 10].

One essential issue in the $hp$-adaptive finite element method is the design of refinement strategy, i.e., to decide which element should be refined and which kind of refinement should be performed. According to approximation theory, $p$-refinement should be performed on elements in which the solution to the partial differential equations is smooth and $h$-refinement should be performed on elements in which the solution is non-smooth [5]. Unfortunately, since the property of the solution is usually unknown, we need to estimate its smoothness using the computed numerical solution and other data. For this purpose many strategies have been proposed and developed. Owens et al. [15, 16] used a priori information of the computational domain and boundary data to determine the location of singularities of the solution, and performed $h$-refinement on elements which had singularities and $p$-refinement elsewhere. Oden et al. [18] introduced the so-called Texas-3-step strategy. Melenk et al. [5] and Heuveline et al. [7] proposed heuristic strategies which made use of the refinement history. Another class of strategies consisted of using error estimators obtained from solving local problems as indicators for guiding the refinement [8, 21]. For other strategies proposed and studied in the literature, we refer to [17, 13, 6, 8, 23, 19].

In this paper, we propose an $hp$-refinement strategy which is based on a posteriori error estimate and estimation of the smoothness of the solution using the reduction rates of the a posteriori error estimate in the refinement history. This strategy is mainly motivated by Melenk et al. [5] and Heuveline et al. [7], it removes the requirement of regular refinement and the dependence on mesh size $h$ in [5, 7], and can be applied to both two and three dimensional elliptic problems.

The layout of the paper is as follows. In §2, the model problem and notations are introduced. In §3, the $hp$-adaptive strategy is deduced in details. In §4, the efficiency of the new strategy is illustrated and compared to some other strategies through two numerical examples. In §5, some concluding remarks are given.

2 Model problem and notations

For a bounded Lipschitz domain $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, the following model problem is considered:

[TABLE]

where $f\in L^{2}(\Omega)$. The problem can be read in the weak form: find $u\in H_{0}^{1}(\Omega)$ such that

[TABLE]

where

[TABLE]

Our goal is to design an $hp$-finite element subspace $V_{hp}\subset H_{0}^{1}(\Omega)$ and to compute a numerical solution $u_{hp}\in V_{hp}$ such that

[TABLE]

and the error meets prescribed tolerance. Here for simplicity of description we will assume $g=0$. In this case $g\neq 0$, the problem can be easily converted to the case $g=0$ with a shift operator.

For the sake of convenience, some notations are introduced here. In the subsequent descriptions, we will denote by $u$ the exact solution of Problem (1), by $u_{hp}$ the numerical solution of the problem with respect to a triangulation $T$ and a finite element space $V_{hp}$ on $T$, and by $e=u-u_{hp}$ the error between the exact solution and the numerical solution. The energy norm, $\|\cdot\|$, is defined as $\|u\|=\sqrt{a(u,u)}$. In the adaptive process, $\varepsilon$ stands for the tolerance, which is the stop criterion, $\eta_{K}$ stands for the error indicator defined on element $K$, and $\eta=(\sum_{K\in T}\eta_{K}^{2})^{1/2}$ is the global error indicator. For a given element $K$, $p_{K}$ and $h_{K}$ denote the degree of the polynomial basis functions on $K$ and the diameter of $K$, respectively. When the element $K$ is divided (refined) into subelements, $c_{K}$ denotes the number of its children. Finally, $N_{d}$ is used to denote the total number of degrees of freedom in the mesh $T$.

3 An $hp$-adaptive strategy

In this section we give our $hp$-adaptive strategy. This strategy is based on the expected error reduction factors of $h$-, $p$-, or $hp$-refinement. The expected error reduction factors are calculated under the assumption that the numerical solution converges algebraically under $h$-refinement and exponentially under $p$-refinement. We will first deduce the expected error reduction factors for various refinement types, then describe the new $hp$-adaptive strategy in details.

First we deduce the expected error reduction factor $\lambda_{h}$ for $h$-refinement. We assume that the optimal convergence rate of the $h$-version of adaptive finite element method is algebraic, which can be written as [2, 3, 4],

[TABLE]

where $p$ is the degree of the piecewise polynomials. Suppose the fine mesh $T_{1}$ is obtained from uniform refinement of the mesh $T$ by dividing each element $K$ into $c_{K}$ subelements. Then the number of degrees of freedom on mesh $T_{1}$ is about $c_{K}N_{d}$.

Suppose we have an appropriate error indicator $\{\eta_{K}\mid K\in T\}$. We make the following hypotheses.

(H1) The error indicator is precise, i.e., there exist constants $C_{1}$ and $C_{2}$ such that,

[TABLE]

(H2) For any element $K$, the error indicators on all its children are equal.

Let $\lambda_{h}$ be the expected error reduction factor for $h$-refinement. By combining (H1) and (H2), we get the following relationship

[TABLE]

Comparing  (7) to  (6), we have

[TABLE]

To improve the efficiency, we use a slightly enlarged $\lambda_{h}$, which is given by

[TABLE]

Next we deduce the expected error reduction factor $\lambda_{p}$ for $p$-refinement. In $p$-refinement the mesh is fixed and the degree of the polynomials is adjusted. On a quasi-uniform mesh with uniform polynomial degree the following error estimation is expected [8, 9]

[TABLE]

where $h$ is the mesh size, $p$ the polynomial degree, $\mu=\min(p,m-1)$, $C$ a constant independent of $h$ and $p$, and $u\in H^{m}(\Omega)$. We make the following hypothesis.

(H3) $\|e\|_{H^{1}(\Omega)}=C\frac{h^{\mu}}{p^{(m-1)}}\|u\|_{H^{m}(\Omega)}$ and $p\geq(m-1)$.

When the degree $p$ is increased by one, by (H3) we have

[TABLE]

Thus the error reduction factor $\lambda_{p}$ is

[TABLE]

$m$ is a positive integer satisfying (H3). In this paper we set $m$ to $p/2+1$. Then we have

[TABLE]

Finally the expected error reduction factor $\lambda_{hp}$ for $hp$-refinement can readily be obtained by combining $\lambda_{h}$ and $\lambda_{p}$, which is given by

[TABLE]

As a widely accepted criterion in adaptive finite element methods, the error should be distributed asymptotically uniformly over all elements [5]. Therefore, elements with large error estimator should be marked for refinement. Here we employ the so-called maximum strategy, which can be described as follows

[TABLE]

where $\alpha\in(0,1)$ is a predetermined parameter.

Our $hp$-adaptive strategy is given below, which is motivated by Heuveline et al. [7] and Melenk et al. [5], using a similar framework. Here $h$-refinement means dividing the element into $c_{K}$ subelements, $p$-refinement means increasing the polynomial degree by 1.

1. Step 1:

Solve the problem on the current mesh $T$ with the current setting of polynomial orders and compute the error indicator $\{\eta_{K}\mid K\in T\}$ and the global error indicator $\eta$. The adaptive process is stopped if $\eta$ is less than or equal to $\varepsilon$ on the current mesh. 2. Step 2:

Mark elements for refinement using maximum strategy. 3. Step 3:

For each marked element $K$:

• (a)

If element $K$ is obtained by $h$-refinement of its parent element $K_{m}$, then check whether the following condition holds

[TABLE]

If yes then mark $K$ for $p$-refinement. Otherwise mark $K$ for $h$-refinement.

• (b)

If element $K$ is obtained by $p$-refinement of its parent element $K_{m}$, then check whether the following condition holds

[TABLE]

If yes then mark $K$ for $p$-refinement. Otherwise mark $K$ for $h$-refinement.

• (c)

If element $K$ is obtained by $hp$-refinement of its parent element $K_{m}$, then check whether the following condition holds

[TABLE]

If yes then mark $K$ for $p$-refinement. Otherwise mark $K$ for $h$-refinement.

• (d)

If element $K$ is not refined in the preceding adaptive step, then mark $K$ for $p$-refinement. 4. Step 4:

Perform $h$-, $p$- or $hp$-refinement as determined by Step 3.111When we perform $h$-refinement additional elements may be refined in order to maintain the conformity of the mesh. 5. Step 5:

Go to Step 1.

The underlying idea behind the above process is that because of the exponential convergence rate of $p$-refinement, it is preferred over $h$-refinement whenever the solution is smooth. If the expected error reduction factor is achieved in the previous refinement, then the solution is considered smooth and $p$-refinement is performed, otherwise $h$-refinement is performed.

Remark: the strategy proposed by Melenk et al. [5] was designed for two dimensional problems. Our strategy is suitable for both two and three dimensional problems and different error reduction factors are deduced. For the strategy proposed by Rannacher et al. [7], the error reduction factor depended on the size of elements. This dependency is removed in this paper.

4 Numerical results

In this section two examples are employed to illustrate the efficiency of the new $hp$-adaptive strategy. These examples are also computed using a traditional $h$-version adaptive finite element method and another existing $hp$-adaptive strategy for comparison.

We have implemented our new $hp$-adaptive strategy using the parallel adaptive finite element toolbox PHG [1]. The computations were performed on the cluster LSSC-III of the State Key Laboratory of Scientific and Engineering Computing, Chinese Academy of Sciences.

In these examples, since bisection refinement is used for $h$-refinement, we have $c_{K}=2$, thus the expected error reduction factors are given by

[TABLE]

[TABLE]

[TABLE]

The error indicator used here is the one introduced by Melenk et al. [5]. Though it was designed for two dimensional problems, it is also valid for three dimensional problems. This error indicator is given by

[TABLE]

where $f_{p_{K}}$ is the $L^{2}(K)$-projection of the function $f$ on the space of polynomials of degree $p_{K}-1$, $h_{f}$ denotes the diameter of the face $f$, $p_{f}=\max(p_{K_{1}},p_{K_{2}})$, where $K_{1}$ and $K_{2}$ are the two elements sharing the face $f$, and $[\cdot]$ denotes the jump of a function across the face $f$.

The parameter $\alpha$ in the maximum strategy is chosen as $0.5$. The linear systems of equations are solved by the PCG (Preconditioned Conjugate Gradient) method with a block Jacobi preconditioner. The initial meshes are generated using NETGEN [26] and the initial polynomial degrees on all elements are set to 2.

For three dimensional Poisson equation the optimal convergence rate is exponential and is expected to be [10]

[TABLE]

where $\gamma$ is a constant.

In the figures the logarithm of the energy error is plotted against $(N_{d})^{1/5}$, and three different strategies are compared. The first one is a traditional $h$-adaptive finite element method, denoted by “HAFEM”. The second one is the $hp$-adaptive strategy introduced in this paper, denoted by “HP/PHG”. The last one is the strategy of Melenk et al., denoted by “HP/MK”.

Example 4.1. In this example, the domain is an $L$-shaped domain given by $\Omega=(-1,1)^{3}\setminus(0,1]\times[-1,0)\times(-1,1)$, and the analytic solution is given by $u=\cos(2\pi x)\cos(2\pi y)\cos(2\pi z)$. The main difficulty in applying high order finite element methods to this problem is that the even and odd derivatives of the solution behave differently at each point in the domain, hence pure $p$-refinement may not improve the numerical solution [25]. The initial mesh is uniform with 144 elements.

The convergence histories of different strategies are shown in Figure 1 and statistics about the final meshes are shown in Table 1. We can observe that the two $hp$ strategies exhibit exponential convergence rate while the $h$-version converges algebraically. We can also observe that the HP/PHG strategy performs better than the HP/MK strategy.

Example 4.2. In this example, the computational domain is given by $\Omega=(-1,1)^{3}\setminus[0,1)^{3}$, and the analytic solution is given by $u={(x^{2}+y^{2}+z^{2})}^{\frac{1}{4}}$, whose gradient has a vertex singularity. The initial mesh is uniform with 172 elements.

The convergence histories and final meshes are shown in Figure 2 and Table 2 respectively. Again for this example, the $h$-version converges algebraically while the two $hp$-versions converge exponentially. Data in Table 2 shows that the performance of our strategy is much better than that of the HP/MK strategy.

5 Conclusion

A simple and easy to implement $hp$-adaptive strategy based on error reduction prediction is proposed. This strategy is suitable for two and three dimensional problems. The efficiency of the strategy is demonstrated through two numerical examples. Although the strategy is discussed with the Poisson equation in this paper, it is applicable to general elliptic problems. It also provides a general framework which can be easily extended to other problems.

Acknowledgments

This work is supported by the 973 Program under the grant 2011CB309703, by China NSF under the grants 11021101 and 11171334, by the 973 Program under the grant 2011CB309701, the China NSF under the grants 11101417 and by the National Magnetic Confinement Fusion Science Program under the grants 2011GB105003.