Analysis of marking criteria for mesh adaptation in Cosserat elasticity

The article is devoted to comparison of finite element marking criteria for adaptive mesh refinement while solving plane Cosserat elasticity problems. The goal is to compare the resulting adaptive meshes obtained with different marking strategies. Mesh refinement and error control is done using the functional type a posteriori error majorant. Implemented algorithms use the zero-order Raviart-Thomas approximation on triangular meshes. Four widely used marking criteria are utilized for mesh adaptation. The comparative analysis is presented for two planestrain problems.


Introduction
The a posteriori error control theory is one of modern directions in computational mathematics for PDE's with applications to mechanical engineering. Various error estimates are used as quantitative measure of discretization errors which take place in all numerical simulations. Corresponding error indicators help to determine subdomains with large errors for further adaptive mesh refinements.
The functional a posteriori error estimate for two-dimensional problems of Cosserat elasticity theory is considered in this paper. The Cosserat continuum [1] generalizes the classical linear elasticity theory (see, for example, [2] and [3] for mathematical statements) adding the rotational degrees of freedom independently of translational and therefore taking into account a material microstructure. Corresponding numerical methods for the Cosserat continuum intensively developed since the XXI century (see, for example, [4][5][6][7][8]). First paper concerning modern a posteriori error control for computed approximations appeared in 1994 [9]. The functional-type error estimates for the Cosserat elasticity were obtained later in [10,11] and implemented in [12][13][14]. In paper [13] first results for adaptive mesh refinement were shown.
For the classical linear elasticity the theory of a posteriori error control has evolved over the decades (see, e.g. [15][16][17][18] and [13] for reviews). The functional approach to a posteriori error control for problems of continuum mechanics developed starting from [19]. For more detailed review of functional-type error estimates for linear and nonlinear problems we refer to [18] and [20][21][22] and papers cited therein. This work continues the research started in [23] for classical linear elasticity and examines the sensitivity of adaptive mesh refinement algorithms with a posteriori error control for Cosserat continuum to element marking criterion choice.

Functional type error estimate and indicator
The investigated functional-type error estimate for the Cosserat linear elasticity in 2D has the following form [11,13] |||u -ũ||| ≤ D(ũ, s*) + R(s*), (1) where vector u contains all components of the exact solution, which is generally unknown, ũ represents the approximations of these components (e.g. the numerical solution, obtained by finite element method), s* is a set of auxiliary variables (tensor and vector fields), and |||...||| denotes the global energy norm of the error. The right-hand side of (1) is called the functional error majorant and is denoted as M. The first term D of the majorant represents errors in constitutive relations. The second term R contains the residuals and a set of meshindependent constants.
The important property of this estimate is exactness -the inequality (1) transforms to equality with the proper choice of auxiliary variables s*, the term R turns to zero and the term D turns to the estimated norm. Thus, a reasonable choice of free variables in the functional-type error estimate allows obtaining accurate guaranteed upper bound of the error norm. Also we note that the estimate is guaranteed, which means that the inequality holds for any arbitrary approximate solution ũ from the corresponding energy space, regardless of the approach used for its calculation. The majorant depends only on known data: the approximate solution, mesh-independent constants, additional variables. All auxiliary fields can be constructed with any finite elements technique suitable for space H(div) -the Hilbert space of square summable vector-functions with square summable divergence.
For mesh adaptation algorithms, the local error indicator is needed besides the global error estimate. The functional M or its first term D can be considered as a sum of values from each finite element and further used for indicating the subdomains with a high error level.

Adaptive algorithm and marking criteria
The mesh adaptation algorithm in general consists of four steps (see, for example, [24,25]) 1. Solve: compute ũ on current finite element mesh; 2. Estimate: compute the global error estimate (in our case -the functional M) and the local error indicator (contributions to M or D from all finite elements); 3. Mark: mark mesh elements with comparatively large local errors by some marking strategy (using error threshold or percentage of the total amount of elements); 4. Refine: divide the marked elements and do local mesh refinements. For marking step, a criterion for selecting the elements for splitting should be chosen. It can be represented by the operator m transforming the values of the error indicator on each element into an array of zeros and ones (zero means that element should not be split, onethe element will be split on the refinement step). We will refer to this operator as a marker. The choice of marker affects the adaptation process and produces different mesh sequences, so the resulting adaptive meshes may differ dramatically. The main goal is to obtain the mesh with as less number of nodes as possible keeping the desired computational accuracy. In this paper we consider the four most commonly used markers: 1. The maximum value marker m1: all elements are refined, for which the value of error indicator i is greater than α·max(i) for i=1,Ne; where i is an element number and Ne is the total number of elements in the mesh. Here α  (0;1) is a treshold level, typically set to 0.5. 2. The mean value marker m2: all elements are refined, for which the value of error indicator i is greater than the mean value over all set. 3. The m3 marker sets the refining of the given number or percentage of elements with the maximum error. Therefore, it is necessary to sort the elements by descending indicator value on each step. 4. The m4 marker is called the bulk criterion: the first k elements with the maximum error are refined, where k is chosen to fulfill the inequality where α  (0;1) is some treshold level. For this marker it is also necessary to sort the elements by descending indicator value. The marker comparison can be done visually: looking at obtained adaptive meshes and locations of elements chosen for refinement on each step. But to have some quantitative measure we use the marker accuracy proposed in [18] where m ref is the marker calculated for the reference error indicator on a same mesh as for m, ||...||1 is the 1-norm of the vector. The value of  allows to estimate how the number and the location of the elements selected with the given marker are close to each other when using the majorant-based indicator and the reference indicator. The reference indicator is computed with a very fine (reference) approximate solution. This is a standard approach for numerical justification of adaptive algorithms.

Numerical results
In this section we consider two examples from [13,14]. In paper [14] some adaptation results were obtained with marker m4. Both test problems are isotropic plane-strain formulations. The zero-order Raviart-Thomas approximation [26] of the auxiliary variables s* is implemented to compute the majorant M. The first term D of the majorant calculated on each element of the mesh is used as an error indicator for the approximate solution ũ.
For comparative analysis the above-mentioned reference indicator is used which is the local values of error energy norm ||| ũref -ũ||| on each element. Here ũref is the reference solution (used instead of unknown exact solution), obtained on a sufficiently refined mesh. It is commonly used to estimate the quality of error estimates and adaptive meshes. In all cases the adaptation starts with a uniform mesh with a relatively small number of nodes, and ends as soon as the given level of relative error e reaches the 2% threshold. Adaptive algorithms are implemented in MATLAB.

Example 1. Square domain with clamped boundary
For the first example we consider a test problem with Dirichlet boundary conditions and x-y symmetry. The domain is a unit square with all boundary edges clamped, the same volume force is applied in both x and y directions. The marker accuracy is computed for uniformly refined triangular mesh with 545 nodes and 1024 elements. Results are collected in Table 1.  Figure 1 shows the elements selected for refinement on a uniform mesh with four markers.
In the first column the plots constructed with the reference indicator are shown in order to estimate whether the elements have been selected correctly. In the second column the results for the majorant-based indicator are presented. We refine 20% of the elements with the maximum error for the m3 marker and for the bulk criterion m4 value  = 0.3 is chosen.
For estimating the quality of adaptive meshes, a reference mesh was constructed by refining 25 elements with the maximum error on each step. Adaptation stops when the relative error becomes less than 2%. This reference mesh can be constructed only for test problems as a numerical experiment, since it will require too many resources in real computations. The characteristics of the reference mesh are collected in the last line of Table 2.
The number of mesh nodes and elements, the relative error e and the majorant's efficiency index Ieff for each marker on the final adaptive meshes are listed in Table 2. The Nst is the number of adaptation steps. The meshes are depicted in Figure 2.  We note that the efficiency index of the majorant is not affected by the marking criterion. It can be seen from Table 2 and Figure 1 that the m2 marker leads to excessive mesh refinement due to choosing too many elements for refinement on each step. The opposite effect is observed for the m1 and m4 markers when choosing a small amount of elements for refinement leads to higher number of algorithm steps.

Example 2. Square domain with a hole
For the second example, we consider the square domain with a circular hole in the center. The left edge is fully clamped and the tensile loading is applied to the right edge. The marker accuracy for triangular mesh with 1147 nodes and 2228 elements is collected in Table 3. In Figure 3 the corresponding elements selected for refining are presented, selected with four markers. For this example, the marker accuracy for m3 and m4 is higher.
The number of mesh nodes, elements, the relative error e and the efficiency index for each marker on the final adaptive meshes and characteristics of the reference mesh are listed in Table 4. The final adaptive meshes are presented in Figure 4. In this case the m3 marker leads to most excessive mesh refinement. The highest number of algorithm steps for this example occurs when the m1 marker is used, however, the obtained mesh is close to the reference one by the number of nodes. The smallest number of nodes was obtained with marker m4.

Conclusions
The comparative analysis of the adaptive algorithms with four commonly used element marking criteria was conducted. The results show that the proper marker choice can significantly reduce the number of adaptation steps and the number of nodes in the resulting mesh. For Cosserat elasticity, the results correlate with the similar research done for classical linear elasticity [23]. We note that the main areas of node condensation for the resulting adaptive meshes do not depend on the marker choice. The bulk criterion m4 showed the most stable results giving the mesh close to the reference one within a reasonable number of adaptation steps. In addition, it should be mentioned that the number of nodes in the resulting mesh and the number of algorithm steps for the markers m1, m3 and m4 depend on their parameters. Tuning these parameters for a specific problem can improve the quality of the obtained adaptive meshes.
Research is supported by the Grant of the President of the Russian Federation MD-1071.2017.1.