A new anisotropic mesh adaptation method based upon hierarchical a posteriori error estimates

TL;DR

This paper introduces a novel anisotropic mesh adaptation method for finite element solutions of elliptic PDEs, utilizing hierarchical a posteriori error estimates to generate efficient meshes with minimal computational cost.

Contribution

The paper presents a new anisotropic mesh adaptation strategy based on hierarchical a posteriori error estimates, improving efficiency and reliability over existing methods.

Findings

Few Gauss-Seidel iterations suffice for accurate error approximation.

The method outperforms strategies using local error estimators or recovered Hessians.

Numerical examples demonstrate effectiveness in heat conduction problems with material jumps.

Abstract

A new anisotropic mesh adaptation strategy for finite element solution of elliptic differential equations is presented. It generates anisotropic adaptive meshes as quasi-uniform ones in some metric space, with the metric tensor being computed based on hierarchical a posteriori error estimates. A global hierarchical error estimate is employed in this study to obtain reliable directional information of the solution. Instead of solving the global error problem exactly, which is costly in general, we solve it iteratively using the symmetric Gauss--Seidel (GS) method. Numerical results show that a few GS iterations are sufficient for obtaining a reasonably good approximation to the error for use in anisotropic mesh adaptation. The new method is compared with several strategies using local error estimators or recovered Hessians. Numerical results are presented for a selection of test examples and a mathematical model for heat conduction in a thermal battery with large orthotropic jumps in the material coefficients.

Tables

Table 1. Table 1 : Heat conduction in a thermal battery: material coefficients and boundary conditions.

Region $k$	$D_x^k$	$D_y^k$	$f^k$
1	25	25	0
2	7	0.8	1
3	5	0.0001	1
4	0.2	0.2	0
5	0.05	0.05	0

Table 2. (a) Material coefficients.

Region $k$	$D_x^k$	$D_y^k$	$f^k$
1	25	25	0
2	7	0.8	1
3	5	0.0001	1
4	0.2	0.2	0
5	0.05	0.05	0

Table 3. (b) Boundary conditions.

Boundary $i$	${\alpha }^i$	$g^i$
1	0	0
2	1	3
3	2	2
4	3	0

Equations

$$(P)\{\begin{array}{l}\text{Find }u\in V\text{ such that}\\ a(u,v)=F(v),\forall v\in V\end{array}$$

$$(P_h)\{\begin{array}{l}\text{Find }{u}_h\in V_h\text{ such that}\\ a(u_h,v_h)=F(v_h),\forall v_h\in V_h.\end{array}$$

$$({\mathcal{T}}_h^{(0)},V_h^{(0)})\to u_h^{(0)}\to M_h^{(0)}\to ({\mathcal{T}}_h^{(1)},V_h^{(1)})\to u_h^{(1)}\to M_h^{(1)}\to \cdots$$

$$\Vert R_h{u}_h-u\Vert \le \beta \Vert u_h-u\Vert ,$$

$$\Vert u-u_h\Vert \le \frac{1}{1-\beta }\Vert R_h{u}_h-u_h\Vert .$$

$$I_h{R}_h{v}_h=v_h\forall v_h\in V_h$$

$$\Vert u-u_h\Vert \le \frac{1}{1-\beta }\Vert R_h{u}_h-u_h\Vert =\frac{1}{1-\beta }\Vert R_h{u}_h-I_h{R}_h{u}_h\Vert .$$

$$\Vert v-I_{h}v\Vert \le C\cdot \mathcal{E}({\mathcal{T}}_h,v),$$

$$\Vert u-u_h\Vert \le \frac{C}{1-\beta }\mathcal{E}({\mathcal{T}}_h,R_h{u}_h).$$

$${\bar{V}}_h=V_h\oplus W_h,$$

$${\bar{I}}_h{w}_h=0,\forall w_h\in W_h$$

$$R_h{u}_h=u_h+z_h$$

$$\Vert R_h{u}_h-{\bar{I}}_h{R}_h{u}_h\Vert =\Vert u_h+z_h-u_h\Vert =\Vert z_h\Vert =\Vert z_h-{\bar{I}}_h{z}_h\Vert .$$

$$\Vert u-u_h\Vert \le \frac{1}{1-\beta }\Vert z_h-{\bar{I}}_h{z}_h\Vert \le \frac{C}{1-\beta }\mathcal{E}({\mathcal{T}}_h,z_h).$$

$$a(e_h,v)$$

$$(E_h)\{\begin{array}{l}\text{Find }{z}_h\in W_h\text{ such that}\\ a(z_h,w_h)=F(w_h)-a(u_h,w_h)\forall w_h\in W_h.\end{array}$$

$$\frac{1}{d}\mathrm{tr}\left({(F_K^{\prime })}^T{M}_K{F}_K^{\prime }\right)=\det {\left({(F_K^{\prime })}^T{M}_K{F}_K^{\prime }\right)}^{\frac{1}{d}},\forall K\in {\mathcal{T}}_h$$

$$M_K=\frac{1}{|K|}{\int }_{K}M(x)dx.$$

$$|K|\sqrt{\det (M_K)}=\frac{{\sigma }_h}{N},\forall K\in {\mathcal{T}}_h$$

$${\sigma }_h=\sum _{K\in {\mathcal{T}}_h}|K|\sqrt{\det (M_K)}.$$

$${\int }_K\sqrt{\det (M_K)}dx=|K|\sqrt{\det (M_K)}.$$

$$\Vert v-I_{h}v{\Vert }_{L^q(K)}^q$$

$$\Vert v-I_{h}v{\Vert }_{L^q(\Omega )}^q\le C\sum _{K\in {\mathcal{T}}_h}|K|{\left(\mathrm{tr}\left({(F_K^{\prime })}^T|H_K|F_K^{\prime }\right)\right)}^q.$$

$$\mathcal{E}({\mathcal{T}}_h,v)=\sum _{K\in {\mathcal{T}}_h}|K|{\left(\mathrm{tr}\left({(F_K^{\prime })}^T|H_K|F_K^{\prime }\right)\right)}^q.$$

$$\mathcal{E}({\mathcal{T}}_h,v)$$

$$\ge d^q\sum _{K\in {\mathcal{T}}_h}|K|{\left(\det \left({(F_K^{\prime })}^T|H_K|F_K^{\prime }\right)\right)}^{\frac{q}{d}}$$

$$=d^q\sum _{K\in {\mathcal{T}}_h}|K{|}^{\frac{d+2q}{d}}{\det (|H_K|)}^{\frac{q}{d}}$$

$$=d^q\sum _{K\in {\mathcal{T}}_h}{\left(|K|{\det (|H_K|)}^{\frac{q}{d+2q}}\right)}^{\frac{d+2q}{d}}$$

$$\ge d^q{N}^{-\frac{2q}{d}}{\left(\sum _{K\in {\mathcal{T}}_h}|K|{\det (|H_K|)}^{\frac{q}{d+2q}}\right)}^{\frac{d+2q}{d}},$$

$$d^q{N}^{-\frac{2q}{d}}{\left({\int }_{\Omega }{\det (|H|)}^{\frac{q}{d+2q}}dx\right)}^{\frac{d+2q}{d}},$$

$$\frac{1}{d}\mathrm{tr}\left({(F_K^{\prime })}^T|H_K|F_K^{\prime }\right)={\det \left({(F_K^{\prime })}^T|H_K|F_K^{\prime }\right)}^{\frac{1}{d}},\forall K\in {\mathcal{T}}_h.$$

$$M_K={\theta }_K|H_K|$$

Full text

This is a preprint of a contibution published by Elsevier Inc. in J. Comput. Phys., 229(6) (2010), pp. 2179–2198.

© 2009. Licensed under CC-BY-NC-ND 4.0 (https://creativecommons.org/licenses/by-nc-nd/4.0).

The final version is available online at https://dx.doi.org/10.1016/j.jcp.2009.11.029.

††journal: Journal of Computational Physics

1 Introduction

Anisotropic mesh adaptation has proved to be a useful tool in numerical solution of partial differential equations (PDEs). This is especially true when problems arising from science and engineering have distinct anisotropic features. The ability to adapt the size, shape, and orientation of mesh elements according to certain quantities of interest can significantly improve the accuracy of the solution and enhance the computational efficiency.

Criteria for an optimal anisotropic triangular mesh were already given by D’Azevedo [1] and Simpson [2] in the early nineties of the last century. A number of algorithms for automatic construction of such meshes have since been developed.

A common approach for generating an anisotropic mesh is based on generation of a quasi-uniform mesh in some metric space. A key component of the approach is the determination of an appropriate metric often based on some type of error estimates. Unfortunately, classic isotropic error estimates do not suit this purpose well because they generally do not take the directional effect of the error or solution derivatives into consideration. This explains the recent interest in anisotropic error estimation; for example, see anisotropic interpolation error estimates by Formaggia and Perotto [3], Huang [4], and Huang and Sun [5]. Such error estimates for numerical solution of PDEs can be found, among others, in works by Apel [6], Kunert [7], Formaggia and Perotto [8], and Picasso [9].

It is worth pointing out that most existing anisotropic error estimates are a priori, requiring information of the exact solution of either the underlying problem or its adjoint, which is typically unavailable in a numerical simulation. A widely-used approach of avoiding this difficulty in practical computation is to replace the information by one recovered from the obtained numerical approximation. A number of recovery techniques can be used for this purpose, such as the gradient recovery technique by Zienkiewicz and Zhu [10, 11] and the technique based on the variational formulation by Dolejší [12]. Zhang and Naga [13] have recently proposed a new algorithm to reconstruct the gradient (which can also be used to reconstruct the Hessian) by fitting a quadratic polynomial to the nodal function values and subsequently differentiating it. It has been shown by Zhang and Naga [13] and by Vallet et al. [14] that the latter is robust and works best among several recovery techniques. Generally speaking, recovery methods work well when exact nodal function values are used but may lose some accuracy when applied to finite element approximations on non-uniform meshes. Nevertheless, the optimality of mesh adaptation based on those recovered approximations can still be proven under suitable conditions, see Vassilevski and Lipnikov [15]. More recently, conditions for asymptotically exact gradient and convergent Hessian recovery from a hierarchical basis error estimator have been given by Ovall [16]. His result is based on superconvergence results by Bank and Xu [17, 18], which require that the mesh be uniform or almost uniform.

The objective of this paper is to study the use of a posteriori error estimates in anisotropic mesh adaptation. Although a posteriori error estimates are frequently used for mesh adaptation, especially for refinement strategies and recently also for construction of equidistributing meshes for numerical solution of two-point boundary value problems by He and Huang [19] as well as in connection with the moving finite element method by Lang et al. [20], up to now only few methods for their use in anisotropic mesh adaptation have been published. For example, Cao et al. [21] studied two a posteriori error estimation strategies for computing scalar monitor functions for use in adaptive mesh movement; Apel et al. [22] investigated a number of a posteriori strategies for computing error gradients used for directional refinement; and Agouzal et al. [23] proposed a new method for computing tensor metrics provided that an edge-based a posteriori error estimate is given. Moreover, Dobrowolski et al. [24] have pointed out that error estimation based on solving local error problems can be inaccurate on anisotropic meshes. This shortcoming of local error estimates can be explained by the fact that they generally do not contain enough directional information of the solution, which is global in nature, and that their accuracy and effectiveness are sensitive to the aspect ratio of elements, which can be large for anisotropic meshes. We thus choose to develop our approach based on error estimation by means of globally defined error problem. To enhance the computational efficiency, we employ an iterative method to obtain a cost-efficient approximation to the solution of the corresponding global linear system. Numerical results show that a few symmetric Gauß–Seidel iterations are sufficient for this purpose. This is not surprising since the approximation is used only in mesh generation and it is often unnecessary to compute the mesh to a very high accuracy as for the solution of the underlying differential equation. Numerical experiments also show that the new approach is comparable in accuracy and efficiency to methods using Hessian recovery. We also test it with a more challenging example: a heat conduction problem for a thermal battery with large and orthotropic jumps in the material coefficients.111A Sandia National Laboratories benchmark problem.

The outline of the paper is as follows. In section 2, the new framework of using a posteriori hierarchical error estimates for anisotropic mesh adaptation in finite element approximation is described. In section 3, the optimal metric tensor based on the interpolation error is developed. Several implementation issues are addressed in section 4. Numerical results obtained with the new approach and with Hessian recovery-based methods are presented in section 5 for a selection of test examples. Numerical results for the heat conduction problem are given in section 6. Finally, section 7 contains conclusions and comments.

2 Model problem and adaptive finite element approximation

In this section, we describe a new framework of using a posteriori hierarchical error estimates for anisotropic mesh adaptation in finite element approximation.

2.1 Model problem and finite element approximation

Consider the boundary value problem of a second-order elliptic differential equation. Assume that the corresponding variational problem is given by

[TABLE]

where $V$ is an appropriate Hilbert space of functions over a domain $\Omega\in\mathbb{R}^{2}$, $a(\cdot,\cdot)$ is a bilinear form defined on $V\times V$, and $F(\cdot)$ is a continuous linear functional on $V$. The finite element approximation $u_{h}$ of $u$ is the solution of the corresponding variational problem on a finite dimensional subspace $V_{h}\subset V$, i.e.,

[TABLE]

If the bilinear form $a(\cdot,\cdot)$ is coercive and continuous on $V$, both variational problems $(P)$ and $(P_{h})$ have unique solutions. The finite dimensional subspace $V_{h}$ is often chosen as a space of piecewise polynomials associated with a given mesh, say $\mathcal{T}_{h}$, on $\Omega$. The variational problem $(P_{h})$ results in a system of $\dim(V_{h})$ linear algebraic equations.

2.2 Adaptive linear finite element solution

In this work we consider a linear finite element method, where $V$ is taken as $H^{1}(\Omega)$ and $V_{h}$ is the space of continuous, piecewise linear functions over $\mathcal{T}_{h}$.

Let $\mathcal{T}_{h}^{(i)}$ ($i=0,1,\dotso$) be an affine family of simplicial meshes on $\Omega$ and $V_{h}^{(i)}$ the corresponding space of continuous, piecewise linear functions. The adaptive solution is the result of an iterative process described as follows.

We start with an initial mesh $\mathcal{T}_{h}^{(0)}$. On every mesh $\mathcal{T}_{h}^{(i)}$ we solve the variational problem $(P_{h})$ with $V_{h}^{(i)}$ and use the obtained approximation $u_{h}^{(i)}$ to compute a new adaptive mesh for the next iteration step. The new mesh $\mathcal{T}_{h}^{(i+1)}$ is generated as an almost uniform mesh in a metric space with a metric tensor $M_{h}^{(i)}$ defined in terms of $u_{h}^{(i)}$. This yields the sequence

[TABLE]

The process is repeated until a good adaptation is achieved. An example of such adaptive meshes is shown in Fig. 1.

Typically, the metric tensor $M_{h}$ depends on the Hessian of the exact solution of the underlying problem [3, 25], which is often unavailable in practical computation. The common approach to avoid this difficulty is to recover an approximate Hessian from the computed solution. We consider here an alternative approach, which uses an a posteriori error estimator for defining and computing $M_{h}$.

2.3 Mesh adaptation based on a posteriori error estimates

Let $R_{h}$ be a reconstruction operator applied to the numerical approximation $u_{h}$. It can be either a recovery process, a smoothing operator, or an operator connected to an a posteriori error estimate. We assume that the reconstruction $R_{h}u_{h}$ is better than $u_{h}$ in the sense that

[TABLE]

for a given norm $\|\cdot\|$, where $0\leq\beta<1$ is a constant.

From the triangle inequality we immediately have

[TABLE]

If the reconstruction $R_{h}$ has the property

[TABLE]

for some interpolation operator $I_{h}$, we can bound the finite element approximation error by the (explicitly computable) interpolation error of the reconstructed function $R_{h}u_{h}$, viz.,

[TABLE]

Moreover, from the interpolation theory we know that the interpolation error for a given function $v$ can be bounded by a term depending on the triangulation $\mathcal{T}_{h}$ and derivatives of $v$, i.e.,

[TABLE]

where $C$ is a constant independent of $\mathcal{T}_{h}$ and $v$. Therefore, we can rewrite 4 as

[TABLE]

In other words, up to a constant, the solution error is bounded by the interpolation error of $R_{h}u_{h}$.

2.4 Hierarchical basis

One possibility to achieve the property 3 is to use the hierarchical decomposition of the finite element space. Let

[TABLE]

where $W_{h}$ is a hierarchical extension of $V_{h}$ to $\bar{V}_{h}$. Each $\bar{v}_{h}\in\bar{V}_{h}$ has a unique representation $\bar{v}_{h}=v_{h}+w_{h}$ with $v_{h}\in V_{h}$ and $w_{h}\in W_{h}$. If an interpolation operator, $\bar{I}_{h}:\bar{V}_{h}\mapsto V_{h}$, can be defined such that

[TABLE]

and if we define $R_{h}$ through

[TABLE]

for some $z_{h}\in W_{h}$, then we shall have the property 3 and the estimate 6. Moreover,

[TABLE]

Consequently, we can estimate the finite element approximation error by evaluating the interpolation error of $z_{h}$, i.e.,

[TABLE]

In the context of a posteriori error estimates, $z_{h}$ is typically taken as a hierarchical basis error estimator.

2.5 A posteriori error estimate based on hierarchical basis

The computation of the error estimator is based on a general framework, details on which can be found among others in the work of Bank and Smith [26] or Deuflhard et al. [27]. The approach is briefly explained as follows.

Let $u_{h}\in V_{h}$ be a linear finite element solution of the variational problem $(P_{h})$ and let $\bar{V}_{h}=V_{h}\oplus W_{h}$, where $W_{h}$ is the linear span of the edge bubble functions. Obviously, $\bar{V}_{h}$ is a subspace of piecewise quadratic functions. Moreover, we can define $\bar{I}_{h}$ as the vertex-based, piecewise linear Lagrange interpolation. This interpolation satisfies 7 since the edge bubble functions vanish at vertices.

Let $e_{h}=u-u_{h}$ be the error of the finite element solution $u_{h}$. Then for all $v\in V$ we have

[TABLE]

The error estimate $z_{h}$ is then defined as the solution of the approximate error problem

[TABLE]

The estimate $z_{h}$ can be viewed as a projection of the true error onto the subspace $W_{h}$. Note that this definition of the error estimate is global and its solution can be costly. Several solution methods will be discussed in section 4.

Once $z_{h}$ is determined, the reconstruction $R_{h}u_{h}$ is derived from 8. Then, if assumption 1 holds, the finite element approximation error can be controlled by minimizing the interpolation error of $z_{h}$, i.e., the right-hand side in 9. In this paper, we construct optimal metric tensors with respect to interpolation error estimates $\mathcal{E}(\mathcal{T}_{h},z_{h})$ for the $L^{2}$ norm. We assume that the reconstruction $R_{h}u_{h}=u_{h}+z_{h}$, where $z_{h}$ is computed from ($E_{h}$), gives a better approximation to $u$ than $u_{h}$, i.e., $\beta<1$ in 1.

3 Metric tensor based on linear interpolation error estimate

3.1 Equidistribution and alignment

Let $\Omega$ be a polyhedral domain in $\mathbb{R}^{d}$ and let $\mathcal{T}_{h}$ be a simplicial triangulation on $\Omega$. For every element $K\in\mathcal{T}_{h}$, there exists an affine invertible mapping $F_{K}\colon\hat{K}\to K$ such that $K=F_{K}(\hat{K})$, where $\hat{K}$ is the reference element. We assume that $\hat{K}$ has been chosen to be equilateral and have a unitary volume. We denote the Jacobian matrix of $F_{K}$ by $F_{K}^{\prime}$ and the number of elements in $\mathcal{T}_{h}$ by $N$.

As mentioned before, we consider an adaptive anisotropic mesh as a uniform mesh in the metric specified by a metric tensor $M$. Such a mesh is referred hereafter to as an $M$-uniform mesh. It can be characterized by shape-orientation and size requirements on mesh elements; see [28].

Alignment condition (i.e., shape-orientation requirement). The elements of an $M$-uniform mesh $\mathcal{T}_{h}$ are equilateral in the metric specified by $M$. This can be expressed as

[TABLE]

where $M_{K}$ is the average of $M$ on element $K$,i.e.,

[TABLE]

The left-hand side term of equality 11 is equal to the arithmetic-mean of the eigenvalues of matrix ${(F_{K}^{\prime})}^{T}M_{K}F_{K}^{\prime}$ while the right-hand side term is equal to their geometric-mean. The arithmetic-mean geometric-mean inequality implies that 11 holds if and only if the eigenvalues of matrix ${(F_{K}^{\prime})}^{T}M_{K}F_{K}^{\prime}$ are all equal. Element $K$ is equilateral in the metric $M_{K}$ when it satisfies 11.

Equidistribution condition (i.e., size requirement). The elements of an $M$-uniform mesh have an equal volume in the metric $M$, i.e.,

[TABLE]

where

[TABLE]

Note that the left-hand side of 12 is equal to the volume of element $K$ in metric $M_{K}$, i.e.,

[TABLE]

3.2 Anisotropic interpolation error bound for piecewise quadratic functions

Elementwise anisotropic interpolation error estimates are developed in [3, 8, 5]. Here, we follow the theory in [5]. Consider the piecewise linear Lagrange interpolation ($k=1$) of a piecewise quadratic function $v$ on an arbitrary mesh $\mathcal{T}_{h}$. The elementwise interpolation error measured in the $L^{q}$ norm ($q\geq 1$) is given by

[TABLE]

where $H_{K}$ is the Hessian of $v$ on the element $K$, $|H_{K}|=\sqrt{H_{K}^{T}H_{K}}$, $C$ is a constant independent of $\mathcal{T}_{h}$ and $v$, and $\mbox{tr}(\cdot)$ denotes the trace of a matrix. Note that $H_{K}$ is constant on $K$ since by assumption $v$ is quadratic on the element. Summing over all elements of $\mathcal{T}_{h}$ provides an upper bound for the global interpolation error

[TABLE]

One may notice that we have used $L^{q}$ norm for the error. As we shall see later (cf. 20), an optimal global error bound in this norm can be obtained for the non-regularized case. In principle, the same procedure also works for other norms or semi-norms particularly the $H^{1}$ semi-norm. However, it is unclear that the interpolation error bounds obtained in [5] for other norms will lead to an optimal global bound for $M$-uniform meshes.

From this, we can set $\mathcal{E}(\mathcal{T}_{h},v)$ in 5 to

[TABLE]

It has a lower bound as

[TABLE]

where we have used the arithmetic-mean geometric-mean inequality in 15 (recalling the trace and determinant of a matrix are equal to the sum and product of its eigenvalues, respectively) and Hölder’s inequality in 16. If $\max_{K\in\mathcal{T}_{h}}\operatorname{diam}(K)\rightarrow 0$, where $\operatorname{diam}(K)$ denotes the diameter of $K$, we see that the asymptotic lower bound on $\mathcal{E}(\mathcal{T}_{h},v)$ is

[TABLE]

which is invariant for all meshes of the same number of elements $N$. Thus, a mesh on which $\mathcal{E}(\mathcal{T}_{h},v)$ attains a lower bound 16 can be considered to be an asymptotically optimal mesh.

3.3 Optimal metric

The optimal metric $M$ is defined such that the interpolation error bound $\mathcal{E}(\mathcal{T}_{h},v)$ defined in 14 attains its lower bound 16 on $M$-uniform meshes of $N$ elements associated with $M$.

We first notice that equality in 15 holds if the $M$-uniform mesh satisfies

[TABLE]

Comparing this with the alignment condition 11, a property satisfied by the $M$-uniform mesh, suggests that $M$ be defined as

[TABLE]

with some scalar function $\theta_{K}$.

Next we notice that equality in 16 holds if the mesh satisfies

[TABLE]

Comparing this to the equidistribution condition 12, another property satisfied by the $M$-uniform mesh, yields

[TABLE]

This condition can be used for determining $\theta_{K}$. Thus, we obtain the optimal metric tensor as

[TABLE]

The interpolation error bound 14 attains its lower bound 16 on any $M$-uniform mesh associated with this metric tensor. From 13 we obtain

[TABLE]

for any $M$-uniform mesh associated with the metric tensor 18. Bound 20 has been obtained in [5] for $q=2$ and obtained and shown to be optimal in [29] for general $q\geq 1$.

The metric tensor defined by 18 is not necessarily positive definite since both $|H_{K}|$ and $\det(|H_{K}|)$ can vanish locally. To avoid this difficulty, the error bound is regularized with a positive parameter $\alpha_{h}$, i.e.,

[TABLE]

Using the same procedure as above, by minimizing the above (regularized) error bound we obtain the optimal metric tensor as

[TABLE]

The regularization parameter plays a role of controlling the intensity of mesh adaptation. Indeed, as $\alpha_{h}\to\infty$, $M_{K}\to I$ and a uniform mesh results. On the other hand, as $\alpha_{h}\to 0$, the mesh adaptation is increasingly reliant on $|H_{K}|$. To balance between these situations, we follow [5] and define $\alpha_{h}$ through the algebraic equation

[TABLE]

or equivalently

[TABLE]

where the factor $2^{\max\left\{1,\frac{dq}{d+2q}\right\}}$ has been used so that lower and upper bounds can be obtained for $\alpha_{h}$; see 24 and its derivation below. With this definition, about half of the mesh elements are concentrated in regions where $\det(M_{K})$ is large [5]. Moreover, $M_{K}$ is invariant under a scaling transformation of $v$.

Equation 23 has a unique solution since its left-hand side is monotonically decreasing with $\alpha_{h}$ increasing (assuming that $|H_{K}|$ is not all zero for all elements of $\mathcal{T}_{h}$), and tends to $+\infty$ (which is greater than the right-hand side) as $\alpha_{h}\to 0$ and $|\Omega|$ (which is less than the right-hand side) as $\alpha_{h}\to\infty$. Moreover, it can be solved using a simple iteration scheme such as the bisection method. Furthermore, lower and upper bounds on $\alpha_{h}$ can be obtained,

[TABLE]

Indeed, from 23 we have

[TABLE]

which leads to the right inequality of 24. On the other hand,

[TABLE]

which gives the left inequality of 24.

The interpolation error bound for a corresponding $M$-uniform mesh can be obtained as follows. From 21 and using the equidistribution and alignment conditions we have

[TABLE]

For $\alpha_{h}$ defined in 23, $\sigma_{h}=2^{\max\left\{1,\frac{dq}{d+2q}\right\}}|\Omega|$. Combining this with 21 we obtain

[TABLE]

In our computation we use the mesh generation software bamg (bidimensional anisotropic mesh generator developed by F. Hecht [30]) to generate new adaptive meshes for a given metric tensor $M$. Note that bamg requires that the metric tensor be further normalized such that all elements have a unitary volume in the metric. Thus, in actual computation we use a normalized metric tensor

[TABLE]

where $N$ is the desired number of mesh elements and

[TABLE]

It is remarked that the metric tensor can also be normalized using a prescribed error level; see [25].

4 Computation of the metric tensor and anisotropic meshes

We discuss here some implementation issues for two-dimensional problems.

The computation typically starts with a coarse regular Delaunay mesh of the domain and a desired number of mesh elements, $N$. For a given triangular mesh $\mathcal{T}_{h}^{(i)}$ at step $i$, we compute the numerical approximation $u_{h}^{(i)}$ with a standard linear finite element method. Based on $u_{h}^{(i)}$ and $\mathcal{T}_{h}^{(i)}$, we then compute $z_{h}^{(i)}$ as an approximation to the solution of the approximate error problem $(E_{h})$. Once $z_{h}^{(i)}$ has been obtained, it is straightforward to compute its elementwise Hessian and define the new metric tensor $M^{(i)}$ according to 22,

[TABLE]

where the error is measured in the $L^{2}$-norm, i.e., $q=2$. A new mesh is generated with bamg according to the metric tensor $\mathcal{M}^{(i)}={\left(\sigma_{h}^{(i)}/N\right)}^{-1}M^{(i)}$. The process is repeated until a good adaptation (see discussion below) is achieved.

4.1 Mesh quality measure

In order to characterize the mesh adaptation quality and to define an appropriate stopping criterion for the mesh adaptation process, we introduce the alignment and equidistribution quality measures [4]

[TABLE]

and

[TABLE]

which characterize how closely the mesh satisfies the alignment and equidistribution conditions 11 and 12, respectively.

Using $M^{(i)}$, $Q_{ali}^{(i)}$, and $Q_{eq}^{(i)}$, the estimate 21 can be reformulated as

[TABLE]

where

[TABLE]

is the overall mesh quality measure and takes into account both the shape and the size of elements. Since $Q_{ali}$ and $Q_{eq}$ appear in $Q_{mesh}$ as a product, their effects are not independent but compensate for each other. As a consequence, the mesh can have a good overall quality when small elements are shaped worse than large elements or well-aligned elements are worse shaped than worse aligned elements. Note that $Q_{ali}$, $Q_{eq}$, $Q_{mesh}$ $\geq 1$; and $Q_{ali}=Q_{eq}=Q_{mesh}=1$ if and only if the underlying mesh is $M$-uniform (cf. 25).

In the following numerical tests, the mesh adaptation process is stopped when

[TABLE]

where $\varepsilon_{mesh}$ is a tolerances chosen as $\varepsilon_{mesh}=0.1$ in our computation.

4.2 Computation of the error estimator

A key component of the procedure is to find the solution $z_{h}$ of problem $(E_{h})$. Note that $(E_{h})$ is a global problem and finding its exact solution can be as costly as for computing a quadratic finite element approximation to the original PDE problem. Three approaches are considered here for solving or approximating $(E_{h})$.

Edge-based error estimator. The expense of the error estimation can be significantly reduced, if the bilinear form $a$ in $(E_{h})$ is replaced by an approximation $\tilde{a}$ that allows a more efficient solution of the resulting linear system. A very efficient approach in two dimensions is to reduce the original problem to a series of local error problems which are defined over two elements sharing a common edge and can be solved efficiently. The approach is equivalent to the application of one Jacobi’s iteration (starting from zero) to the linear system resulting from the global error problem, i.e to the replacement of the stiffness matrix resulting from $(E_{h})$ by its diagonal. This approach has been successfully used in finite element computations [27, 31, 20]. Moreover, it has been shown [27] that such an error estimator is spectrally nearly equivalent to the original one under suitable conditions.

Despite its success in isotropic mesh adaptation, the approach does not seem to work well for anisotropic mesh adaptation. This may be explained by the fact that estimators based on local error problems generally depend on the aspect ratio of elements and can become inaccurate when the aspect ratio is large, a case that is often true for anisotropic meshes. Moreover, such estimators may not contain enough directional information of the solution which is global in nature and essential to the success of anisotropic mesh adaptation.

Node-based error estimator. This approach is similar to the edge-based error estimator, with the error estimator being obtained by solving a series of local error problems defined on node patches with homogeneous Dirichlet boundary conditions.

Inexact solution of the full error problem. In this approach the full error problem is kept but only an approximation to its exact solution is sought and used for the computation of the metric tensor. In our experiments, a few symmetric Gauß–Seidel iterations are employed to obtain such an approximation. In the following computation, Gauß–Seidel iterations are repeated until the relative difference of the old and the new approximations is under a given tolerance GS-RTOL.

It is noted that globally defined error estimators have the advantages that they are often independent of element aspect ratio and contain more directional information of the solution. Moreover, it is known [24] that the full hierarchical basis error estimator is efficient and reliable for anisotropic meshes.

Numerical comparison among these approaches is given in the next section.

5 Numerical examples

In this section, we present some numerical results for a selection of two-dimensional problems with an anisotropic behaviour. We first compare different approaches in solving the error problem $(E_{h})$ and then the new method with some common Hessian recovery methods. At the end of the section, we give further examples to demonstrate the ability of the method to generate appropriate anisotropic meshes.

Convergence is illustrated by plotting the finite element solution error against the number of elements. We use the $L^{2}$-norm for the error because the monitor function $M_{K}$ is optimized for this norm. For the inexact solution of the full error problem, $\text{\emph{GS-RTOL}}=0.01$ is chosen as a relative tolerance for the iterative Gauss–Seidel approximation.

5.1 A first example

Consider the boundary value problem

[TABLE]

with $\Omega=(0,1)\times(0,1)$. The right-hand side $f$ and the Dirichlet boundary conditions are chosen such that the exact solution is given by

[TABLE]

The solution exhibits a strong anisotropic behaviour and describes the interaction between a boundary layer along the $x$-axis and a shock wave along the line $y=x-0.5$. A solution plot is given in Fig. 1(a).

Reduced vs. full error estimators. As mentioned in section 4, on anisotropic meshes, there can be a significant difference in accuracy between estimators obtained by solving localized error problems and those obtained by means of a globally defined error problem. In our first test, we investigate the influence of the three error estimators described in the previous section on mesh adaptivity.

Results for the error of the adaptive solution against the number of elements are presented in Fig. 2.

As expected, the full error estimator works best, leading to a smaller error than those obtained with local error estimators. The node-based error estimator works better than the edge-based error estimator, mainly because it involves more elements and, in this sense, is more global.

The same observation can also be made from Fig. 3, where adaptive meshes obtained with the error estimators are shown. For these mesh example, the desired number of mesh elements $N$ in the normalized metric tensor given by 26 has been set to $600$. All methods produce correct mesh concentrations, although mesh alignment and orientation are different. In the mesh controlled by the full error estimator elements near the boundary layer and the shock wave are very thin, have a large aspect ratio222Aspect ratio is longest edge divided by shortest altitude. An equilateral triangle has an aspect ratio of $\sqrt{3}/2\approx 0.9$. of up to $46.9$, and are properly aligned with the fronts of the shock wave and the boundary layer (Fig. 3(c)). On the other hand, the elements of meshes controlled by reduced error estimators have rather moderate aspect ratios of $12.8$ and $14.3$ and are less anisotropic (Figs. 3(a) and 3(b)).

The accuracy of the corresponding finite element solutions is different, too. The mesh controlled by the full error estimator leads to a solution error $\|e\|_{L^{2}}=$1.6\text{\times}{10}^{-3}, less then one half of $\|e\|_{L^{2}}=$3.7\text{\times}{10}^{-3}, the error obtained using the node-based error estimator, and about one third of $\|e\|_{L^{2}}=$5.0\text{\times}{10}^{-3}$$, the error achieved with the edge-based error estimator.

These results are in good agreement with the comments made in section 4 that the full error estimator will do a better job than reduced ones for anisotropic mesh adaptation. Reduced error estimators are able to capture the distribution of the magnitude of the true error and yield a good mesh concentration. However, they fail to produce proper mesh alignment, i.e., they does not contain enough information for proper shape and orientation adaptation.

Effects of the number of Gauß–Seidel iterations. We now investigate how many iterations are sufficient for obtaining a valuable approximation to the error equation. Figure 4(a) presents results for different iteration numbers to compute the full error estimator. As one can see, a few iterations are sufficient for obtaining an approximation good enough for mesh adaptation. The convergence lines are very close to each other. The exact solution of the error problem leads to a smaller error, but the difference is hardly visible. Three steps of the symmetric Gauß–Seidel method produce an almost optimal mesh for this example.

Comparison to Hessian recovery methods. Two Hessian recovery methods are considered for comparison purpose.

Quadratic least squares fitting. This method was recently developed by Zhang and Naga [13] and proved to be robust and reliable. It computes a local quadratic fitting to function values or their approximations at some neighboring points and obtains a Hessian approximation by differentiating the polynomial twice.

Variational formulation. This approach recovers the Hessian, which does not exists in the classical sense for piecewise linear functions, by means of a variational formulation [12]. Precisely, let $\phi_{i}\in V_{h}$ be the piecewise linear basis function at node $(x_{i},y_{i})$. Then the nodal approximation to the second-order derivative $u_{xx}$ of a function $u$ at $(x_{i},y_{i})$ is defined as

[TABLE]

The same approach is used to approximate $u_{xy}$ and $u_{yy}$.

Figure 4(b) shows the error against the number of elements for each method. For comparison purpose, results obtained using the analytical Hessian are also included. All methods provide almost the same results. Particularly, the method based on the global estimator with three Gauß–Seidel iterations is comparable to the recovery-based methods.

It is worth noting that although the quadratic least squares fitting is generally more accurate and robust than the variational method, both produce basically the same adaptive mesh. This seems to confirm the conjecture that highly accurate Hessian recovery is not necessary for good mesh adaptation.

5.2 Further examples

We consider two boundary value problems in the form 27 with now the right-hand side $f$ and the Dirichlet boundary condition being chosen such that the exact solution is given by the following functions:

[TABLE]

The first function represents a shock wave along the line $y=1.25-x$ while the second models a boundary layer near the coordinate axes.

We compare the error for finite element solutions obtained with the global error estimator and the quadratic least squares Hessian recovery. Results for the quasi-uniform (regular Delaunay) mesh and the edge-based error estimator are also given. Figures 5 and 6 show the results.

As in section 5.1, we can see that mesh adaptation significantly reduces the finite element error compared to a quasi-uniform mesh having the same number of elements. The mesh based on the edge-based error estimator provides a good mesh concentration and is clearly better than a quasi-uniform one, but it is almost isotropic and inferior to a mesh obtained with the use of the full error estimator. Again, one can observe that the elements of the meshes obtained by means of the full error estimator and the quadratic least squares fitting are properly aligned with the shock wave and the boundary layers. Thus, the new method produces results comparable to those obtained with recovery-based methods.

5.3 Discontinuous gradients

Next, we consider problems whose solution has a discontinuous gradient along a certain interface in the domain. This situation arises in elliptic problems with discontinuous coefficients in the diffusion term such as heat conduction problems with jumps in material coefficients. Difficulties when using gradient recovery methods for such problems were already pointed out in [32], and this is true for the Hessian recovery as well: if the numerical approximation is accurate enough, we should expect a discontinuity in its gradient and its Hessian. Since most Hessian recovery methods employ some sort of averaging over a certain region, they can be very inaccurate near discontinuities. This issue can readily be observed in the following simple example.

Let $\Omega=(0,1)\times(0,1)$. Consider the boundary value problem

[TABLE]

where

[TABLE]

and the Dirichlet boundary condition is chosen such that the exact solution is given by

[TABLE]

The solution has a gradient jump of magnitude $\alpha$ across the line $x=0.5$, but is continuous on $\Omega$ and linear in each of the subdomains. We take $\alpha=6$ in our computation.

We first consider the situation where the mesh does not contain the information of the interface. In this situation at least part of the interface does not consist of edges. In order to match the sharp bend in the solution along the interface, the adaptive mesh should exhibit a strong concentration of elements around $x=0.5$ oriented along the interface. In this test, the quadratic least squares and the global hierarchical basis error estimator both succeed in providing an appropriate mesh adaptation and, again, deliver comparable results (Fig. 7(b)).

The situation is different if the interface is present in the mesh. In this case, the analytical solution $u$ belongs to the corresponding finite element space and, consequently, the numerical approximation computed by means of the linear finite element method is exact (Fig. 7(a), right). Hence, no adaptation is required and the proper mesh should be a uniform mesh. Now, consider the mesh adaptation using the quadratic least squares Hessian recovery. Because of the sharp bend in the solution, the recovered Hessian should be very large, $\approx\mathcal{O}(1/h)$, near $x=0.5$, but zero elsewhere, because the solution is linear in each of the subdomains. This should lead to an excessive over-adaptation near the interface. On the other hand, we expect no adaptation for the hierarchical basis error estimator in this case because the numerical solution is exact and, consequently, the error estimator is zero everywhere in $\Omega$. A quasi-uniform mesh should result. Figure 7(c) presents mesh examples. We see that the adaptation by means of the Hessian recovery (left) leads to a strong element concentration along the interface line, as predicted, whereas the mesh based on the hierarchical error estimator (right) is almost uniform.

We also expect a similar behaviour of these methods for general problems exhibiting gradient jumps or similar discontinuities along internal interfaces. Thus, for such problems, it can be of advantage to use the a posteriori error estimator for effective mesh adaptation because of the more efficient employment of given degrees of freedom.

6 Heat conduction in a thermal battery

In this section, we consider heat conduction in a thermal battery with large orthotropic jumps in the material coefficients. The mathematical model considered here is taken from [32, 33] and described by

[TABLE]

where $\Omega=(0,8.4)\times(0,24)$ and

[TABLE]

The data for each material $k$ and for each of the four sides $i$ of the boundary starting with the left-hand side boundary and ordering them clockwise are given in Table 1.

The analytical solution for this problem is unavailable. The geometry and the contour and surface plots of a finite element approximation are given in Fig. 8.

We compare the quadratic least squares Hessian recovery and the full error estimator. For this example we found that three steps of the symmetrical Gauß–Seidel method were not sufficient for a full mesh adaptation and increased the number to seven, which proved to be enough to achieve at least a comparable error estimate as the one obtained with quadratic least squares Hessian recovery.

Fig. 9 shows global error estimates (obtained by solving exactly the approximate error problem $(E_{h})$) for finite element solutions on adaptive meshes controlled by the full error estimate or Hessian recovery and having all or no predefined interface edges. (The interface consists of edges when a mesh has all predefined interface edges.)

Typical adaptive meshes with predefined interface edges for both methods are shown in Figure 10.

The results are in good agreement with those in section 5.3. When the interface edges are not present in the mesh, both methods provide similar results. On the other hand, when the mesh contains all the information of the interface, the quadratic least squares Hessian recovery produces a mesh with strong element concentration near all internal interfaces (Fig. 10(a)), whereas the full error estimator leads to a mesh (cf. Fig. 10(b)) that has higher element concentration in the corners of the regions, has a proper element orientation near the interfaces between the regions $2$ and $3$, and is almost uniform in regions where the solution is nearly linear (cf. Fig. 8(c) for the surface plot of a computed solution).

Meshes without predefined interface edges are quite similar to those in the example with discontinuous gradients (sections 5.3 and 7(b)). The interfaces are recognized by the both methods and the obtained adaptive meshes are dense near the interfaces.

Once again, the numerical results for this example show that a recovery method can lead to over-concentration of elements. The new method, on the other hand, produces only necessary concentration and is also able to catch the directional information of the solution required for proper element alignment. This example also demonstrates that the new method can be successfully used for problems with jumping coefficients and strong anisotropic features.

7 Conclusions and comments

In the previous sections, we have presented a mesh adaptation method based on hierarchical basis error estimates and shown that anisotropic mesh adaptation can be successfully controlled by a posteriori error estimators. Numerical results have shown that the new method is fully comparable in accuracy with commonly used Hessian recovery-based methods and can be more efficient for some examples by producing only necessary element concentration.

A key idea in the new approach is the use of the full hierarchical error estimator for reliable directional information of the solution. To avoid the expensive exact solution of the global error problem, we employed only a few steps of the symmetric Gauß–Seidel iteration for the efficient solution of the resulting linear system. Numerical results have shown that this is sufficient for obtaining an approximation to the error good enough for the purpose of mesh adaptation.

Acknowledgments

The work was partially supported by the German Research Foundation (DFG) under grants SFB568/3 and SPP1276 (MetStroem) and by the National Science Foundation (USA) under grants DMS-0410545 and DMS-0712935.

The authors are grateful to the anonymous referees for their valuable comments.