Error analysis for quadtree-type mesh coarsening algorithms adapted to pixelized heterogeneous microstructures

For pixel-based microstructure representations we propose two different variants of adaptive, quadtree-based mesh coarsening algorithms that serve the purpose of a preprocessor for finite element analyses in the context of numerical homogenization. High resolution is preserved at interfaces for accuracy, coarse-graining in the defect-free interior of phases for efficiency. Error analysis is carried out on the micro scale by error estimation which itself is assessed by true error computation. Modified stress recovery schemes for an error indicator are proposed which overcome the deficits of the standard superconvergent recovery scheme for nodal stress computation in cases of interfaces with stiffness jump. By virtue of error analysis the improved efficiency by the reduction of unknowns is put into relation to the increase of the discretization error and thereby sets a rational basis for decisions on favorable meshes having the best trade-off between accuracy and efficiency as underpinned by various examples.


Introduction
Pixel-and voxel-based representation of reconstructed microstructures obtained from tomographic imaging methods (scanning transmission electron microscopy tomography (STEM), focused ion beam tomography (FIB) and X-ray computed tomography (XCT)) finds nowadays an ever more increasing interest not only in research but also in industrial applications [23,25,39]. Reconstructed microstructures can directly be analyzed for characteristic quantities such as volume fractions of phases and pores, size distributions of pores, phases and grains, interface and surface areas, constrictivity, tortuosity, and many more. Furthermore, the pixel/voxel meshes can be used as input for numerical analysis in general, for finite element simulations in particular, in either case for the prediction of effective physical properties in silico. Corresponding voxel to mesh conversions were driven by the representation of bone [22,26], for a general overview see [33].
B Bernhard Eidel bernhard.eidel@uni-siegen.de 1 Heisenberg-Group, Institute of Mechanics, Department Mechanical Engineering, University Siegen, Paul-Bonatz-Str. [9][10][11]57068 Siegen, Germany The corresponding highly resolved, uniform discretizations are well-suited to accurately describe the geometry of interfaces and defects in microstructures and, consequently, to capture the physical processes in these regions of interest. For the defect-free interior of phases and grains with their smooth field properties at slow variations, the high resolution is typically a wasteful luxury and coarse-graining adequate. The present work proposes two different variants of adaptive algorithms for quadtree-based mesh coarsening of the original, uniform, pixel-based meshes that serve the purpose of a preprocessor for consecutive finite element analyses. The criterion for mesh coarsening is enabled for pixels sharing the same phase or grain. We showcase the applicability in various examples and measure the improved efficiency in terms of the reduction of unknowns and a corresponding speed-up in computational homogenization.
Quadtree-/octree-based mesh coarsening of heterogeneous, multiphase materials is frequently used, see e.g. [24,28,29,34,41] for application in homogenization and fracture, in the context of the Scaled Boundary FEM [19,45,46] and many more. In some cases the computational results for the coarse-grained and the uniform discretizations have been checked by visual inspection or by a comparison of global errors, in case of analytical solutions by error computation.
But what is still missing is a thorough error analysis including spatial distributions of errors.
The present work fills this gap by error estimation which itself is assessed by true error computation, the latter based on overkill solutions. In the context of error estimation we propose two different, modified stress recovery schemes which overcome the deficits of the standard, superconvergent scheme of Zienkiewicz and Zhu for the case of multimaterialinterfaces with stiffness jumps and corresponding stress jumps. We show that the standard scheme spoils error estimation in that case and that the modified recovery leads to a very good agreement of a posteriori estimation with exact error computation.
By virtue of error analysis we can weigh the pres and cons of adaptive mesh coarsening; the improved efficiency by the reduction of unknowns is set into relation to the increased discretization error. Hence, the quantitative error analysis sets a rational basis for deciding on best trade-offs between accuracy and efficiency for the adaptively coarse-grained meshes.
In the present work the pixelized microstructure images are understood as representative area elements (RAE) each of a microheterogeneous multiphase material. They are used in a computational framework for the micro-to-macro scale transition based on mathematical homogenization [5,8,51]. This scale transition is carried out by a two-scale finite element method referred to as Finite Element Heterogeneous Multiscale Method (FE-HMM) introduced as a particular case of the general HMM in [54]. For more comprehensive descriptions of the method see [2,3,53]. A formulation of FE-HMM for linear elastic solid mechanics based on [1] with a conceptual and numerical comparison to the FE 2 method is described in [13], a short precursor in [12]. The numerical analysis for coupling conditions fulfilling the Hill-Mandel postulate of micro-macro energy densities [20,21] is presented in [17]. For an extension to geometrical nonlinearity and hyperelastic materials see [14]. FE-HMM's methodological twin is the FE 2 method, [16,27,31,32,40], with overviews e.g. in [18,42,50].
By the embedding of the pixelized microstructures as RAEs into a macro continuum through homogenization, another relevant aspect of error analysis is how the micro discretization error propagated to the macro scale influences macroscopic properties, here in terms of the effective elasticity tensor and macroscopic deformations.
For the mentioned identity of FE-HMM and FE 2 the achieved results of the present work are equally valid for both methods.

Mesh coarsening
The present section introduces two different algorithms for quadtree-based mesh coarsening. Mesh coarsening is applied to pixelized microstructures and works as a preprocessor generating input data for consecutive finite element analyses. Since each pixel exhibits a unique color code, different phases can be distinguished without ambiguity.

Quadtree-based mesh coarsening
The general procedure of mesh coarsening shall be outlined. As the name "quadtree" already indicates, the origin of all computations is a micro mesh with a uniform triangulation T 0 , meaning that each node in the interior of the micro mesh has exactly four neighboring nodes and, consequently, each node is shared by four elements. All elements of the micro mesh initially exhibit the same element side length denoted by h. Figure 1a shows an example for a quadtree micro mesh with elements of side length h. There are four nodes in the interior of the mesh, each with four neighboring nodes and each adjacent to four elements. In quadtree-based mesh coarsening the set of four elements adjacent to one micro node can be merged to one new element of side length 2h. One option for suchlike coarsening is shown in Fig. 1b. Here, the elements around the top left of the inner nodes are coarsened to one new element. Similarly, the elements around each of the other internal nodes could have been coarsened.
Mesh coarsening reduces the number of elements and thereby the number of degrees of freedom (ndof). The original uniform mesh contains 32 degrees of freedom (16 nodes, each with 2 degrees of freedom). Mesh coarsening removes three nodes from the uniform mesh (a) which results in the new mesh (b). Furthermore there are two nodes, which do not have four neighboring nodes any more, because they lie on edges of the coarsened element. To preserve continuity at the element edge, these nodes have to be coupled to the coarsened element by a constraint. These nodes are referred to as "hanging nodes", their constraints read as (1) Fig. 1 Quadtree: a original, uniform quadtree mesh, b quadtree mesh with one coarsened element with the two master nodes at the end of the edge to which the hanging node is connected. Due to the constraint the hanging nodes do not have any degrees of freedom left. The constraints are realized in analogy to the micro-macro coupling by introducing further Lagrange multipliers. For details about the micro-macro coupling see [17]. This rationale to merge four elements having a node in common into one coarser element shall be used in an efficient preprocessor for finite element analyses. Algorithms are needed to decide, whether an element is target of coarsening or not.
Finite element mesh generation based on quadtrees was introduced by Yerry and Shephard [55], for comprehensive overviews of original works we refer to [43,44].

Coarsening algorithms
For pixelized microstructures with two or more phases, the elements in the interior of the phases should be marked for coarsening, while the elements at the phase boundaries should preserve their initial, fine discretization. The coarsening criterion is defined as follows. The algorithm contains an if-statement, which distinguishes whether an element is marked for coarsening or not. The if-statement contains two conditions, which have to be fulfilled to allow the element to be marked for coarsening.

Algorithm 1 Mesh coarsening criterion
The first condition screens for each micro element, if one of its nodes belongs to a phase boundary. In the case that the element is directly adjacent to a phase boundary, it is not marked for coarsening. Otherwise the element is located inside a phase and target of coarsening.
Typically, a coarsening step is executed more than once to efficiently reduce the number of elements and thereby the number of degrees of freedom. This may lead to 2 n −1 hanging nodes on an element edge after n coarsening steps. To prevent multiple hanging nodes on an element edge in additional coarsening steps, the second condition in Algorithm 1 is introduced which scans, if one of the element nodes is already involved in a hanging node constraint. If the element should not be further coarsened, the number of hanging nodes on one edge is limited to one. In the context of quadtree mesh-ing the case of at most one node in a side is referred to as a balanced leaf.
The mesh coarsening following Algorithm 1 might lead to steep gradients in element size at phase boundaries with a densely packed cluster of kinematical constraints, which artificially stiffens these regions. Less restrictions at boundaries imply less steep gradients in element size. For that goal the conditions from Algorithm 1 are slightly modified. Instead of scanning if one of the element's nodes is either a node at the phase boundary or is used for the definition of a hanging node constraint, the modified Algorithm 2 scans, if one of the element's nodes is either adjacent to an element at a phase boundary or is adjacent to an element which is used for the definition of a hanging node constraint. Figure 2a shows a uniform mesh with elements belonging to two different phases, a yellow and a green phase. All elements marked with a black cross in Fig. 2b are marked for coarsening by both Algorithms 1 and 2. The elements marked with red crosses are only marked for coarsening following Algorithm 1. The once coarsened meshes with crosses marking the elements for the second coarsening step are displayed in Fig. 2c following Algorithm 1 and Fig. 2d following Algorithm 2. In both meshes only few elements are marked for further coarsening, following Algorithm 1 at least the four elements in the upper left corner may be coarsened to one bigger element, displayed in Fig. 2e.
The mesh coarsening following Algorithm 1 leads to a steep gradient of element size at phase boundaries, which results in strong restrictions due to the hanging node constraints. If the red marked element in Fig. 2e is considered, each of its nodes is obviously included in a hanging node constraint.
Algorithm 2 leads to a much smaller gradient in element size than Algorithm 1. It is obvious that less elements vanish by coarsening, such that the number of degrees of freedom will be somewhat larger for the mesh in (d) compared to the mesh in (c), but the reduced number of kinematical constraints leads to better, i.e. less stiff results.
The two coarsening algorithms will be investigated for several examples in Sect. 4, they will be referred to as hard coarsening (following Algorithm 1) and soft coarsening (following Algorithm 2).

Fig. 2
Quadtree mesh coarsening at material interfaces: a original mesh, b original mesh with marked elements following Algorithms 1 (black and red crosses) and 2 (black crosses only), c once coarsened mesh following Algorithm 1 with elements marked for next coarsen-ing step following Algorithm 1, d once coarsened mesh with marked elements for further coarsening following algorithm 2 and e resulting twice coarsened mesh following Algorithm 1. (Color figure online)

Error computation and error estimation
The aim of the present work is to assess adaptive mesh coarsening as a preprocessor applied to pixelized microstructure images. To compare the results of the coarsened micro meshes to the results of the original, uniform micro mesh, the corresponding micro discretization errors will be investigated. Here, error analysis is based on error computation obtained through a reference solution for a superfine discretization (h → 0) and based on a posteriori error estimation, where the latter will be validated by the former.
For error computation and error estimation we will use the energy norm. The energy norm reads for the solution on the micro domain where A is the linear elasticity tensor on the micro scale showing heterogeneity on the micro domain for the multiphase characteristics of microstructure in our investigations.

Error computation
The integrals for error calculation in the energy norm (2) are approximated by numerical integration of Gauss-Legendre. The computations are carried out on micro element level of the discretization for the reference solution. For the error in the energy norm it follows In this case the results of the FE-HMM solution on the coarse meshes are projected onto the finer grid of the reference solution. Figure 3 provides a sketch of the projection from a coarse discretization onto a finer reference discretization for one quadrature point of an element of the reference mesh. Therein, the quantities of the coarse mesh are projected onto the quadrature points of the reference solution such that the error of the quantities of interest can be calculated, e.g. for the error of stress and strain in the energy norm according to (4).

Stress computation for error estimation
In (engineering) practice, error computation is prohibitive. Instead, an a posteriori error estimation is carried out for the particular discretization in use. The class of reconstructionbased error estimation relies on the computation of improved stresses (and strains for the error in the energy norm) at the nodes of a finite element; the integral deviation of suchlike reconstructed stress distribution from the distribution of inner-element stresses at its quadrature points in a suitable norm (typically energy norm or L 2 -norm of stress) yields the elemental error estimate called error indicator [56][57][58]. If bounds cannot be proved, it is referred to as refinement indicator instead of error indicator. The deviation can be understood as a residual (vanishing for a characteristic element size h → 0) that sets the link from reconstruction-based to residual-based error estimation [6]. This relation has been established for particular settings, see for example [56] with the proof of equivalence by Rank in the appendix of that reference; for a comprehensive overview of error estimation in finite element analysis we refer to [4].
The particular aim here is to apply the reconstructionbased error indicator of Zienkiewicz-Zhu in an appropriate format to multiphase microstructures, the material class of interest in the present work.

The standard superconvergent patch recovery (SPR)
The point of departure is the error estimation of Zienkiewicz and Zhu which relies on superconvergence of stress and strain. Superconvergence in the present context refers to the property of strain and stress to converge in the same order as the displacements do O(h q+1 ) with the polynomial order q of the shape functions. This is an exception to the nominal order q for strain which follows from differentiation for strain calculation. A necessary condition for superconvergence in quadrilaterals is their rectangular shape. In that case along with sufficient regularity, superconvergence can be measured for q=1 in the element center, for q=2 in the 2×2 Gauss points; these points are referred to as Barlow-points [7]. Notice that pixelized images used here as initial, uniform discretization as well as their quadtree coarsened offsprings each exhibit a quadratic shape and therefor fulfill the necessary condition for superconvergence. In [57] a procedure was introduced that transfers the superconvergence property from the Barlow points in the element interior to element nodes; it is referred to as "superconvergent patch recovery" (SPR). These recovered superconvergent nodal values enter the error indicator and thereby largely determine its accuracy; the error distribution is the prerequisite for adaptive mesh refinement [58].
By applying the SPR to calculate nodal stresses, a continuous stress distribution is generated. For finite elements with C 0 -continuity the displacement field is continuous, stress and strain is not. Instead, they exhibit jumps at nodes, which is even more pronounced, if a particular node is right at the interface of different materials having a large stiffness mismatch. The original SPR does not account in any form for phase information, so that there are some modifications needed to make it applicable to these structures.
Since the generation of classical patches in the format of SPR is not possible for multiphase structures in any case, a different scheme is introduced for reconstructing nodal stress and strain from their counterparts in the interior of the element.
For ready reference, the rationale of the SPR is briefly reiterated for linear shape functions q=1 in a 2D setting, where strain and stress are calculated each in the center of rectangular elements. These values are transferred by a least-square procedure to the finite element node in the direct neighborhood. Elements having such a node in common are referred to as patch in the recovery procedure.
Stresses on the patch are prescribed component-wise by with, for the case of linear shape functions, Vector P contains polynomial terms of bi-linear shape functions for n dim = 2. For the determination of the unknown vector a the function has to be minimized. Therein, (x i , y i ) are the coordinates of the superconvergent points, n is the number of superconvergent points of the total patch and σ h * (x i , y i ) are the stresses in these superconvergent points. Minimization of F(a) implies that a fulfills the condition which can be solved for a Recovery of nodal stresses using standard SPR and modified SPR: a standard SPR for a patch with only one material phase and b standard SPR applied to a two phase material patch ignoring the interface, c modified SPR applied to a two-phase material with explicit account of the element phase by introducing two sets of stress at interface nodes Stresses in the central node of the patch can be recovered by inserting its nodal coordinates (x N , y N ) into the P-vector in (5).

The modified superconvergent patch recovery
The standard SPR does not account in any form for different phases and interfaces in patch generation. For all nodes right at phase boundaries the corresponding standard patch is made up of elements belonging to different material phases. Figure 4a gives a sketch of the standard SPR for a patch with elements from one and the same material phase, Fig. 4b gives the same sketch for a node at the boundary between a stiff (green) and a soft (yellow) phase. The stiffness jump at the phase boundary induces a stress jump, an effect which is the more pronounced, the larger the stiffness mismatch. In the example of Fig. 4b the stiff (green) phase exhibits considerable larger stress values than the softer (yellow) one. Fitting a bi-linear interpolation function through stress and strain at the support of the patch results in a reconstructed stress at the central node which is an average of the values of different magnitudes and thus ignores -i.e. falsifies-the true mechanical phenomenon of a discrete stress jump. In the context of error estimation based on intentionally improved nodal stresses in comparison to stress in the element interior, suchlike falsified nodal stresses yield a spoiled support for interpolation, a fact that was already recognized in [35].
A solution to the problem is the assignment of two sets of stress for nodes at interfaces; it accounts for the mechanics and, as it will turn out, improves error estimation.
In order to apply this modified SPR to the central node at the phase boundary from Fig. 4b, two patches have to be  Figure 4c displays the patches and the stresses following from each patch. Fitting the bilinear interpolation function (blue) through the stresses of the supports of the stiffer phase (green) obviously leads to a higher stress value for the central node at the phase boundary than fitting the bi-linear interpolation function (red) through the supports of the weaker phase (yellow). By introducing the duplex stress values, both sets of stress are available for further use.
Expanding the patch into the phase interior is in perfect analogy to the treatment of edge and corner nodes at the mesh boundary for a single-phase material [17].
Nevertheless, there are several cases of element arrangements at interfaces which prevent the generation of standard patches. One case is the scattered element phase distribution which is characterized by a non-straight interface line as illustrated in Fig. 5a. The node marked in blue has only one single adjacent element from the yellow phase. If the patch is expanded to the interior of the phase, there are still only three elements for constructing the node's patch. The yellow element from the upper right of the displayed section of the micro mesh is not a candidate for patch-construction, since there would be three points on a straight line, making matrix A in (9) singular.
There are further situations conceivable, where even less than three elements are available for patch generation. Imagine the case of one (or two connected) elements, which are placed in between elements from another phase. In this case the patch has to be build up by one (two) elements. A similar situation arises for hanging nodes, since there are always exactly three elements surrounding them, see Fig. 5b.
To treat patches with less than four elements the number of terms in vector P in Eq. (6) 1 may have to be reduced, since otherwise the number of superconvergent points is not sufficient to determine the parameters of vector a in (6) 2 .
For interface nodes at the corner of two patches belonging to the same phase, nodal stress is computed by averaging the values obtained from the involved patches according to (14).

Averaged stresses and strains based on "intra-element patches"
The main disadvantage of the SPR is the fact, that -especially for adaptively refined meshes and meshes with multiple phases-regular patches based on four surrounding elements can not be properly constructed for all nodes.
To avoid too heavy computational efforts in generating sufficiently regular patches, a different but similar method is proposed and will be investigated. For this method stress and strain from the standard quadrature points are used instead of stress and strain from superconvergent points. Using these stresses and strains and the element shape functions, nodal element stresses (σ node,ele ) and strains (ε node,ele ) can be calculated.
For each quadrature point with coordinates x qp,i with i = 1, . . . , n qp it must hold Written up in vector-matrix-representation Eq. (11) reads as Nodal element stresses and strains for each element are obtained by The general procedure of using the shape functions to extrapolate stresses and strains from quadrature points to the element nodes-as done in (11)-(13)-is standard. Nonstandard is that in the averaging for nodal values according to (14) only those contributions of the same material phase are considered instead of the contribution of all elements with n ele,phase the number of adjacent elements from the same material phase. This procedure is quite similar to the patch recovery method, since it is basically the same as fitting a bi-linear polynomial through the four standard quadrature points thus giving rise to the term "intra-element patch". By doing so each node obtains from the adjacent elements a set of elementwise stress and strain. Figure 6a shows a patch with two different phases, where for each of its elements the bi-linear shape functions are fitted through the quadrature points (a). The central node exhibits four stress values, one from each of the elements. The duplex stress related to the stiffer (green) phase is averaged from the three nodal stresses following from the bi-linear shape functions (blue) fitted through the quadrature points of the stiffer elements. Analogously, the duplex stress related to the weaker (yellow) phase is directly derived from the bi-linear shape function (red) fitted through the quadrature points of the one element belonging to the weaker phase. Figure 6b displays the resulting stress, the central node now exhibits exactly two sets of stresses, the duplex stresses following from the two different phases.
By doing so the costs of patch construction vanish, whereas nodal stress and strain are still calculated via bilinear polynomials.

Error estimator and effectivity index
As described above, recovery-based error estimation is based on improved nodal stresses and strains, denoted by σ and ε . They are computed by the SPR or by the averaged nodal stresses and strains based on intra-element patches.
The estimated error in the energy norm on the micro domain reads as Averaging method for a node attached to four elements with different phases: a the central node obtains from all adjacent elements a set of elementwise stress and strain, b the duplex nodal values result from averaging stresses and strains each of the same phase according to (14) Here the superconvergent nodal stresses and strains have to be evaluated in the quadrature points x h i . The procedure of calculating quantities in the quadrature points from nodal values is executed using the shape functions, analogous to (11).
Compared to the error computation based on a nominally exact reference solution according to (4), the error estimation is clearly much cheaper, since the integration of the error is carried out on the same mesh with triangulation T h instead of on another, reference triangulation T ref with h → 0, which needs an extra, expensive computation.
The quality of the error estimator is typically assessed by the so-called effectivity index θ which is defined as the ratio of the estimated error ē to the true error e For consistency the effectivity index must tend to unity as the exact error tends to zero.

Numerical examples
In this section a thorough error analysis is carried out for both the original, uniform discretization of pixelized images and the discretizations following from consecutive mesh coarsening.
Since microstructures are considered representative volume elements (here in the 2D setting: representative area elements RAE) in a two-scale finite element method for computational homogenization such as FE-HMM and FE 2 , all quantities in the present numerical analysis refer to the micro scale such that the attribute "micro" can be silently omitted.
In the analysis we will focus on  To investigate the discretization error on one specific RAE the corresponding macroscopic displacement field has to be kept constant for all micro discretizations, otherwise there would be an influence from the macro quantities. Here the macroscopic displacement field obtained from the original, uniform mesh in the RAE will be used and only the postprocessing is executed for all micro meshes.
Remark For linear elasticity it is possible to precompute the effective elasticity tensor for it is constant. Then the macroscopic problem can be solved by a standard, single-scale finite element method. The error analysis on the micro scale however, which is the main topic of the present work, requires a separate consideration of the fully resolved microstructure on the RAE level. Moreover, linear elasticity shall be understood as a simple prototype model for the more general case of hyperelasticity, where the tangent moduli are deformation-dependent and therefore do not allow for calculations beforehand.

Cross
In the first example a microstructure of a cross-shaped, stiff inclusion in a softer square matrix is considered. For a periodic metamaterial of tetragonal symmetry with similar unit cells, elastic material parameters of a micromorphic mate- rial model were identified by nonstandard homogenization [36], which captured the observed band-gap properties fairly well [9]. The error calculations will be executed on different reference meshes to investigate (i) the convergence of the calculated micro discretization error and (ii) the discretization which is sufficiently fine to serve for a reference solution. Furthermore, (iii) the estimated errors employing different reconstruction schemes are compared with and assessed by the calculated errors which yields their effectivity index as defined in (17). Figure 8a shows the cross microstructure with its two different phases, silicon carbide SiC (red) and diamond (white), Fig. 8b shows the micro mesh after three coarsening steps employing the soft coarsening method. Table 1 provides the number of degrees of freedom of the original, uniform mesh (128 × 128 elements) and the meshes resulting from the coarsening steps, where coarsening step 0 refers to the original mesh. After three coarsening steps the number of degrees of freedom is reduced to less than 12% compared with the initial mesh. Furthermore the number of constrained degrees of freedom is given. The initial mesh is uniform and thus does not exhibit any deactivated degrees of freedom, after mesh coarsening the ratio of constrained ndof to all ndof increases from 5.0% for the first coarsened mesh to 18.5% for the third coarsened mesh.

Mesh coarsening
Here and in the following the number of degrees of freedom in the tables equals the total number of degrees of freedom of the micro mesh minus the number of deactivated degrees of freedom of the hanging nodes. Since different coupling conditions on the RAE imply additional but different reductions of unknowns, they are not yet subtracted from the tabulated figures to allow for a streamlined comparison. Table 2 shows the errors in the energy norm for the original, uniform mesh and the meshes from the three coarsening steps based on reference solutions on different discretizations. To achieve reliable results, the discretization of the reference solution should at least be eight times finer per dimension of space than the discretization of the original, uniform mesh.

Computed versus estimated micro discretization error
For future use we validate error estimation by the true error computation and thereby distinguish between different schemes of error estimation/stress computation, (i) the standard SPR blind for phase distributions, (ii) the modified SPR accounting for interfaces, and (iii) the averaging of elementwise quantities belonging to the same material phase. Table 3 displays the results for the above three different methods (i)-(iii). As expected, the standard SPR, blind for interfaces inside the patch, leads to gross errors, an overestimation of the true error by more than 60% for the uniform mesh.
The estimated errors based on the modified SPR and on the averaging of elementwise quantities exhibit an excellent agreement with the true errors as indicated by effectivity indices θ close to one. The results suggest that error estimation is applicable to multiphase microstructures having interfaces along with a stiffness jump, if the two phases are strictly distinguished in the calculation of nodal stresses and strains.

Error distribution
For the investigation of error distributions on the micro domain the relative elementwise micro discretization error is analyzed. This relative error is the ratio of the computed or estimated error of an element to the element's energy norm. By doing so the influence of the element size is eliminated which is especially important for the non-uniform meshes. Figure 9 shows the distribution of the estimated error computed with a micro mesh of 1024 × 1024 elements with different color scaling. The original scaling in the left indicates that the maximum relative discretization error of about 12% is located in and confined to the cross' corners. The cross' edges are hardly visible for the errors virtually vanishing-similar to the interior of the two phases. In order to closer investigate the micro discretization errors, the range of the color scaling is changed to [0−1%]. The error distribution with reduced colorbar range underpins that there are only small spots with relative errors of 1% or above at the cross' corners while even the edges exhibit an error of less than 1%.

Remark
Since there are only a few, small spots with errors of above 1%, the overall discretization error is small. For efficient mesh coarsening it is preferable to have a contrast in error distribution, major ones all over phase boundaries and minor ones in the phase interior. In that case, element coarsening in the phase interior leads to an increase of these minor errors, but since the dominant errors at phase boundaries remain unaltered, the overall error does not strongly increase. In this particular example however, the error at the phase boundary is not dominant enough to allow for efficient mesh coarsening. The estimated overall micro discretization error on the micro mesh with 1024 × 1024 elements is 0.0022 [F L], the estimated micro discretization error of the interface elements is 0.0014 [F L]. This is not sufficient for effective mesh coarsening and leads to an already strongly increased error between the original, uniform mesh and the first coarsened mesh, see Tables 2 and 3.

Conclusions
For the excellent accuracy of error estimation in the present example we will discard the expensive, and in engineering practice hardly feasible error computations in the following. From the computational point of view the averaging technique for nodal stress/strain computation is the simplest one of the three techniques, since no patch generation around each node is required. For its poor results ignoring interfaces with their stress jumps, the standard SPR disqualifies for multiphase error estimation. For that reason we continue our analysis in the following with a focus on the modified SPR and the error estimation based on intra-element patches along with averaging.

DFG-Heisenberg
The second example is a two-phase microstructure showing a portrait of the physicist Werner Heisenberg 2 as shown in Fig. 10a. The white phase exhibits the elasticity properties of SiC, the blue phase those of diamond. The original, uniform micro mesh exhibits 896 × 896 elements. In this example the different coarsening algorithms from Sect. 2 are compared, on the one hand the hard mesh coarsening which leads to large gradients in element size at phase boundaries and on the other hand the soft mesh coarsening, where the element size gradient is reduced. The resultant mesh after five soft coarsening steps is displayed in Fig. 10b.

Mesh coarsening and discretization error
For its very fine initial discretization compared to the previous example, the DFG-Heisenberg microstructure is subject to more coarsening steps, but restricted to five to avoid too large elements at the RAE boundaries. Figure 11 shows the error increase for mesh coarsening over the corresponding ndof for the soft and the hard coarsening techniques; the estimated error is obtained by the averaging method and the modified SPR. Both mesh coarsening algorithms yield an almost constant increase of the estimated errors from step to step; since, however, only the first two steps considerably reduce the ndof in contrast to the consecutive ones, the gain in efficiency versus the accuracy loss pays off only in these first two steps.
The diagram indicates that the second mesh obtained by hard coarsening and the fourth mesh obtained by soft coarsening coincide almost quantitatively for the same number of unknowns in their total error. For a comparison of the two coarsening algorithms a focus is on these two micro meshes in the following. Table 4 shows the number of degrees of freedom of the original, uniform micro mesh and its coarsened counterparts. Especially during the first two coarsening steps the number of degrees of freedom is strongly reduced while the last two coarsening steps reduce them less than 1%. The soft coarsening algorithm reduces the number of degrees of freedom to less than 12%, the hard mesh coarsening even to less than 7%. The ratio of deactivated degrees of freedom is 4.5% (soft) and 5.2% (hard) for the first coarsened mesh and increases to 16.3% (soft) and 28.6% (hard) for the fifth coarsened mesh.
Furthermore, Table 4 shows the estimated errors in the energy norm for both error estimation techniques and both mesh coarsening algorithms. Similar to the first example the estimated errors by the modified SPR are larger than those of the averaging method.
Soft mesh coarsening after its fifth step reduces the number of unknowns by about 88% for an error increase of 27%, when the averaging technique is used, and by 19% for patch recovery. Hence, the considerable efficiency gain need not be paid by a comparable accuracy loss. The hard mesh coarsening leads after five coarsening steps even to a reduction of unknowns of more than 93%, which results however in a major error increase; namely of 56% for the modified SPR and 70% for averaging.

Distributions of error and micro strain
The distribution of relative element errors in the energy norm is the first quantity of interest in the comparison of Fig. 12a-c. For (a) the major errors are located at the phase boundaries, while they vanish in the phase interior. Only at those interfaces of the letters "DFG", which are aligned with the xand y-axes, the errors vanish likewise, which is equally true for cases (b) and (c). For (b) there is an increased error in the phases' interior, while the major error is still at phase boundaries. The maximum elementwise error is in quantitative agreement with the one of the initial, uniform discretization. The error distribution for case (c) exhibits minor errors in the phases' interiors compared to case (b) of the fourth soft coarsened micro mesh and resembles therein the original, uniform mesh. In exchange, the errors at the phase boundaries are not as sharp as for both other investigated meshes, resulting in a maximum elementwise error which is slightly decreased.
Complementary to the error plots the high accuracy in the interior of the phases can be visually underpinned by the contour plots of strain, here in terms of component ε x x as displayed in Fig. 12d-f. The results reveal that strain ε x x obtained for the soft coarsened mesh is in good agreement with the uniform mesh in its general distribution and in particular in the maximum, while the hard coarsened mesh shows some minor deviations which are attributed to the increased number of hanging node constraints close to the phase boundaries.

Conclusions
For the DFG-Heisenberg microstructure a strong reduction of unknowns is obtained in the first coarsening steps, while the fourth and fifth coarsening steps do not pay off any more in view of the modest efficiency gain and the rather strong accuracy loss. For the soft coarsened mesh the distributions of discretization error and strain are in quantitative agreement with the uniform mesh almost everywhere, while the hard coarsened micro mesh exhibits some minor deviations.

Circle
The next example is a stiff circular inclusion in a soft matrix. The microstructure is chosen for its high regularity since it does not contain any sharp corners, where one phase might sting into the other one. While we have used periodic coupling so far, mesh coarsening in the present example shall be analyzed for different coupling conditions on the RAE boundary, for the cases of Dirichlet, Neumann and periodic coupling. A detailed error analysis for different coupling  conditions along with (almost) uniform discretizations is presented in [17], for the 3d counterpart with adaptive resolution see [30].
In this example we combine the averaging method for error estimation with the soft mesh coarsening algorithm. For the good agreement of both estimation methods (i.e. averaging with modified SPR) at least in the first coarsening steps of the previous examples, we can focus on one of them. Furthermore, the soft mesh coarsening algorithm is chosen for its smaller errors in its first coarsening steps compared to hard coarsening. Figure 13 shows the circular inclusion microstructure with its two different phases and the fifth soft coarsened mesh; the red phase is chosen to be SiC, the yellow phase is diamond.

Mesh coarsening and micro discretization error
The original, uniform micro mesh of the circular inclusion problem is discretized by 1280×1280 elements. The number of coarsening steps is again restricted to 5 steps. Table 5 shows the number of degrees of freedom for the original, uniform microstructure and after the single coars-ening steps with their reduction factors. The fact that the inclusion is rather small in both its size and its boundary length enables a very efficient mesh coarsening; in the fifth coarsening step the total number of degrees of freedom is reduced to less than 1.5% of the original number of degrees of freedom. The ratio of deactivated degrees of freedom is 0.6% for the first coarsened mesh and increases to 17.4% for the fifth coarsened mesh.
The estimated micro discretization errors are displayed in Table 5 for different micro-macro coupling conditions. Especially the calculations with the last two coarsened micro meshes show major deviations to the original, uniform mesh, but the number of micro elements is decreased by more than 98% compared to the original, uniform mesh.
The diagram in Fig. 14a reveals the impact of different coupling conditions on the total error for mesh coarsening. The corresponding estimated errors all show the same characteristic trend for mesh coarsening. Dirichlet coupling leads to the smallest error and has the least number of remaining degrees of freedom in the system of equations while Neumann coupling leads to the biggest estimated error and has the highest number of remaining degrees of freedom.

Convergence analysis
For microstructures in pixelized discretizations the questions arise, whether existing corrugated interfaces are an artefact of a limited pixel resolution or whether they are real; moreover, whether the re-entrant corners are of harm for the regularity of the problem as indicated by a reduced convergence order.
Standard finite element discretizations with their alignment at interfaces tend to make smooth what is possibly non-smooth, pixel-based discretizations tend to introduce corrugated interface representation including singularities due to re-entrant corners.
The first question must be answered separately for each case, the second question is addressed in the following for the matrix-inclusion problem; its interface is defined to be smooth. A pixel-based discretization and a standard finite element discretization are shown in the insets of Fig. 14b, c with a zoom into the interface. A convergence analysis is executed by uniformly refining the initial meshes. The convergence results in Fig. 14b, c are in almost quantitative agreement for the two cases; they fall somewhat behind the nominal order of q=1. In the present setting of a smooth interface the pixelized interface representation with its re-entrant corners does not spoil convergence.
This result can not be generalized. Imagine the limiting case of a coarse pixel resolution where the circular inclusion degenerates into a square. In that case the regularity of the BVP is reduced where the loss of order depends on the relative inclusion size and stiffness mismatch between the phases as shown in [13,17].
Remark Why do total errors depend on boundary conditions? All diagrams in Fig. 14 show maximal errors for Neumann coupling, smallest for Dirichlet, PBC is in between. Why is that and why in this order? Notice that for heterogeneous matter the RAE boundaries typically show for constant traction conditions larger small-scale displacement variations than for the "long-wave" PBC kinematics and much more than the straight, linear displacement condition. So the fulfillment of constant tractions is most demanding for a fine spatial resolution among these three. Vice versa, for a given discretization the errors at the boundaries are expected to be largest for Neumann, smallest for Dirichlet, for PBC in between. This rationale is underpinned by the simulation results in Fig. 15.
Although each example in the present work exhibits this order in error size, it can not be generalized as counterexamples indicate [17].

Distributions of error and micro strain
Again, the distribution of the micro discretization error on the micro domain will be investigated. Since the first examples including the present one have cleared that the last coarsening steps do not pay off, we will stick to the third coarsened mesh and compare the results for different coupling conditions to the results of the original, uniform mesh.
The error distributions as well as the normal strain component ε x x in Table 6 show a very good agreement of the coarsened discretization with the uniform one in both the maximum elementwise figures and the distribution.

Conclusions
The high regularity of the circular inclusion and especially the small length of the boundary between the two phases leads to a very efficient mesh coarsening; after three coarsening steps the number of micro elements is reduced by more than 97% along with small discretization errors and good agreement in strain. Furthermore, different coupling conditions on the micro domain boundary coincide qualitatively in a considerable reduction of unknowns in the first coarsening steps while the discretization error does moderately increase. Again, the last coarsening steps do not pay off for all coupling conditions. The softer the coupling, the larger the error, which is true for the initial, uniform mesh and for the error evolution in coarse-graining. Figure 16 shows microscopy pictures already converted into binaries of a Diamond/SiC thin film microstructure, which was fabricated by chemical vapor deposition (CVD). For the computation of effective elasticity properties through homogenization we refer to [15].

Diamond/SiC microstructure
For the computations an area element (marked with red frame in Fig. 16) is chosen which is statistically representative for the basic population of the full image in Fig. 16a with respect to the phase fraction ratio (white phase 22.6% for the full image, 21.1% for the RAE) and the effective elasticity properties (deviations in the coefficients throughout less than 4.5%).
Again we employ averaging in error estimation along with the soft mesh coarsening algorithm. Notice in Fig. 16 that the discretization of the microstructure is relatively fine in the upper and right edge of the RAE compared to its lower and left edges. This is due to the choice-here on purpose-of a sub-optimal initial discretization having 1098 elements per edge where 1098 exhibits 2 in its prime factorization only once. All error data in 10 −2 [F L] Table 8 Diamond/SiC microstructure: distribution of (lines 1 and 2) the relative micro error ||ē|| A( e ) /||u|| A( e ) and of (line 3 and 4) the strain component ε xx (in 10 −2 ) for the original, uniform mesh and the 5th soft coarsened mesh with Dirichlet, periodic and Neumann coupling

Mesh coarsening and micro discretization error
The original, uniform micro mesh of the Diamond/SiC microstructure is discretized with 1098 × 1098 elements. The number of coarsening steps is again restricted to 5 steps. The diagram in Fig. 17 reveals that, similar to the previous examples, the estimated error increases rather moderately during the first coarsening steps, but strongly in the last ones. For an explanation for the dependency of errors on coupling conditions and their order in magnitude see Sect. 4.3.2. Table 7 underpins the graphical results by explicit figures; again, the mesh coarsening reduces the number of unknowns after five coarsening steps to less than 10% of the original, uniform mesh. The estimated errors exhibit -similar to the previous examples-a modest increase for the first coarsened meshes and a strong increase afterwards.
The ratio of deactivated degrees of freedom is 2.8% for the first coarsened mesh and increases to 18.8% for the fifth coarsened mesh.
The estimated errors in Table 7 underpin that, similar to the previous examples, the errors show modest increase for the first coarsened meshes and a strong increase afterwards.

Distributions of error and micro strain
Since the first examples have shown that the last coarsening steps do not pay off, the focus is on the results for the third coarsened mesh in the comparison with the original, uniform mesh for different coupling conditions. The error distributions of the coarsened meshes in Table 8 exhibit quantitative agreement in the maximal error figures with the uniform, fine discretization. Only in the interior of phases having undergone a considerable mesh coarsening the errors are somewhat increased.
The normal strain component ε x x of the diamond/SiC microstructure on one micro domain is displayed in Table 8 for the different micro-macro coupling conditions for the original, uniform mesh and for the third soft coarsened mesh. For each of the coupling conditions a very good agreement of strain distributions on the micro domain is achieved, although the coarsened micro mesh (line 4) contain less than 8% of the degrees of freedom of its original, uniform counterpart.

How micro discretization errors influence macro scale quantities
Micro discretizations do not only influence the accuracy on the micro level but also have an impact on the macroscopic stiffness through the homogenized elasticity tensor. In order to investigate and quantify this influence, two different, but related macro scale measures are employed; the maximum nodal displacement u max and the coefficients of the homogenized elasticity tensor A 0 , for the computation of the latter see [1,13]. The results indicate that for all coupling conditions the maximum nodal displacement is hardly influenced by mesh coarsening on the micro scale; even for the fifth coarsening step showing large micro errors the macro displacement error is below 0.015%.
The normed entries of the homogenized elasticity tensor displayed in Fig. 19 indicate that mesh coarsening does hardly influence the macroscopic stiffness. Even for the fifth coarsening step the maximal error of the homogenized coefficients is in the range of 0.02%. Notice that the very good agreement of the percentaged errors in A 0 and in u max hints at the fact that both the L 2 -norm of the macro displacement error and the error in the coefficients of the homogenized elasticity tensor converge in O(h 2q ) see [2,13,17].

Conclusions
The application of the mesh coarsening to a real world microstructure proves that an efficient mesh coarsening inside of the single phases is possible without significantly falsifying the results. Furthermore, the example has demonstrated that different discretizations on opposite edges of the RAE are not of harm for accuracy. Macro level quantities indicate that the considered mesh coarsening on the micro scale has virtually no impact on the macro level.   Table 9 reveals that the first two coarsening steps already reduce the number of unknowns by 87%, while the last two steps to hardly more than 4%. Results for the fifth coarsened mesh exhibit an error increase of 21% while the number of unknowns is reduced by more than 90%. Figure 21 displays the estimated micro error for mesh coarsening and the computation time to solve the micro problem. Similar to the previous examples the error shows only minor increase during the first coarsening steps and increases stronger afterwards. As a consequence, the computation time is strongly reduced in the first two coarsening steps to 45% and 18% respectively, while additional coarsening steps ren-der only a minor speedup. The computation time of the fifth coarsened mesh is 11% of the original computation time.

Distributions of error and micro strain
Based on the results of previous examples the analysis of the micro discretization errors has a focus on the third coarsened mesh and and its comparison with the original, uniform mesh for PBC. The error distribution of the coarsened mesh in Fig. 22 agrees very well with that one of the uniform, fine mesh; in the phase interior the larger elements from mesh coarsening induce somewhat larger errors. The same good agreement is observed for the normal strain component ε x x as displayed in Fig. 22. Again, the element size mismatch on opposite RAE edges has no visible influence on accuracy.

Conclusions
The application of the mesh coarsening to a microstructure with more than two different phases leads to similar results as the mesh coarsening for microstructures with only two different phases.

Summary and conclusions
The aim of the present work was the numerical analysis of quadtree-based adaptive mesh coarsening of pixelized microstructures, which is used as a preprocessor for micro meshes in two-scale finite element frameworks for homogenization such as FE-HMM and FE 2 . The basic assumption turned out to be justified in all examples, that interfaces are the mechanical hotspots in multi-material systems (here due to a stiffness mismatch; in a thermomechanical context [52] due to a mismatch in the temperature expansion coefficients, etc.) and therefore shall preserve high spatial resolution, while the phase interior with rather weakly varying field properties can undergo mesh coarsening.
Two different coarsening algorithms were investigated, one with a steep gradient of element size at phase boundaries, a second, modified one with a smaller gradient. To put the improved efficiency by mesh coarsening into relation to the accuracy loss, the corresponding micro discretization error was evaluated.
For that aim the well-known superconvergent patch recovery (SPR) was modified in its stress recovery to fit the requirements of microstructures with multiple material phases. To avoid the numerical effort of the SPR, a scheme based on an intra-element patch for the recovery of nodal stress and strain along with averaging of contributions from elements of the same phase was proposed and compared to the modified SPR.
The simulation results for various examples strongly support the following conclusions.
(i) The gradient of element size at a phase boundary should not be chosen too steep. The example in Sect. 4.2 has shown that for a hard and a soft coarsened mesh of virtually the same number of degrees of freedom, the soft coarsened one shows better agreement with the uniform, fine discretization in the distributions of error and strain. (ii) There is no need to execute expensive calculations of reference solutions on superfine meshes. As the example from Sect. 4.1 has shown, both error estimation schemes -modified SPR and averaging technique-lead to sufficiently accurate estimations of the true micro error. (iii) It does not pay off to execute as many coarsening steps as possible. All examples from Sect. 4 have shown, that only the first coarsening steps lead to major reductions of the degrees of freedom without a strong increase in the micro discretization error. The number of favorable coarsening steps was two or three for our examples with a number of degrees of freedom in between 1.6 and 3.3 million for the initial, uniform meshes. (iv) A microstructure with high regularity -meaning no sharp corners-has shown an outstanding relation between reduction of degrees of freedom and increase of micro discretization error. For the circular inclusion from Sect. 4.3 a reduction in the number of degrees of freedom to 8% of the original figure even lead to a micro discretization error which is increased by 8%. To achieve such outstanding results for less regular microstructures, the coarsening algorithms may be further developed, for example by keeping the gradient in element size even smaller at sharp corners. (v) The pixelized, hence non-smooth representation of smooth interfaces is not necessarily of harm for the regularity of the BVP. For the problem of a circular inclusion in a matrix the comparison with a smooth, standard discretization with elements oriented to the interface has shown the same convergence order for both cases. In general, however, re-entrant corners inherent to pixelized meshes can cause order reduction depending on relative pixel size and stiffness mismatch. (vi) Quadtree-based adaptive remeshing of pixelized microstructures has turned out to be an appropriate preprocessor for efficient finite element analyses. The examples from Sect. 4 have shown that already two coarsening steps are sufficient to reduce the number of degrees of freedom to considerably more than 80%, rather about 90%. The increase of the corresponding discretization error compared to the initial uniform mesh is not more than about 15%. (vii) While the micro scale exhibits requirements for accuracy and efficiency in its own right, the impact of micro discretization errors on the macro scale in a computational homogenization framework was analyzed as well. Here, the effective macroscopic response in terms of coefficients of the homogenized elasticity tensor and macroscopic displacements has turned out to be virtually insensitive even to considerable micro errors induced by mesh coarsening. The results even foster the considerable efficiency benefits through adaptive mesh coarsening in view of the vanishing accuracy losses.
(viii) All the above results are true for different conditions applied to RAE boundaries, for Dirichlet, periodic and Neumann coupling. For their errors it holds, the softer the condition, the larger the errors.

Concluding remarks
A comprehensive accuracy analysis must cover all involved aspects, from the image resolution along with possibly blurred interfaces, over segmentation and the assumed, corrugated or sharp characteristics of interface, to the type of discretization and the solution method in use.
In the present work we did not consider the question which initial pixel size is favorable in accuracy while feasible in numerical costs cf. [37,47]. A relatively coarse initial mesh will introduce a modeling error, where ensuing adaptive coarse graining will add a discretization error. The finite cell method shifts to some extent the efforts from meshing to numerical quadrature, but notwithstanding must make an assumption of the interface shape [10,11,38,48,49].
In the present two-scale FEM for homogenization, errorcontrolled, microstructure-adapted mesh coarsening is the key to enable the usage of highly resolved microstructure representations as a starting point with the benefit of considerably reduced computational costs for minor accu-racy losses. The extension to 3D mesh coarsening of voxel data is expected not only to intensify the efficiency benefits of the 2D scenario. It is expected that it will enable microstructure analyses in high resolution at all, which are otherwise, for the uniform, voxelized discretization, prevented by memory-limitations, if no parallelized domaindecomposition methods are used.