A simple anisotropic mesh-refinement strategy for triangular elements in 2D

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.


Introduction
Designing the optimal finite element mesh for a certain problem is a key issue in order to make inexpensive and accurate computations.Standard h-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.

ISRN Applied Mathematics
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.

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 u ∈ V such that where u is the exact weak solution, a •; • is a semilinear form linear in the second argument , l • is a linear functional, and V 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 V h ⊂ V.In the Finite Element Method, we construct V h as piecewise polynomials on N el elements Ω e , constituting the domain Ω N elem e 1 Ω e .We thus solve the problem.Find the approximate solution u h ∈ V h such that The difference in exact and approximate solution can now be defined as e def u − u h .
We will now introduce the appropriate secant form of a •; • .Upon using the Gateaux derivative with respect to the first argument of a, we construct the secant form where u s def u h se, s ∈ 0, 1 .The form a S •, •; v, w is nonsymmetrical in general but bilinear in v and w.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 E u, u that can be any Gateaux-differentiable functional that satisfies the condition E v, v 0 for any function v. Hence, E •, • may be a norm of e a special case.Using the Gateaux-derivative of E •, • with respect to the first argument, we define the secant form of E u, u h as follows: which is a linear functional in w.We then obtain the relation We are now in the position to introduce the variational format of the dual problem.
Having computed ϕ, we can use 2.8 , 2.9 , and 2.4 to formulate the exact error representation as follows: for any ρ h ∈ V h .Here, the arbitrary function ρ h can be subtracted due to the Galerkinorthogonality 2.6 .Typically, we may choose ρ h as some projection πϕ of ϕ onto V h , for example, the nodal interpolant, the L 2 -projection, or the FE-interpolant solution .

Approximation of Error
In general, a S and E S depend on the solution u, 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 ϕ ∈ V such that where the approximations Remark 2.1.In the case of a linear problem, the secant form a S is equal to the semilinear form of the variational format and therefore 12 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 ϕ ∈ V. Following the method described in 17 , we compute the enhanced dual solution in an enriched solution space V V h ⊕ ΔV.
First, we note that the dual problem can be solved using the same FE-mesh as for the primal problem, that is, ϕ h ∈ V h is the solution of Based on the solution of the dual problem on the primary FE-mesh, ϕ h ∈ V h , the enhancement is solved globally or locally resulting in the enhancement Δϕ ∈ ΔV.Substituting the enhanced dual solution, ϕ ϕ h Δϕ, into 2.10 and choosing πϕ ϕ h , results in the following, approximate, error estimate: 2.14

Error Indicators Contributions from Discretization
Based on the error estimation above, we wish to use an adaptive h-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.

Relating Error Indicators to Edge Refinement
Usually, error indicators are obtained through restricting the error estimator 2.10 or 2.14 to each element Ω e , whereby E u, u h 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 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 {Γ f } N edge f 1 .Thus, we want to create a contribution scheme like where is the contribution to the chosen measure E •, • due to the interpolation error related to the edge Γ f , that is, the error that arises from the absence of further basis functions in the primal solution space V h .
First, we write the error contributions as a sum over the added basis functions.The enhancement space ΔV is spanned by the basis {ΔN i } N Δ i 1 of hierarchical polynomials with corresponding coefficients ΔΦ i , such that Δϕ N Δ i 1 ΔN i ΔΦ i .Thus, 2.14 can be rewritten as 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 η edge f must be reduced when edge Γ f is subdivided.

Anisotropic Distribution of Error Indicators to Edges
We now turn to the case of using piecewise linear approximation for the primal problem i.e., V h 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.
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 η Δ i to that edge.Thus the coefficients c fi for all i such that ΔN i is an edge bubble are where f * i indicates the edge on which the bubble ΔN i 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 f ei for all i such that ΔN i is an internal bubble function are where F i indicates edges surrounding the element on which the bubble ΔN i lives and |Γ j | is the length of edge number j.
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 Δϕ 0 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.

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.
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.

Refinement Strategy
The adaptive refinement is set to perform mesh-refinements until a certain tolerance is met, that is, until < TOL.When the error is larger than the preset tolerance, edges with indicators, η edge f 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, f, with contributions with a suitably chosen constant C ∈ 0, 1 .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 illconditioned 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.
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 C 0.3 is used.No shape-control of the elements is applied or needed in the presented results.

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 u 2 E ≈ 3.51 • 10 −2 .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.

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 10 −4 1 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 u 2 E ≈ 0.96.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.

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.
h ; Δϕ , 3.1 where R e •; • is the restriction of R •, • to Ω e .The contribution η elem e is related to the size of element Ω e , and in the adaptive procedure we hope to reduce η elem e by decreasing the local mesh-size.

Figure 1 :Figure 2 :
Figure 1:The conventional linear basis functions on a triangular element in 2D used for the FEapproximation.

Figure 3 :
Figure 3: Illustration of the similarity in adding a quadratic basis function a and subdividing the element with remaining linear approximation b .

Figure 4 :
Figure 4:Local subdivision of triangular elements with given edges to refine. 1, 2, or 3 edges can be marked for subdivision indicated with circles, " • " or " • " .In the case that two edges are marked for refinement the sub-triangulation depends on which of the two edges has the highest error indicator marked with filled circle, " • " .

Figure 5 :
Figure 5: Solution for the Poisson problem, illustrated by the isolines of the solution.

2 E /||u|| 2 EFigure 6 :Figure 7 : 2 E / u 2 E
Figure 6: Convergence of error with refinement Ndof no. of degrees of freedom using different meshadaption strategies for the Poisson problem.

Figure 8 : 2 E /||u|| 2 EFigure 9 :
Figure 8: Solution for the reaction-diffusion problem, illustrated by the isolines for the solution.Note the pronounced boundary layers produced by the Dirichlet boundary conditions.