#### Abstract

A novel approach to anisotropic mesh-refinement of linear triangular elements in the finite element method is proposed. Using the method of solving a hierarchically refined dual problem in the a posteriori error analysis, indicators for edge refinement, rather than element refinement, can be obtained. A complete adaptive strategy is presented and illustrated by some numerical examples.

#### 1. Introduction

Designing the optimal finite element mesh for a certain problem is a key issue in order to make inexpensive and accurate computations. Standard refinement strategies rely on adapting the mesh size, that is, the element diameters, in different areas of the computational domain. The refinement of the elements is, typically, carried out with the goal of conserving or “improving” aspect ratios with some uniform shapes, confer to [1, 2]. However, for many problems, the need for resolution is very different in different directions. For instance, in case of boundary layers, much higher resolution is needed perpendicular to the boundary than in the parallel direction. The reason for this is the anisotropic nature of the error with respect to generation and transport. Traditionally, these effects are dealt with a priori when constructing the initial mesh. If the nature of the error is not known a priori, or in order to save time generating meshes, one wishes to adapt the mesh with respect to shape and orientation a posteriori together with the mesh-size adaptation. The important issue is then how to obtain indicators for such a refinement. Methods used for anisotropic mesh-refinement rely on estimation of error in different directions, for instance, in terms of interpolation errors or Hessian approximations, confer to [3–6]. For application to fluid mechanics, see [7–9]. In this paper, we propose a simple novel approach to relate regular, isotropic, error indicators to edges of triangular elements, whereby anisotropic refinement is achieved by subdividing edges with large error contributions.

Developments in establishing strategies for goal-oriented error computations rely heavily on the idea of solving a dual problem, confer to [10–12]. Much of current research is based on guaranteed bounds, confer to [13, 14]. However, resorting to approximate error estimation is often appealing, confer to [15], and —in the case of general nonlinear problems—necessary. The combination of goal-oriented error estimation and anisotropic mesh refinement is discussed by Richter [16].

Introducing an enhanced dual solution in terms of hierarchically added basis, as introduced in [17], standard element indicators can be replaced by edge indicators in a straight-forward manner. In particular, our approach enables anisotropic mesh-refinement without adding cost or complexity to the error estimation. Indicating the edges to be refined also gives the advantage that element refinement becomes a one-step local procedure, since the edges to refine are indicated for the mesh globally.

In this paper we first recall the general framework for a posteriori error computation for an abstract problem. Secondly, we discuss the relations between the a posteriori error estimation and error indicators related to the discretization. Next, a strategy for subdividing the mesh for given error indicators is suggested. Finally, we study the procedure on two test problems on a unit square, comparing anisotropic mesh-refinement with a (nearly) isotropic strategy. The Poisson problem illustrates the behavior for a problem of isotropic nature, while we use a reaction-diffusion problem with low conductivity to induce strong anisotropy along boundaries.

#### 2. A Posteriori Error Estimation

##### 2.1. Derivation of Exact Error Representation

For completeness, we recall the general strategy in establishing the error representation for an abstract problem as described by for example, Eriksson et al. [10] and Becker and Rannacher [11]. We thus consider the general nonlinear PDE in the space-time domain . Considering a standard Galerkin-scheme with homogeneous Dirichlet-conditions, the standard variational format of the problem is the following: Find such that where is the exact (weak) solution, is a semilinear form (linear in the second argument), is a linear functional, and is the appropriate Sobolev space (with the appropriate regularity and boundary/initial conditions).

In order to define the corresponding finite element solution, we introduce the approximation space . In the Finite Element Method, we construct as piecewise polynomials on elements , constituting the domain . We thus solve the problem. Find the approximate solution such that The difference in exact and approximate solution can now be defined as .

We will now introduce the appropriate secant form of . Upon using the Gateaux derivative with respect to the first argument of , we construct the secant form where , . The form is nonsymmetrical (in general) but bilinear in and . Using the fundamental theorem of calculus, we obtain the difference where we introduced the weak residual

From (2.2) we obtain the Galerkin orthogonality for general nonlinear equations: which will be used in the error analysis, confer to below.

We will now turn to the error analysis. In order to retain maximal generality in the formulation, we select the appropriate goal-oriented error measure that can be any (Gateaux-differentiable) functional that satisfies the condition for any function . Hence, may be a norm of a special case. Using the Gateaux-derivative of with respect to the first argument, we define the secant form of as follows: which is a linear functional in . We then obtain the relation

We are now in the position to introduce the variational format of the dual problem. Find such that

Having computed , we can use (2.8), (2.9), and (2.4) to formulate the * exact* error representation as follows:
for any . Here, the arbitrary function can be subtracted due to the Galerkin-orthogonality (2.6). Typically, we may choose as some projection of onto , for example, the nodal interpolant, the -projection, or the FE-interpolant (solution).

##### 2.2. Approximation of Error

In general, and depend on the solution , which is unknown. Hence, approximations are needed in the definition of the dual problem. For instance, the most straightforward assumption is to replace the secant forms with the corresponding tangent forms around the primal FE-solution, confer to [17]. Disregarding the nonlinear effects, we henceforth assume the approximate dual problem. Find such that where the approximations and were used.

*Remark 2.1. *In the case of a linear problem, the secant form is equal to the semilinear form of the variational format and therefore
which is linear in both arguments.

The exact error representation is given by (2.10). In order to approximate the error, we introduce the approximation of as . Following the method described in [17], we compute the enhanced dual solution in an enriched solution space .

First, we note that the dual problem can be solved using the same FE-mesh as for the primal problem, that is, is the solution of
Based on the solution of the dual problem on the primary FE-mesh, , the enhancement is solved globally or locally resulting in the enhancement . Substituting the enhanced dual solution, , into (2.10) and choosing , results in the following, approximate, error estimate:

#### 3. Error Indicators Contributions from Discretization

Based on the error estimation above, we wish to use an adaptive -refinement strategy. Thus, reliable indicators are needed, holding information of where and how the error is generated (in terms of lacking discretization). We will now turn to the special case of a 2D spatial (boundary value) problem, for which the FE-mesh is made up by linear triangular elements.

##### 3.1. Relating Error Indicators to Edge Refinement

Usually, error indicators are obtained through restricting the error estimator (2.10) or (2.14) to each * element *, whereby can be estimated or bounded by a sum of contributions from each element. Based on the error contribution, it is then decided upon whether or not to refine the specific element. Such a contribution scheme would appear as
where is the restriction of to . The contribution is related to the size of element , and in the adaptive procedure we hope to reduce by decreasing the local mesh-size.

We now propose, instead of computing the contributions of error from each element, to derive the error estimate as a sum of contributions from each element edge . Thus, we want to create a contribution scheme like where is the contribution to the chosen measure due to the interpolation error related to the edge , that is, the error that arises from the absence of further basis functions in the primal solution space .

First, we write the error contributions as a sum over the added basis functions. The enhancement space is spanned by the basis of hierarchical polynomials with corresponding coefficients , such that . Thus, (2.14) can be rewritten as since is linear in its second argument. The next step is to relate to . We propose to introduce coefficients such that In order for (3.2) and (3.3) to coincide, it is clear that the coefficients must satisfy the relation In addition to this, we need to construct the coefficients so that the error contribution is related to the edge in a reasonable way, that is, the error corresponding to must be reduced when edge is subdivided.

##### 3.2. Anisotropic Distribution of Error Indicators to Edges

We now turn to the case of using piecewise linear approximation for the primal problem (i.e., ) and enriching this with hierarchical basis functions, one edge bubble on each edge (full quadratic basis) and one internal bubble on each element (cubic parasitic term). The linear bases functions are depicted in Figure 1, while the hierarchical basis are depicted in Figure 2 (quadratic bases) and Figure 3 (cubic base), for one element.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

In view of the similarity of the quadratic basis function and the piecewise linear basis function that corresponds to a subdivision of that edge, confer to Figure 2, it is natural to relate the entire contribution to that edge. Thus the coefficients for all such that is an edge bubble are where indicates the edge on which the bubble lives.

Concerning the internal bubble function, that has its support inside an element, we distribute its error contribution over the surrounding edges based on fractions of the circumference. Thus the coefficients for all such that is an internal bubble function are where indicates edges surrounding the element on which the bubble lives and is the length of edge number .

*Remark 3.1. *In our experience, the use of a fully quadratic basis (one order higher in convergence) is in many cases sufficient to obtain sharp error estimate, confer to [17]. The rationale for using the, higher order, internal bubble functions as well as the quadratic edge functions is to obtain error indicators for * all* edges. If only edge-based hierarchical basis functions are added, no error contributions will be obtained for edges on Dirichlet boundaries, since means that the nodal values, , are zero as well. Thus, if higher resolution is needed along a Dirichlet boundary, it cannot be detected without adding internal functions in the hierarchical basis.

#### 4. Mesh Refinement and Adaptivity

##### 4.1. Subdivision of Triangular Elements in 2D

In the case of element indicators, one important issue is to control the * quality* of the refined mesh. This is in most cases done adopting strategies to either (i) conserve the initial shape of elements, confer to [1], or (ii) to thrive at making new elements as isotropic (uniform) as possible in each subdivision, confer to [2].

In the proposed method, the refinement of an element is governed by which edge/edges of the element is marked for division due to high error contribution. In this fashion, we wish to construct the optimal mesh, depending on the problem, with respect to both size, aspect ratio, and orientation of elements. If the nature of the error (generation and transport) is anisotropic, the algorithm will generate an anisotropic mesh, while for an isotropic behavior of the error, an isotropic mesh will develop during remeshing.

Turning again to the case of linear triangular elements, the refinement of an element is based on the edges to be subdivided, 0, 1, 2, or 3. In the case that 1 edge is indicated for subdivision, two new elements are created. If all three edges are marked for subdivision, an isotropic refinement of the old element into 4 new elements is performed. In the case that 2 of the edges are indicated, the old triangle is subdivided into three new triangles. Two possible subdivisions are possible, based on which of the two error indicators (on marked edges) is the highest. The edge corresponding to the highest error indicator is made vertex for all three triangles, thus enriching interpolation quality in the element primarily along that edge. The possible subdivisions of a triangle are illustrated in Figure 4.

**(a)**

**(b)**

**(c)**

**(d)**

Note that using the proposed method, each element is refined independently of its neighbors once the subdivided edges are indicated on a * global* level. Thus neither transitional elements nor iterative procedures need to be adopted.

##### 4.2. Refinement Strategy

The adaptive refinement is set to perform mesh-refinements until a certain tolerance is met, that is, until . When the error is larger than the preset tolerance, edges with indicators, larger than a certain value are refined. This value, the local tolerance, is chosen as a fraction of the largest single error contribution. We thus subdivide all edges, , with contributions with a suitably chosen constant . By multiplying with the sign of the estimate we may deal with negative errors, that is, we only refine edges with contributions increasing the absolute value of the error.

*Remark 4.1. *Considering the linear triangular elements, we map the basis functions from a reference element onto the actual element in in standard fashion. When an extremely anisotropic mesh evolves, problems can occur in the element formulation due to an ill-conditioned Jacobian. In order to circumvent this problem, we can restrain the refinement in such a fashion that we always refine the longest edge of a triangle that has become too anisotropic. We can, for instance, control the condition number of the Jacobian or the maximum ratio for the lengths of the edges.

#### 5. Numerical Examples

We will investigate two scalar problems defined in 2D on the unit square, . The first problem, the Poisson problem, has isotropic character, while the second problem, that of reaction diffusion, is constructed such that it contains (anisotropic) boundary layers. For each of the two problems, we compare the proposed anisotropic refinement strategy to those of isotropic and uniform mesh-refinement. For isotropic mesh-refinement, element indicators are used for a procedure of longest-edge splitting, as proposed by Rivara [2]. In the presented examples, a refinement ratio of is used. No shape-control of the elements is applied (or needed) in the presented results.

##### 5.1. Poisson Problem

We study the model problem The symmetric problem has a corresponding energy, which is the square of the energy norm, defined by For illustration, we choose the error measure for adaptation to be the energy of the error. This results in the well-known property that the dual solution and the error function are the same. The energy for the solution of this problem is computed as . The solution is shown in Figure 5. In Figure 6, we study the convergence of error with mesh-refinement using uniform, isotropic and anisotropic strategies. Finally, in Figure 7, adapted meshes for isotropic and anisotropic mesh-refinement are shown.

**(a)**

**(b)**

##### 5.2. Reaction-Diffusion Problem

We study the model problem
The symmetric problem has a corresponding energy
In the example, the diffusion-coefficient is chosen to be to provoke appearance of boundary layers. For illustration, we choose the error measure for adaptation to be the energy of the error. The energy for the solution of this problem is computed as . The solution is shown in Figure 8. In Figure 9, we study the convergence of error with mesh-refinement using uniform, isotropic and anisotropic strategies. Finally, adapted meshes for isotropic, and anisotropic mesh-refinement are shown in Figure 10. We note that with anisotropic mesh-refinement the mesh develops differently in the four corners of the domain. This is due to the indeterminacy of choosing how to refine an element, for which two edges of * equal* size error contributions are marked for subdivision, confer to Figure 4.

**(a)**

**(b)**

#### 6. Conclusions and Outlook

The proposed method shows promising behavior in the following respects: (i) it is more efficient than an isotropic refinement strategy for problems of anisotropic nature, (ii) its behavior is similar to that of the isotropic strategy when the nature of the problem is isotropic. Thus the method is truly adaptive, in the sense that the mesh becomes anisotropic if and only if, the nature of the error demands it, without adding any extra cost to the computation. Furthermore, the refinement procedure is extremely simple, since all element subdivisions are completely independent of their neighbors once the edges to subdivide are indicated.

For future work, the method seems appealing to extend to 3D due to (i) the efficiency in terms of degrees of freedom and (ii) the ease of implementation. If a corresponding method is to be used for quadrilateral (or hexahedral) elements, additional constraints must be applied, in order to control the legality of the iso-parametric mapping.