1 Introduction
Topology optimization (TO) can be described as an approach that optimally distributes material in a specified domain under a set of constraints, such that the performance function of the structure achieves a maximum [1]. In the past two decades, TO has widely been used in various academic and industrial disciplines. For a survey on the latest developments in TO as well as its recent applications, see the review papers by Sigmund and Maute [2], van Dijk et al. [3], and Deaton and Grandhi [4].
Typically, in popular densitybased TO, the domain is discretized into a finite set of elements and a density value is associated with every finite element [1]. The density of an element indicates the volume fraction of that element filled with a certain amount of material, and can vary from 0 (void) to 1 (solid). These density values are optimized during the course of optimization. Since in traditional approaches, density is assumed to be constant inside an element, a large number of finite elements as well as associated design variables are needed to obtain a well defined design with the desired structural features and boundary resolution, especially for threedimensional (3D) problems [5]. The computational costs associated with TO are mainly determined by the used finite element analysis (FEA) and associated sensitivity analysis, which limits the number of elements and consequently the design resolution.
With the growing popularity of TO, a clear need exists for improved methods that can deliver high quality results at the lowest computational cost. Various approaches have been proposed in the past to reduce the computational costs associated with solving largescale TO problems [6, 7, 8, 9, 10, 11, 12]. These focused mainly on improving the efficiency of solving the FEA systems of equations. Another possibility that has been explored in the existing literature is to modify the way the FEA system is defined in the first place through the use of adaptive FEA formulations. Popular adaptive FEA approaches are refinement and refinement [13]. However, the standard formulations for these methods use FEA based error criteria for adaptation of the mesh. These by themselves are not well suited for TO, as they do not take the need for refinement based on design considerations into account [14]. In the final designs obtained from TO, it is desirable that the material distribution is clearly defined. Thus, the refinement criterion used in TO should depend on the material distribution as well.
Maute and Ramm [15] proposed an adaptive mesh refinement (AMR) approach which involved optimizing the topology of the design followed by approximating the boundaries using cubic or Beźier splines. After every cycle of TO, shape optimization was performed followed by remeshing of the solid domain. The whole process was repeated over a series of cycles and the new mesh generated at the end of each cycle was used as the domain for the TO problem of the next cycle. Van Keulen and Hinton [16]
for the first time combined the TO with an FEA error based refinement strategy. The recovery of material, in their approach, was controlled by the stress level in the adjacent elements and mesh densities were determined using (a) the standard ZienkiewiczZhu error estimator and (b) the shortest distance to the materialvoid boundary. Both these approaches involved remeshing the whole domain at the end of each cycle, which was computationally expensive.
Costa and Alves [17] presented an AMR strategy which involved refining only the solid material region. For TO problems, intermediate densities are found to be prevalent near the boundaries. On the assumption that refinement of these regions can reduce the intermediate densities, Stainko [18] proposed to refine the region only around the materialvoid boundary. Bruggi and Verani [14] progressed in the direction of the work proposed by [16]
, and proposed a goalbased AMR strategy that properly guides the progression of refinement and coarsening in TO. For refinement/coarsening, a dualweighted residual based FEA indicator as well as a heuristic densitygradient based indicator were used. While most of these methods helped to achieve the required
adaptivity in TO, the fixed choice of density values for refinement at every cycle of TO led to excessive numbers of elements getting refined, thereby leading to undesired increase in computational costs. Gupta et al. [19] proposed a heuristic scheme to control the refinement/coarsening bounds at every cycle of TO. The proposed scheme was combined with refinement and very clear material descriptions with low gray regions were obtained. Other adaptive formulations involving refinement or a similar approach include adaptive refinement of polygonal elements [20, 21], combining a continuous density field representation with adaptive mesh refinement [22] and efficient TO based on adaptive quadtree structures [23].Another possible way to reduce FEA costs is the adaptive refinement, as stated earlier, where the mesh topology remains the same. Additionally, for smooth problems, the accuracy of refinement is dramatically higher than that of refinement for the same computational costs [13]. Increasing the polynomial order of the shape functions gives an exponential rate of convergence. Other advantages of refinement are its robustness against locking effects and high aspect ratios [24]. However, due to the fact that the conventional TO approaches assume an elementwiseconstant density distribution, using higherorder shape functions inside a finite element is not an efficient approach. Although it reduces the FEA error to some extent, it cannot improve the material definition within the element.
The recently proposed Finite Cell Method (FCM) offers new perspectives to overcome this limitation [25]. FCM is an FEbased modeling approach where the analysis mesh is decoupled from the material distribution domain and higher order shape functions are used [24]. This approach can handle a materialvoid boundary within an element through the use of appropriate integration schemes. Recently, a similar approach was proposed by Nguyen et al. [26] for TO, termed as multiresolution topology optimization (MTO), where the analysis and design meshes are decoupled. Here, design mesh denotes the distribution of the design points which are used to generate the material distribution. The density values associated with these points serve as optimization parameters for TO. In MTO, a coarse analysis mesh was used and inside every finite element, a large number of design points were positioned. This allowed a high resolution density distribution inside every finite element, unlike an elementwiseconstant density distribution as in standard TO approaches. In spite of using low order shape functions and coarse elements, the method is still capable of generating high resolution structures, albeit with reduced analysis accuracy. To increase this accuracy, recently a version of MTO has been proposed, where the potential of higher order polynomial shape functions has been investigated in the context of MTO [27]. Other approaches based on a similar concept were further presented in [28, 29]. Note that in [26] and other research papers thereafter, the term ‘multiresolution’ refers to allowing the possibility for multiple different design resolutions for the same choice of analysis resolution. In line with these works, we also refer to our formulation as an MTO approach.
It is important to note that although the design and analysis meshes can be decoupled, the iterative updates of the design variables in TO are based on the analysis results. In a recent study, we showed that for a given finite element mesh and polynomial order of FE shape functions, there exists an upper bound on the number of design variables that can be used in TO [30]. A density resolution beyond this threshold cannot be fully captured by the FEA and can lead to issues such as nonuniqueness. For certain cases, it can also lead to poorly performing designs. Thus, when using large numbers of design points inside an element, both for analysis accuracy as well as wellposedness of the TO problem, higher order shape functions and corresponding numerical integration schemes need to be chosen.
Parvizian et al. [25] proposed a TO strategy based on FCM where a coarse analysis mesh with high order shape functions as well as a high order numerical integration scheme is used. Although expected to give more reliable results, FCMbased TO may not necessarily satisfy the bounds proposed in [30], which implies it might still be prone to numerical issues. Groen et al. [31] presented results related to rigorous numerical investigations of FCMbased TO. Their observations show close resemblance with those in [30]. Also, the authors showed that using FCMbased TO, remarkable speedups of more than 3 and 60folds for 2D and 3D problems, respectively, could be obtained over the traditional TO approach. However, for certain configurations of FCMbased TO, it is possible that the design consists of ‘QRpatterns’, comprising disconnected or loosely connected material parts which cannot be correctly modeled by the employed modeling scheme [32]. Use of density filtering with a sufficient filter radius was found to suppress the QRpattern artifacts [27, 30, 31], but has the undesired consequence of reducing the design resolution. Applying refinement was also found to reduce the issue, but rapidly raises the computational cost.
Hereafter, we use the term MTO to refer to all the TO approaches (including FCMbased TO) where the design and analysis discretizations are decoupled. The goal of MTO approaches is to obtain high resolution, high quality designs at low analysis costs. Possible ways to increase resolution versus cost could include using a finely discretized density mesh, reducing the filter size, using shape functions of low polynomial order to describe the state field, etc. However, each of these approaches has certain limitations which can adversely affect the analysis accuracy. Using too many density cells and low polynomial order shape functions can lead to nonuniqueness in the design field and result in numerical instability [30]. Reducing the filter size can lead to the formation of QRpatterns, which are numerical artefacts and can affect the model accuracy [31, 32]. Using higher order shape functions can circumvent these problems, however, the analysis related costs are significantly increased. Due to this, the advantage of MTO over the traditional TO approach could be lost. In an MTO setting, this requires considering adaptivity both of the analysis and the design, which thus far has not been explored.
In this work, we present an adaptive MTO approach that enables a better balance between resolution and computational costs. Local adaptation is applied to both the analysis and the design description, which allows computational effort to be concentrated in regions of interest. Moreover, the adaptivity allows rigorous prevention of QRpattern artefacts. We coin the term ‘adaptivity’, an adaptive multiresolution TO scheme where both the design resolution and FE polynomial order can be locally adapted based on certain refinement/coarsening criteria. Here, the symbol ‘’ should not be confused with the one in  adaptivity, where it refers to domain decomposition and mesh overlaying [33]. It is assumed that computational costs are the limiting factor, and that the manufacturingimposed length scale is at or below the smallest lengthscale that can be reached by the adaptive TO process. Our approach can obtain high resolution representations of the material field at significantly lower computational costs compared to nonadaptive MTO approaches. At the same time, by jointly adapting design and FE discretization, we ensure that the bounds proposed in [30] are satisfied and instability issues are avoided. For refinement/coarsening after every TO cycle, analysis error, correctness of the design as well as the error associated with QRpatterns are used. For this purpose, we also propose a novel indicator. Various numerical tests are conducted to analyze the capabilities of the method as well as its robustness. The scope of this paper is restricted to linear elastostatic problems and the material is assumed to be isotropic, however, the method is expected to be applicable to a wider range of problems.
In the following section, theory of multiresolution TO is presented followed by discussions related to choice of design distribution, polynomial orders and numerical integration schemes. Section 3 subsequently presents the theory and formulation for the proposed adaptivity approach. The applicability of this method is presented on a set of numerical examples (Section 4), and discussion and related conclusions are stated in Section 5 and 6, respectively.
2 Multiresolution Topology Optimization
2.1 Domain and variable definitions
In this work, we propose an adaptive MTO formulation based on selective refinement/coarsening of the design as well as analysis domains. First a conceptual description is provided, whereas the mathematical formulation follows in Section 2.2. The proposed approach uses three meshes: design mesh, background mesh (comprising density cells) and analysis mesh. The analysis mesh is used to determine the solution of the physics at hand (e.g. displacement field) and the design mesh represents the distribution of design points in the domain. For simplicity, we use a structured mesh setting, as often used in topology optimization. In an adaptive setting, the analysis resolution and distribution of design points in the domain can be nonuniform. The background mesh is added to provide a convenient link between the analysis and design meshes. More details related to the role of the background mesh follow later in this section.
For practical implementation, we introduce the notion of MTO elements. An MTO element comprises a finite element, a set of design points and an overlapping background element comprising a regular grid of density cells. They all occupy the same spatial domain, and this ordered arrangement is chosen to simplify implementation in an existing FE framework. For example, Fig. 1 shows the schematic representation of a Q2/d8 MTO element using a Q2 (bilinear quadrilateral) finite element and consisting of 8 design points distributed nonuniformly in the domain. The overlapping background element comprises density cells. A density design variable is associated with each design point. During optimization, these density variables are updated at every iteration based on the response functions and the corresponding design sensitivities.
To generate suitably uniform distributions of design points within an element for any number of design variables, a variant of the
means clustering method is used [34, 35]. This approach divides the design domain into segments (clusters) with roughly equal areas. The design points are assumed to be located at the centroids of these clusters. For selfcontainment, the details of the method are discussed in Appendix A. We use this approach to obtain an approximately uniform distribution of any given number of design points in the MTO element domain. The achievable resolution limit of the design depends on the spacing between the design points. For a given number of design points and without a priori knowledge of the optimal design, a uniform distribution allows the best possible resolution. Note here that the proposed adaptive MTO approach is independent of the choice of methodology for the distribution of design points, and any other method to distribute points in a domain can be applied, including a set of predefined patterns.The aligned background mesh consists of a uniform grid of equallysized density cells in the whole domain, such that a certain number of these cells overlap with every finite element. For these density cells, the respective finite element is referred as the parent analysis cell. For example, in Fig. 1, density cells overlap with the parent Q2 finite element (analysis cell). The density is defined at the centroid of every density cell and is assumed to be constant inside it. This density is obtained from the design mesh through a localized projection operation.
The density inside any density cell of the background mesh is calculated using projection (as shown in Fig. 1, defined in detail in Section 2.2), and only those design points are used which lie within the same MTO element. The role of the localized projection is to define density values in all the density cells of the respective MTO element. The projection is restricted to the considered MTO element for two reasons: (i) to minimize the associated loss in design resolution of MTO elements adjacent to other MTO elements with fewer design points and (ii) to enable elementlevel implementation. While choosing the local projection radius , it needs to be ensured that the density inside each density cell can be defined. The mathematical details related to choosing this projection radius are provided in Section 2.2. An example is presented in Fig. 2, which shows a domain of MTO elements, each comprising a Q1 finite element and density cells. As can be seen, the distribution of design points can be nonuniform. The four MTO elements from topleft to bottomright consist of 4, 9, 3, and 7 design points, respectively. In the bottomright MTO element shown in Fig. 2, a partial projection circle can be seen, which is due to the fact that the projection is restricted to within this MTO element. Mathematical details related to projection are provided in Section 2.2.
The stiffness matrix for every MTO element is obtained by numerical integration using a Gaussian quadrature scheme. For this purpose, the stiffness matrix contribution at the integration point needs to be known, which in turn requires knowing the density value at that point. This density value, referred further as ‘projected density’, is obtained through a projection on the background mesh, denoted by (Fig. 1). Fig. 3 illustrates how these density values are computed. It shows a mesh of MTO elements, comprising Q1 finite elements and the corresponding background domain with density cells per element. Here, ‘Q1’ refers to quadrilateral finite elements with shape functions of polynomial order 1. Similar to the approach described in [26, 27, 31], the projected densities are computed using a distanceweighted projection of design densities found in the neighborhood of a certain radius over the background mesh. In this work, density filtering is used for the projection [36].
The use of the background mesh facilitates adaptivity, i.e. the use of different numbers of design points in adjacent elements. In the absence of the background mesh, the nonuniform design field when directly projected on the analysis mesh, can lead to irregular boundary features which are not desired. The design variables are not directly linked to the density cells of the background mesh, because it would not allow an adaptive formulation anymore. Moreover, such a formulation would significantly increase the number of design variables and would lead to nonuniqueness related issues [30]. The background mesh provides the flexibility of having a reference discretization independent of the number of design variables. Moreover, it simplifies the numerical integration required for the stiffness matrix.
2.2 Mathematical formulation
In this paper, the applicability of a adaptive MTO approach is demonstrated on mechanical problems of two different types: minimum compliance and compliant mechanism.
For the chosen problems, the problem statement for TO can be expressed as
(1)  
where, denotes the objective functional, and , and
denote the global stiffness matrix, displacement vector and load vector, respectively. The vector
is chosen based on the type of problem and will be discussed in Section 4.1. The volume constraint restricts the total volume fraction of the given material to be less than certain predefined volume .Next, the details related to various steps associated with the proposed multiresolution modeling approach are described. The matrix in Eq. 2.2 is obtained from the global assembly of the element stiffness matrices , which can be expressed as
(2) 
where and denote the straindisplacement matrix and constitutive matrix, respectively, and is the number of integration points. More details related to the choice of numerical integration are discussed in Appendix B. The subscript refers to the integration point and denotes the respective integration weight. The construction of the
matrix depends on the choice of the material interpolation model as well as the material itself. In this work, solid isotropic material interpolation (SIMP) model
[1] is used such that(3) 
where is the Young’s modulus of the solid material and is a very small value (typically ) used to avoid singularity of the system stiffness matrix. Also, denotes the density calculated at the integration point, is the penalization power and denotes constitutive matrix normalized by the Young’s modulus.
The densities at the integration points are calculated by projecting density values from the density cells in the background mesh (Fig 3). For this purpose, we employ a linear projection approach for based on the density filtering method which is widely used in TO [36]. Mathematically, it can be stated as
(4) 
where refers to density values for the cells contained in the background mesh with their centers lying within a distance from the corresponding integration point (Fig. 3), and their number is denoted by . Here, terms reduce linearly with distance from the integration point, i.e.,
(5) 
where denotes the Euclidean distance operator.
As stated in Section 2.1, the background mesh densities are calculated using the projection from the design mesh to the background mesh. For the MTO element, the density of the density cell is given as
(6) 
where, refers to the density value associated with the design point in the design domain contained within the MTO element, and lying within a distance from the centroid of its density cell. The number of such design points is denoted by , and is the radius of the projection for the element (Fig. 2). Here, is defined as
(7) 
As stated earlier, the projection radius needs to be chosen such that it is as small as possible, however, large enough to define densities for all the density cells that correspond to the respective element. Here, we define it as
(8) 
where denotes problem dimension, and is the edgelength of the MTO element. The operator denotes ceiling function which rounds the contained floatingpoint number to the nearest greater integer value. The term refers to edgelength of the density cells. Next, to obtain a projection length slightly larger than the diagonal, we multiply by . Note that Eq. 8 has been obtained empirically through observations based on various design distributions obtained using the means clustering approach. For other approaches of choosing the locations of design points, where for any value of , the distance between the design points can be provided mathematically, it is possible that even lower values of work. Lower values of can help to further reduce the loss in design resolution caused due to the choice of localized projection , and this could be a potential direction for future research.
Fig. 3(a) shows an example of a cantilever beam subjected to a point load, which we will use to illustrate the MTO concept. The domain is discretized using finite elements. For each MTO element, 225 design points, distributed in a square grid of , are used to represent the design field. The polynomial order of the shape functions is chosen to be 10. The choice of shape functions is made in a way that the elementlevel uniqueness bounds defined in [30] are not violated. As per the uniqueness bound, the number of design points influencing any finite element cannot be greater than the number of deformation modes of that element, With equal to 10, the number of deformation modes is 239, which is greater than 225. With and equal to 10 and 225, respectively, the MTO elements are referred as Q10/d225 type elements. For this example, the projection radius is set to 0.13 times the elementlength, which is equivalent to the size of 2 density cells.
Fig. 3(b) shows the optimized design obtained using the outlined MTO configuration. Clearly, the employed MTO approach allows the definition of higher resolution material features on relatively coarser MTO elements. However, in Fig. 3(b), there are parts of the domain where even lowerorder elements and lower design resolution are sufficient. For example, there are totally void MTO elements, where even linear shape functions with only one design point can be used. Clearly, the computational time of the MTO approach can be reduced by exploiting this fact in an efficient way, and in the next section, we propose an approach to do this.
3 dpadaptivity
3.1 General description of the method
We present here a adaptive version of the MTO method which is capable of enhancing further the ratio between the design resolution and analysis cost compared to nonadaptive MTO. The proposed MTO method efficiently distributes the design variables and locally adapts (increases/decreases) the polynomial order of the shape functions. A threepart refinement criterion is defined to select the cells to be refined/coarsened. Note that although the term ‘refinement’ is more commonly used throughout this paper, we implicitly refer to coarsening (reducing the values of and ) as well. Here, ‘refined’ cells are those where additional design points are inserted, or the polynomial order of the shape functions is increased, or both. Similarly, ‘coarsened’ cells are the ones where the design resolution (number of design points) is reduced, or the analysis resolution (shape function order) is reduced, or both. With an adaptive formulation, fewer design variables as well as analysis nodes are used, which provides a computational advantage over the conventional MTO method.
At the start of adaptive MTO, a cycle of TO is performed, using a certain initial uniform design and FEdiscretization. A ‘TO cycle’ refers to the entire process from starting with an initial design and optimizing it over a number of iterations (or up to a certain stopping threshold) to reaching an improved design. During a TO cycle, the shape function order and design points of all elements remain fixed. In the optimized design, refinement and coarsening zones are subsequently identified based on an integrated criterion comprising an analysis errorbased indicator, a densitybased indicator, and a QRbased indicator. Here, QRerror refers to the error due to the incapability of the chosen shape function in modeling the displacement field arising from a highresolution density representation allowed within that element [32]. More details related to these indicators are discussed in Section 3.2.
All steps from analyzing the design for refinement to updating the and values for the whole domain, constitute one cycle of adaptivity. The general structure of a adaptive MTO cycle is as follows:

Perform optimization of an MTO problem with fixed and values.

Adapt values based on analysis error indicator.

Adapt and values based on densitybased criterion.

Update values to reduce QRerrors in every element.
With the new adapted mesh, the next cycle of TO is performed. Section 3.3 below describes each of the above steps in detail.
3.2 Refinement criteria
In this section, the details related to the three indicators used in our refinement criterion are provided. As stated earlier, although the term ‘refinement’ is frequently used, we implicitly refer to ‘coarsening’ as well in our adaptive approach. Note that although here certain choices have been made for the refinement indicators, the adaptive scheme in itself is not dependent on the choice of refinement indicator, and can be coupled with other appropriate indicators as well.
3.2.1 Analysisbased refinement indicator
For the purpose of analyzing the modeling related error, the Kelly error estimator has been used [37]. This error indicator analyzes the jump in the gradient of the solution across any face (edge in 2D) of adjacent elements. The error for any element is calculated in a relative sense by integrating the error in the gradient jump across all faces of the respective element. Based on the relative error estimate, only a certain fraction of the MTO elements is selected for updating the orders of the polynomials (). This error estimator can also be understood as a gradient recovery estimator, for details on this aspect, see [38].
There are two reasons to choose the Kelly error estimator instead of more sophisticated recent approaches, e.g., goaloriented error estimators [14, 39]. The analysis error comprises primarily of two components: element residual and edge residual [39]. Element residual refers to the error in approximating the gradient field within the element, and edge residual denotes the jumps in gradient across the element edges. The element residual is being taken into account through the QRerror analysis. Thus, the analysis indicator needs to only look at the edge residual term. Moreover, our approach requires only a relative error estimate and not the exact error itself. The use of Kelly error estimator suffices both these requirements. Also, this error estimator is simple to implement and the associated computational costs are negligible.
For the purpose of ranking the elements for adaptivity based on the Kelly error estimator, the analysis residual error vector needs to be defined. For the MTO element, can be computed as:
(9) 
where, refers to a face (edge in 2D) of the element and operator denotes the jump in the argument across face . Also, denotes the set of all faces of the element. The constant term is set to , where is the element diagonal and denotes the maximum among the polynomial degrees of the adjacent elements [40]. The residual errors are ranked, and the top 10% and bottom 5% of the elements are selected for increasing and decreasing the values, respectively.
For illustration purposes, we perform a partial adaptive MTO run on the problem shown in Fig. 3(a). Fig. 4(a) shows the optimized cantilever beam design obtained for this problem after one TO cycle. The design has been obtained on a mesh of Q2 finite elements with design points per element. The optimized design clearly shows typical artefacts (QRpatterns) of disconnected structural features. Fig. 4(b) shows the distribution of polynomial shape function orders obtained from adaptivity controlled by only the analysisbased refinement indicator. It is observed that coarsening (reduction in ) has mainly occurred in the void cells which are far from materialvoid boundaries. This is because the jumps in displacement gradients across the edges for these elements are zero. For refinement (increase in ), the elements at the boundary have been preferred.
3.2.2 Densitybased refinement indicator
The densitybased refinement indicator aims at adaptively choosing MTO elements for refinement/coarsening in way that over a number of cycles, the intermediate densities are reduced, and a crisp and highresolution boundary representation is obtained. For this purpose, the refinement indicator proposed in [19] is adapted for our problem and discussed here. This indicator chooses a certain element for refinement/coarsening based on the density value inside that element. For every cycle of adaptivity, refinement (coarsening) density intervals are defined and associated elements are flagged. We adopt this indicator to regulate the number of design points in each MTO element, based on spatial design information specified by the density values of the voxels of the background mesh. The way this indicator affects the number of design variables is discussed in Section 3.3, here we focus on the definition of the indicator itself.
Fig. 6 shows the refinement () and coarsening ( or ) intervals as a function of adaptive cycle. Unlike the other refinement indicators, here the refinement (coarsening) bounds are chosen not to remain constant. Rather, following [19], the range of density values to be chosen for every adaptive cycle increases. Based on the chosen stopping criterion used for every cycle of TO, it is possible that significant parts of the designs obtained during initial cycles consist of intermediate density values. In such scenarios, selecting all gray (intermediate density) elements for refinement can lead to excessive refinement during the initial cycles, which in turn leads to undesired increase in computational burden. Due to the adaptive nature of the refinement indicator proposed in [19], such problems can be avoided.
To start, the densitybased refinement indicator for the MTO element is set to 0. To update , we iterate over all the density cells of the MTO element and consider the sum of individual refinement or coarsening contributions of these cells. Let denote the number of density cells contained within the background mesh associated with the MTO element. Then is updated as follows:

Iterate over from 1 to :

Let the density of the voxel be denoted by .

if ,
set . 
if ,
set . 
if ,
set . 
if ,
set .

Here, the average density is defined using the expression . The variables , , and are the bounds used to characterize the refinement and coarsening zones as shown in Fig. 6, and are defined as follows:
(10)  
(11)  
(12)  
(13) 
Here, denotes the adaptive cycle index, and and are tuning parameters chosen here to be 0.2 and 0.8, respectively.
The tuning parameters and are independent of the index of the adaptive cycle. However, is sensitive to the rate at which the design converges. As stated earlier, our method assumes that the design has sufficiently converged at the end of every optimization cycle. For different problems as well as different mesh resolutions, the amount of gray region may vary at this point. For problems where the designs of initial cycles of the adaptive MTO process are significantly gray, lower values of are recommended. This allows the density range for refinement to expand slowly over a span of cycles. Similarly, for rapidly converging designs, a larger value of is more efficient. Automated adjustment of these parameters could be considered, however, it has not been used in this study.
Fig. 7 shows the shape function field and the design field obtained for the optimized cantilever beam design shown in Fig. 3(a). The shape function field (Fig. 6(a)) denotes the polynomial order of the shape functions used in every finite element. The design field (Fig. 6(b)) denotes the number of design points used in every analysis element. These distributions have been obtained based on adaptive refinement and coarsening controlled by only the densitybased refinement indicator. From Fig. 7, it is seen that the materialvoid boundaries where the intermediate densities are prominent, have primarily been refined. Coarsening occurs in void parts of the domain.
3.2.3 QRerror indicator
In an MTO scheme, it is possible that the employed shape functions cannot accurately model the displacement field arising due to the allowed high order density representations. As stated earlier, this error arising in an MTO setting due to inappropriate modeling is referred to as QRerror. A closedform condition to predict this QRerror is currently not known. Groen et al. [31] proposed a method to estimate the average error for the whole domain by determining a reference solution using a refined uniform mesh, and evaluating the obtained MTO solution against it. In the context of adaptivity, QRerrors must be quantified at element level. We have proposed a method in [32], where an approximation to the QRerror can be obtained for any element through a comparison with a reference solution obtained by local refinement. In this work, we use this costeffective local QRerror indicator proposed in [32]. Once a sufficiently converged design has been obtained from a TO cycle, the QRerror is determined by evaluating the effect of local refinement, as follows.
Let , and denote the element stiffness matrix, displacement solution and internal load vector for the MTO element. Here, denotes the polynomial degree of the shape functions used in this element. Let denote the displacement solution obtained for the element using shape functions of polynomial order . Note that will be obtained by solving the elementlevel system . Here, nodal load is formed by integrating the product of the interpolated original load field and the refined shape functions.
To obtain a unique solution for
, sufficient boundary conditions need to be imposed. Thus, degrees of freedom (DOFs) equal to the number of rigid body modes (3 for 2D) need to be fixed. For this purpose, the displacement solution at 3 DOFs of
is copied directly from for the DOFs which overlap, and the solution at the rest of the DOFs is obtained through solving the finite element system. Once has been obtained, the QRerror can be computed as(14) 
where refers to elementlevel strain energy for the finite element using shape functions of order . Thus, and have been used. This strainenergybased criterion (Eq. 14) has been found to work well for the cases shown in this paper.
Fig. 7(a) and 7(b) show an optimized design obtained after first cycle of MTO run, and the corresponding error distribution obtained using the QRerror indicator for the problem shown in Fig. 3(a). Since the elementlevel test for QRerror is very conservative, it predicts higher error values compared to the actual fullscale TO problem [32]. Thus, to avoid undesired excessive increase in the values of , we restrict the increment of by only 1 per adaptive cycle based on the QRerror test. Also, to avoid excessive spatial refinement per adaptive cycle, only the cells with error value larger than 0.9 are adaptively refined. The elements flagged for refinement are shown in Fig. 7(c). It is observed that the regions where the QRpatterns exist, have been flagged for refinement. Moreover, elements at the material boundaries, which are partially void or solid, also show high value of QRerror and are flagged.
An interesting observation in Fig. 7(b) is that the elements which are completely void or solid also show QRerror values in the range 0.30.5. Although significant, the QRerror values in this range are relatively smaller than other parts of the domain and these elements do not get flagged for refinement. The reason for substantial QRerror values in these regions is the use of low order shape functions. For low values of , the displacement solution for even a uniform density field may not be accurately modeled. When solving elementlevel FE problems with low shape function orders and , it is observed that the modeling accuracy significantly improves when is increased. Due to this, nonzero large values of are recorded in solid and void parts as well.
3.3 adaptivity algorithm
The different steps of adaptivity have briefly been introduced in Section 3.1. After treating the three indicators involved, here we discuss each of these steps in more details. Once a TO cycle has been completed, the optimized design is analyzed using the composite refinement criterion, and the following steps are carried out.

Once a cycle of TO run is completed, get the optimized design for adaptivity.

Perform adaptivity based on analysis error criterion.

Update values for the whole analysis mesh (discussed in Section 3.2.1), where is the analysis error indicator value for the MTO element.

Sort in ascending order such that a corresponding ordered set is obtained.

Set the refine/coarsen flag of the element to 1 for the first fraction of the MTO elements in , and , for the last fraction of the elements. Here, and 1 denote that the cell has been flagged for coarsening (decrease in value) and refinement (increase in value), respectively. For no refinement/coarsening, is set to 0.

Increase/decrease values based on flag .


Refine/coarsen and values based on densitybased refinement criterion.

Update values for the whole domain (discussed in Section 3.2.2), where is the densitybased refinement indicator value for the MTO element.

Sort in ascending order such that a corresponding ordered set is obtained.

Update values by iterating over from 1 to :

For the first fraction of the elements in , do:

if , set . This helps to identify that the current element has been checked for coarsening. Since cannot be lower than 1, no coarsening is performed.

if and , set .


For the last fraction of the elements in , do:

if or , set . This means that if the element has been coarsened or left untreated based on the analysis indicator above, then refine it.



Reduce the difference of values between adjacent elements to a maximum of 2 at this point. This is achieved by iterating through the whole domain () times, where and are the maximum and minimum values of in the domain. At every check, the correction is done by raising the lower value of .

Update the designfield ( values) by iterating over from 1 to :

if , set = 1. This situation occurs when , and the densitybased indicator flags the cell for further coarsening.

if , set equal to the elementlevel upper bound for the element (based on [30]). Thus, , where denotes the number of rigid body modes for that element.


Update the background mesh

Find maximum number of design variables per MTO element ().

Find first perfect square (cube in 3D) number () greater than .

Set the number of density cells per MTO element equal to .

Update projection radius for every MTO element (Eq. 7).



Update values to reduce the QRerror in every MTO element.

Iterate over from 1 to , do:

Calculate the QRerror for the cell (discussed in Section 3.2.3).

Update for the element, if QRerror is greater than a certain error tolerance .


The adaptive MTO cycle is complete once the domain has been adaptively refined based on the three indicators. With the new refined mesh, the next cycle of TO is performed.
4 Numerical tests
4.1 Definition of test problems
To demonstrate the applicability and effectiveness of adaptivity, two test problems of minimum compliance and one compliant mechanism problem are considered [31]. In this paper, only 2D problems are studied, whereas an extension to a 3D setting is a part of future work. Young’s modulus is set to 1 Nm, , and the SIMP penalization factor is set to 3. The domain in each case is discretized using an initial mesh of MTO elements, comprising quadrilateral finite elements with shape functions of polynomial order 2 and design points per element. The radius is set to 0.3, where is the edgelength of any MTO element in the mesh. As a stopping criterion for all the test cases used in this paper, the optimization process for the cycle is terminated when the change in objective value between two consecutive iterations is less than . Here, denotes the minimum required change in objective value between two consecutive iterations of the first MTO cycle, below which the optimization process terminates. For the subsequent cycles, the minimum required change in objective value is reduced by a factor of at every MTO cycle. The adaptive stopping criterion used here allows to control the extent of design convergence per cycle. For the numerical examples used in this paper, and are set to 0.04 and 0.6, respectively, and these values have been found to work well. Based on this, the first () and second () optimization cycles are terminated if the minimum changes in objective value are less than 0.04 and 0.024, respectively.
To validate the accuracy of the MTO modeling of the design, we use the method proposed in [31], where the obtained design is compared with a reference solution. For the reference solution, we discretize the domain using a highresolution traditional TO mesh with elementwise constant densities. In this paper, the reference mesh comprises 320 160 finite elements and the polynomial order of the involved shape functions is set to 3. With this mesh configuration, the resolution of the reference domain is equal to the highest density resolution that has been used in the MTO problem.
For the first test problem, compliance needs to be minimized for a Michell beam cantilever subjected to a point load (Fig. 3(a)). For this case, N and = 1 m. Three variants of this problem are used with maximum allowed material volume fractions set to 0.45, 0.2 and 0.1, to study the capability of the method in low volume fraction problems on coarse meshes. For the other problems used in this paper, only one volume constraint of 0.45 is considered. The second test problem is that of compliance minimization for a cantilever beam subjected to a distributed load (Fig. 8(a)), and it is ensured that the load is consistently distributed over the various cycles of adaptivity. Here, and m. The distributed load tends to generate a lot of fine structures locally, and the resultant design was earlier found to be prone to QR artefacts [31], which makes it an interesting problem. For both these problems, the objective functional of Eq. 2.2 with . The third case is a compliant mechanism problem where a force inverter needs to be designed, such that for a point load at one end, the displacement at the other end is maximized (Fig. 8(b)). Here, spring stiffnesses and are set to 1 Nm and 0.001 Nm, respectively. For the force inverter, in Eq. 2.2 is a vector of zeros with 1 contained at the DOF where needs to be maximized. Thus, . The flexure hinges that are formed in this compliant mechanism problem will have subelement resolution, and this aspect makes also this problem an interesting test for our method.
4.2 Results
Here, we discuss the results obtained for the three test problems using a adaptive MTO scheme. To provide an understanding of the computational advantage of the proposed method, a comparison of CPU times is performed for the designs obtained using the proposed method as well as those obtained using the conventional MTO scheme discussed in [31]. Groen et al. [31] have shown that by using the MTO approach, the computational time can already be reduced by factors of up to 2.9 and 32 for 2D and 3D problems, respectively, compared to the traditional TO approach. In this paper, we demonstrate the potential of adaptive MTO schemes for 2D problems, and for this purpose, we will compare its performance with the nonadaptive MTO scheme, implemented in the same framework and evaluated on the same computing hardware.
4.2.1 Compliance minimization for point load
Problem  Definition  Speedup  DOFs  
Minimum compliance  point load  0.45  4.5  0.98  0.98  22935  17262 
0.20  8.3  0.93  0.98  20056  15096  
0.10  10.0  1.03  0.96  19590  15186  
distributed load  0.45  4.6  0.98  1.0  22636  16932  
Compliant mechanism    0.45  6.2  1.01  1.0  23375  17516 
This case refers to a maximization problem, where a value higher than 1 denotes that the adaptive MTO approach  
performed better over the nonadaptive MTO scheme. 
Fig 10 shows two optimized cantilever designs obtained for the problem shown in Fig. 3(a). The first design (Fig. 9(a)) has been obtained using the traditional nonadaptive MTO scheme, and the other (Fig. 9(b)) by our adaptive approach. For the two cases, the maximum allowed material volume fraction is set to 0.45. Visually, the designs differ only slightly. Table 1 provides the details on various parameters related to MTO cases for the two optimized designs. The first remarkable observation regarding the adaptive MTO result is the reduced computational cost. Adding the adaptive framework to the existing MTO allows a reduction in computational cost by a factor of 4.5. This reduction in cost is mainly due to the reduced number of design variables and free DOFs used in the adaptive MTO case. While the uniformly refined mesh used in MTO comprises 51200 design points and 40400 free DOFs, only 22935 design points and 17262 free DOFs are used in the final (4) cycle of the adaptive MTO run, i.e. a reduction by over 50%. The free DOFs and number of design variables used in the earlier cycles are even lower (Table 2).
Another reason that accounts for the speedup is the reduced number of iterations required in the final cycle of the adaptive method under the same stopping criterion as used for the nonadaptive MTO method. The convergence of the TO process is significantly affected by the choice of the initial design [41]. In our approach, each preceding cycle, after refinement/coarsening, provides a high quality initial design for the next one. Since the design converges significantly in the first 3 cycles itself using less refined meshes, only 18 iterations are needed in the final cycle, while the nonadaptive MTO scheme uses a total of 56 iterations. Table 2 provides the details related to the adaptive MTO run for this case. It is observed that Cycles 1 and 2 use a higher number of iterations. However, since the number of design variables and free DOFs are lower during these cycles, the associated computational cost is not very high.
Cycle  DOFs  Iterations  

1  6560  12800  67  0.86 
2  7204  10646  34  0.97 
3  12818  16256  17  0.98 
4  17262  22935  18  0.98 
In terms of performance, the cantilever design obtained from the adaptive approach slightly outperforms the design obtained using nonadaptive MTO. The obtained performance ratio is equal to 0.98, where and denote the compliance objective values obtained using the proposed method and nonadaptive MTO, respectively. From Table 2, it is observed that the global solution accuracy , where and refer to the objective values reported using adaptive MTO and that evaluated using the reference mesh, respectively. Since solution accuracy is close to 1, it is implied that the final optimized design is correct and free from artefacts. Moreover, we see that with every cycle of refinement, the global solution accuracy has improved. Thus, the adaptive MTO method allows to obtain designs with a desired analysis accuracy.
Fig. 11 shows the distributions of shape function order and design points as well as the optimized designs for 4 cycles of the adaptive MTO run of this case. It can be seen that refinement mainly occurs near the edges of the structure, and coarsening occurs as desired in solid and void parts. The optimized design in Cycle 1 consists of disconnected features, which are primarily the QRpatterns arising from the limitations of low order polynomial shape functions in those parts of the design [32]. Over the next cycles, refinement occurs in those regions and the QRpatterns are eliminated. Since the design points are distributed in the domain using means clustering without symmetry constraints, the distribution of design points itself can be asymmetrical, which in Cycle 2 leads to an asymmetrical design. An example of such asymmetry can be observed in the optimized design of Cycle 2, which gradually disappears over the next cycle.
In general, TO problems involving lower volume fractions of material are more difficult in terms of convergence. Moreover, for problems involving low volume fractions of material, a significant part of the domain comprises voids, and in turn does not require a fine mesh resolution. Clearly, for such scenarios, adaptivity could be potentially beneficial. To investigate this, we study two additional cases of the point load cantilever beam involving lower values of .
Fig. 12 shows the optimized designs for using conventional MTO (Fig. 11(a)) and adaptive method (Fig. 11(b)), respectively. For , the computational time advantage has increased to a factor of 8.3. Also, it is seen that the design obtained using the nonadaptive MTO method differs significantly from the result of adaptivity. Moreover, in terms of performance, the design obtained using adaptivity is relatively less compliant. The ratio is equal to 0.93. The compliance accuracy of the design obtained using the proposed method is found to be 0.98.
As another test case for lower volume fractions, the point load cantilever problem is examined with . Fig. 13 shows the optimized designs for this volume fraction obtained using the conventional MTO method and adaptive MTO, respectively. It is observed that for this volume fraction, the relative reduction in computational cost is even higher. Compared to the conventional MTO, a speedup of 10 times is observed. The increase in speedup is mainly due to the reduced number of free DOFs and design points, and the lower number of iterations required for convergence compared to the nonadaptive MTO. For this case, it is observed that is 1.03, which implies that the design obtained using adaptivity is slightly inferior to that obtained using the nonadaptive version. The analysis accuracy is also slightly lower than in the previous cases, with .
An understanding on the convergence of the adaptive MTO process for can be obtained from Fig. 14. In the first cycle, the design distribution and shape function orders are uniform for the whole mesh. Similar to the case of , it is observed that QRpatterns are formed here as well, which are removed by refinement in later cycles. Compared to Fig. 11, it is observed that only a small part of the domain gets refined. Because of the low volume fraction of material used, a significant part of the domain comprises mainly of void regions, which do not require refinement. For the nonadaptive as well as the adaptive versions of MTO, it is observed that the convergence of the optimization problem slows down significantly when very low material volume fractions are used. For example, for the same error tolerance, the number of iterations required in the final cycle of adaptive method for 0.45 and 0.10 are 18 and 82, respectively. Our observations on the effect of material volume fraction on the convergence of TO process align with the results reported in [42], where similar results have been obtained over a set of numerical experiments.
4.2.2 Compliance minimization for distributed load
For the cantilever beam subjected to a distributed load (Fig. 8(a)), is set to 0.45. Fig. 15 shows the optimized designs obtained using a uniform MTO mesh (Fig. 14(a)) and the adaptive approach (Fig. 14(b)). The information on the two runs is listed in Table 1. As in the case of the point load cantilever, the designs obtained using the nonadaptive and adaptive variants of MTO are very similar. In terms of performance, a speedup of 4.6 times is observed, and the accuracy of the obtained solution is close to 1. The obtained value is 0.98, which implies that the adaptive MTO found a slightly stiffer design.
For both the designs, there exists a small region near the top right boundary which comprises intermediate densities and is not improved even with refinement. With adaptive MTO, this region is more prominent. Among the possible reasons, one explanation could be that the distributed load applied on the upper boundary of the domain requires support material in those parts. In the absence of material near the upper boundary, the load point can get disconnected, which leads to a high overall compliance value for the structure. We observe that the optimizer is not inclined towards adding much solid material in these parts of the domain. Due to this, gray regions are formed, representing fine structural features beyond the design resolution. These intermediate densities can be suppressed by the use of methods such as modified Heaviside projection as has been demonstrated in [31], or simply by adding a solid nondesign region at the top surface.
Using a stronger penalization on the intermediate densities at the later cycles of MTO has also been found to help in reducing the gray areas. Fig. 16 shows two optimized designs for this cantilever problem obtained using adaptive penalization schemes. For the first case (Fig. 15(a)), the initial value of is 3 and it is increased by 1 at every cycle. For the second case (Fig. 15(b)), the increment is by 2 at every cycle. It is observed that with stronger penalization on the intermediate densities, the gray regions are significantly reduced.
To obtain an understanding on how the design evolves over 4 cycles of adaptive refinement, see Fig. 17. Due to the low order of the shape function used in Cycle 1, QRpatterns are observed here. Similar to the previous cases, adaptive refinement in the affected regions helps to remove these artefacts. For Cycle 4, only 16 iterations are needed when using the adaptive method, while the conventional MTO method uses 54 iterations in total. Also, the number of design points and DOFs used in the last cycle of the adaptive MTO are lower than in the conventional MTO method. Together, these two factors make the adaptive MTO method 4.6 times faster in this case.
4.2.3 Force inverter compliant mechanism
To demonstrate the applicability of adaptivity on topology optimization of compliant mechanisms, it is applied to the force inverter problem shown in Fig. 8(b). The allowed volume fraction is set to 0.45 and the goal of the problem is to distribute the material in a way that the displacement is maximized. Fig. 18 shows the optimized designs obtained using conventional MTO (Fig. 17(a)) and the adaptive method (Fig. 17(b)). As in the previous cases, the two designs are very similar. Details related to the MTO runs are reported in Table 1. It is observed that the objective ratio is 1.01. Since this is a maximization problem, a value of higher than 1 denotes that the design obtained using adaptive MTO performs better. is equal to 1.0, which means that the solution is as accurate as the reference solution.
Fig. 19 shows the distribution of design points and shape function orders, as well as the optimized designs for each cycle of adaptivity. Similar to the other cases discussed in this paper, QRpatterns are observed in the results of the first cycle. Nevertheless, the overall material distribution after Cycle 1 already corresponds to the final solution. The QRpatterns eventually disappear in the subsequent cycles due to adaptive refinement of the domain. Refinement primarily occurs in regions where intermediate densities are prominent, and coarsening mainly occurs in the void and solid parts of the domain.
5 Discussions
The primary goal of using an MTO scheme is to obtain a highresolution design at a relatively low computational cost. MTO decouples the design and analysis meshes in way that even for the choice of a coarse analysis mesh, a highresolution density field can be obtained. The potential of MTO has already been demonstrated in [31, 27]. However, there are a few aspects of MTO (e.g. computational cost, QRpatterns) where scope of improvement existed. The adaptive approach presented in this paper addresses these aspects and further enhances the capability of the MTO method.
This paper has mainly been focused on presenting the rationale and detailed formulation of the method. To demonstrate the applicability of adaptive MTO, 2D mechanical test problems have been considered in this study. Intended future work includes exploring the application of the proposed method on problems involving other physics as well as in 3D settings. In [31], it has been shown that MTO can bring a speedup of up to 32 folds over the traditional TO scheme. The improvement in 3D is significantly higher than that observed in 2D. As adaptive MTO reduces the DOFs compared to the conventional MTO method, it is certainly expected to pay off even more in 3D. To really understand the value of the adaptive approach for 3D problems, this hypothesis needs to be tested, and this is a part of our future work.
A preliminary investigation related to the application of adaptive MTO on linear conduction (thermal/electrical) problems with loads distributed throughout the domain, revealed that this approach could bring only limited improvements in speed (less than twofolds) for this problem class. The primary reason is that for this type of problems, the optimized design comprises fine features, dendritic in nature, which spread all across the domain. For example, Fig. 19(a) shows an optimized design obtained for a linear thermal conduction problem using the traditional TO approach. A mesh of elements was used and was set to 1.5 times the length of the element. The material volume fraction was set to 0.3. Details related to the definition of the problem can be found in [44]. It is seen that the optimized design has very few extended void areas, and most of the domain consists of fine material branches. Due to this, the majority of the domain gets refined at every adaptive cycle, which eventually reduces the relative advantage of adaptive MTO method over its nonadaptive variant.
Fig. 19(b) shows an optimized solar cell front metallization design obtained using the traditional TO approach on a mesh of finite elements and set to 1.5 times the element edge length [43]. This design has been obtained by solving a nonlinear electrical conduction problem, and only 45% of the domain is filled with material. For this case, it is seen that significant parts of the domain consists of void regions, which can be easily modeled with low values of and . Clearly, for such cases, the adaptive approach can be used to significantly reduce the associated computational costs. From the two examples of conduction problems discussed here, it is clear that adaptivity could certainly have a potential value for problems where designs feature extended void regions.
To demonstrate the concept of adaptivity, a composite indicator has been formulated in this paper. This indicator consists of an analysis error indicator, a densitybased indicator and a QRindicator. Although certain choices have been made for these indicators, the presented methodology itself is independent of these choices. Either of these indicators can be replaced with other alternatives that exist in the literature. For example, the Kelly estimator used as an analysis indicator in this work can be replaced with other analysisbased refinement indicators, e.g., goaloriented error indicator [45]. Such choices can provide a better control over the absolute error, accordingly helping to make a better choice of mesh resolution and solution accuracy. However, it is important that the tuning parameters associated with the chosen indicators are properly set so that issues related to excessive refinement are avoided. An addition to consider is a limit on, e.g., the increase in DOFs and/or design variables at a given adaptive cycle.
For the analysis indicator discussed in this paper, the top 10% and bottom 5% of the elements corresponding to are chosen for refinement and coarsening, respectively. There is no particular motivation to choose these cutoffs. For problems where the design domain has prominent regions with large jump across the element edges, it is recommended to allow more cells to be refined, so as to reduce the error in fewer cycles. For the densitybased indicator, both and are set to 1.0 for the current study. This ensures that all the elements with are refined and all elements with are coarsened. The reason to set these parameters to 1.0 is that the stopping criterion chosen in this paper allows the design to converge sufficiently at every MTO cycle. Due to this, the intermediate densities are reduced. However, if fewer iterations are permitted per MTO cycle, it is advisable to set and to values less than 1, in order to avoid excessive refinement and coarsening. The tuning of all these metaparameters forms an optimization problem in itself, and as adaptive design approaches become more sophisticated, setting such parameters can become highly nontrivial and timeconsuming. For the present study, no extensive parameter tuning was performed, yet already significant performance gains are observed. We see opportunities for future research in further adaptive and intelligent tuning strategies of metaparameters during the adaptive optimization itself, to take this burden away from the user.
For the MTO method, adaptivity serves as an addon where the design distribution and shape function orders are adapted at every cycle of refinement based on a predefined criterion. However, there are additional aspects of MTO which can be adapted to gain further improvements in accuracy and associated computational cost. Among others, appropriately adapting the filter radius could lead to further improvements. In the context of adaptive refinement, the impact of adaptive filter radius has been explored in [19]. For MTO, this aspect has briefly been discussed in [28]. However, the advantage of using an adaptive filter radius in an MTO setting remains an open question and needs to be explored. Also, based on our numerical tests, here we chose to use an adaptive stopping criterion, such that the cutoff value for minimum change in objective value between successive iterations is relaxed by a factor of 0.6 at every adaptive cycle. However, further investigations are needed to decide how this aspect can be adapted in the most efficient way.
Additional directions associated with adaptivity exist that could be investigated for further improvement of the methodology. For example, currently the number of design variables is set to the maximum allowed value based on the elementlevel upper bound described in [30]. However, it is still an open question whether this is the most appropriate way to refine the design field. Moreover, for the problems presented in this paper, we observed that for the chosen setting, violating the systemlevel bounds (also derived in [30]) did not have any detrimental impacts. Hence, we decided to not incorporate the systemlevel bounds in the method. However, for more complex problems, where the objective is very sensitive to small design changes, the systemlevel bounds might have to be enforced.
To wrap up the discussions, there are several research aspects that can be explored in the context of adaptive MTO. This work lays the foundation for an adaptive MTO scheme that is mathematical reliable as well as computationally efficient. It is hoped that with further research along the directions outlined above, the proposed approach can be improved further.
6 Conclusions
Multiresolution topology optimization (MTO) methods decouple the analysis and design discretizations, such that high resolution design representations are permitted on relatively coarse analysis meshes. In this paper, the first adaptive variant of the MTO scheme, namely adaptive MTO, has been presented. Through several 2D numerical examples, it has been demonstrated that the proposed method can obtain highly optimized designs at significantly lower computational cost than in conventional MTO, and high analysis accuracy. Moreover, undesired features such as intermediate densities and QRpatterns can be significantly reduced in the resulting designs, and a desired analysis accuracy can be enforced. A particularly interesting application of this adaptive MTO method is for TO problems involving low material volume fractions. The speedup over conventional MTO was found to increase with decreasing material volume fraction. It has been shown that for test cases with a 10% maximum relative volume, 10fold speedup can be obtained over the conventional MTO scheme in 2D, when the adaptive MTO method is used. For 3D problems, even higher speedups are expected.
Clearly, the proposed adaptive approach improves on the conventional MTO method by tacking some of the issues associated with it. For future work, we aim at exploring the application of adaptive approach for problems involving different physics and threedimensional problems. However, based on the results presented in this study, it can already be argued that the proposed approach could serve as an important methodology to obtain high resolution designs at an attractive computational cost.
Acknowledgements
This work is part of the Industrial Partnership Programme (IPP) Computational Sciences for Energy Research of the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organisation for Scientific Research (NWO). This research programme is cofinanced by Shell Global Solutions International B.V. The numerical examples presented in this paper are implemented using deal.II [46], an opensource C++ package for adaptive finite element analysis. We would like to thank the developers of this software. We would also like to thank Rajit Ranjan for his help.
Appendices
Appendix A means clustering
means clustering is a cluster analysis technique popularly used in data mining
[34]. It aims to partition observations into clusters such that the observations in each cluster tend to be close to each other. Note that although the problem is computationally difficult, there are various heuristic techniques that can quickly obtain a locally optimal solution.This technique can be used to choose locations of design points within a finite element (FE) and one of the primary advantages of this method is that it is easily applicable to various finite elements differing in geometry. Synonymous to the observations required in means clustering, a large number of uniformly distributed random points are chosen within the FE using Mersenne twister pseudorandom number generator [47]. Given that design points’ locations need be to chosen in the FE, we choose . Next, an initial set of points is chosen in the FE using means++ cluster center initialization algorithm [35]. These points serve as the initial means for the observations.
Let denote the initial locations of design points, then the following two steps are iteratively performed to optimize these locations:

Assignment step: Each observation is assigned to exactly one out of clusters based on the shortest Euclidean distance. Thus, during the iteration, is assigned to the cluster, if
(15) 
Update step: The new centroids of each of the clusters then become the new locations of the design points. The centroids are calculated as follows:
(16)
The two steps are repeated until locally optimal clustermeans are obtained. Note that for every number of design points, these distributions are generated once, and stored for use during optimization.
Fig. 21 shows the initial and optimized distributions of 40 design points in a Qtype FE. The optimized design distribution has been obtained using the means clustering algorithm. Clearly, in the optimized design field, the design points are more uniformly distributed and away from the boundaries of the element.
Appendix B Numerical integration scheme
The element stiffness matrix needs to be accurately integrated for every finite element. For the traditional TO using Q1 elements with elementwise constant densities, a Gauss quadrature rule is sufficient. However, for more complex density fields and higher order shape functions, more advanced ways of integration are needed to obtain correct . One of the possibilities is to use higher order integration schemes. A drawback of this approach is that a solidvoid boundary may not be correctly modeled. However, the associated error is very small, and with higher order integration schemes, numerically correct designs are obtained using MTO.
The density inside every voxel in the background mesh is constant. Thus, a composite integration scheme can also be used, where the voxelcontributions to the stiffness matrix are evaluated first, and these are then summed together to obtain the element stiffness matrix [31]. Since density is assumed to be constant inside each voxel, the choice of integration scheme depends on the polynomial order of the shape functions only. The advantage of this scheme is that the solidvoid boundaries are aligned with the edges of the voxels, due to which the stiffness matrix can be accurately integrated.
The composite integration, in general, is superior over the traditional integration scheme which is based on higher order Gauss quadrature rule. However, since in TO the design changes during the course of optimization, significant amount of information related to the stiffness matrices needs to be precomputed to use it in an adaptive MTO formulation. To avoid this excessive storage issue and to reduce the additional computational costs related to assembling the stiffness matrix at each iteration of MTO, we prefer to use the traditional Gauss quadrature rule with higher number of integration points.
Gauss quadrature rule  
1  0  1  4  2  
4  2  1  4  4  
9  3  2  9  7  
16  5  3  16  11  
25  6  3  16  12  
36  7  4  25  15  
49  9  5  36  19  
64  10  5  36  20  
Table 3 lists the minimum Gauss quadrature rule needed to accurately integrate the element stiffness matrix for several different density fields and polynomial shape functions. Here, only quadrilateral finite elements are considered. Based on the number of design points, a polynomial design field is constructed, and based on the shape functions, the order of element stiffness matrix is determined.
References
 [1] M. P. Bendsøe. Optimal shape design as a material distribution problem. Struct. O., 1(4):193–202, 1989.
 [2] O. Sigmund and K. Maute. Topology optimization approaches: A comparative review. Structural and Multidisciplinary Optimization, 48:1031–1055, 2013.
 [3] N. P. van Dijk, K. Maute, M. Langelaar, and F. Van Keulen. Levelset methods for structural topology optimization: a review. Struct. Multidiscip. O., 48:437–472, 2013.
 [4] J. D. Deaton and R. V. Grandhi. A survey of structural and multidisciplinary continuum topology optimization: post 2000). Structural and Multidisciplinary Optimization, 49:1–38, 2014.
 [5] N. Aage and B. S. Lazarov. Parallel framework for topology optimization using the methods of moving asymptotes. Struct. Multidiscip. O., 47(4):493–505, 2013.
 [6] S. Wang, E. Sturler, and G. H. Paulino. Largescale topology optimization using preconditioned Krylov subspace methods with recycling. Int. J. Numer. Meth. Eng., 12:2441–2468, 2007.
 [7] K. Suresh. Efficient generation of largescale paretooptimal topologies. Struct. Multidiscip. O., 47:49–61, 2013.
 [8] O. Amir, N. Aage, and B. S. Lazarov. On multigridCG for efficient topology optimization. Struct. Multidiscip. O., 49(5):815–829, 2014.
 [9] Z. Xia, Y. Wang, Q. Wang, and C. Mei. GPU parallel strategy for parameterized LSMbased topology optimization using isogeometric analysis. Struct. Multidiscip. O., 56:413–434, 2017.
 [10] O. Amir, M. P. Bendsøe, and O. Sigmund. Approximate reanalysis in topology optimization. International Journal for Numerical Methods in Engineering, 78(12):1474–1491, 2009.
 [11] O. Amir. Revisiting approximate reanalysis in topology optimization: on the advantages of recycled preconditioning in a minimum weight procedure. Struct. Multidiscip. O., 51(1):41–57, 2015.
 [12] J. MartinezFrutos and D. HerreroPerez. Efficient matrixfree GPU implementation of Fixed Grid Finite Element Analysis. Finite Elem. Anal. Des., 104:61–71, 2015.
 [13] I. Babuška and B. Q. Guo. The h, p and hp version of the finite element method: basis theory and applications. Adv. Eng. Softw., 15:159–174, 1992.
 [14] Matteo Bruggi and Marco Verani. A fully adaptive topology optimization algorithm with goaloriented error control. Comput. and Struct., 89(1516):1481–1493, 2011.
 [15] K. Maute and E. Ramm. Adaptive topology optimization. Struct. O., 10(2):100–112, 1995.
 [16] F. van Keulen and E. Hinton. Topology design of plate and shell structures using the hard kill method. In B.H.V. Topping, editor, Advances in structural Engineering Optimization, pages 177–188. Civil Comp. Press, Edinburgh, 1996. presented at the Third International Conference in Computational Structural Technology in Budapest, 2123 August, 1996.
 [17] Joao Carlos Costa and Marcelo Krajnc Alves. Layout optimization with hadaptivity of structures. International Journal for Numerical Methods in Engineering, 58(1):83–102, 2003.
 [18] R. Stainko. An adaptive multilevel approach to the minimal compliance problem in topology optimization. Communications in Numerical Methods in Engineering, 22(2):109–118, 2006.
 [19] D. K. Gupta, M. Langelaar, and F. van Keulen. Combined mesh and penalization adaptivity based topology optimization. In Proceedings, 57th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, AIAA, pages 1–12, 2016.
 [20] H. NguyenXuan. A polytreebased adaptive polynomial finite element method for topology optimization. Int. J. Numer. Meth. Engg., 110:972–1000, 2017.
 [21] T. Y. S. Hoshina, I. F. M. Menezes, and A. Periera. A simple adaptive mesh refinement scheme for topology optimization using polygonal meshes. J. Braz. Soc. Mech. Sci. & Eng., 40(348):1–17, 2018.
 [22] A. B. Lambe and A. Czekanski. Topology optimization using a continuous density field and adaptive mesh refinement. Int. J. Numer. Meth. Eng., 113(3):357–373, 2018.
 [23] J. Wu. Continuous optimization of adaptive quadtree structures. Comput. Aided Des., 102:72–82, 2018.
 [24] J. Parvizian, A. Duster, and E. Rank. Finite cell method: h and pextension for embedded domain problems in solid mechanics. Comput. Mech., 41:121–133, 2007.
 [25] J. Parvizian, A. Duster, and E. Rank. Topology optimization using the finite cell method. Optim. Eng., 13, 2012.
 [26] Tam H. Nguyen, Glaucio H. Paulino, Junho Song, and Chau H. Le. A computational paradigm for multiresolution topology optimization (MTOP). Structural and Multidisciplinary Optimization, 41(4):525–539, 2010.
 [27] T. H. Nguyen, C. H. Le, and J. F. Hajjar. Topology optimization using the pversion of the finite element method. Structural and Multidisciplinary Optimization, 56(3):571–586, 2017.
 [28] Tam H Nguyen, Glaucio H Paulino, Junho Song, and Chau H Le. Improving multiresolution topology optimization via multiple discretizations. International Journal for Numerical Methods in Engineering, 92:507–530, 2012.
 [29] Y. Wang, J. He, and Z. Kang. An adaptive refinement approach for topology optimization based on separated density field description. Computers and Structures, 117:10–22, 2013.
 [30] D. K. Gupta, G. .J. van der Veen, A. Aragón, M. Langelaar, and F. van Keulen. Bounds for decoupled design and analysis discretizations in topology optimization. Int. J. Numer. Meth. .Eng., 111(1):88–100, 2017.
 [31] J. P. Groen, M. Langelaar, O. Sigmund, and M. Reuss. Higherorder multiresolution topology optimization using the finite cell method. Int. J. Numer. Meth. .Eng., 110(10):903–920, 2017.
 [32] D. K. Gupta, M. Langelaar, and F. van Keulen. QRpatterns: Artefacts in Multiresolution Topology Optimization. Struct. Multidiscip. O., 58(4):1335–1350, 2018.
 [33] E. Rank. Adaptive remeshing and hp domain decomposition. Comput. Method. Appl. M., 101:299–313, 1992.
 [34] D. Mackay. Information Theory, Inference and Learning Algorithms. Cambridge University Press, 2003.
 [35] D. Arthur and S. Vassilvitskii. Kmeans++: The Advantages of Careful Seedings. In SODA ‘07: Procc 18th Annual ACMSIAM Symp, Discrete Algorithms, pages 1027–1035, 2007.
 [36] T. E. Bruns and D. A. Tortorelli. Topology optimization of nonlinear elastic structures and compliant mechanisms. Comput. Method. Appl. M., 190:3443–3459, 2001.
 [37] D. W. Kelly, J. R. Gago, O. C. Zienkiewicz, and I. Babuska. A posteriori error analysis and adaptive proceses in the finite element method. Part I  Error analysis. Int. J. Numer. Meth. Eng., 19:1593–1619, 1983.
 [38] M. Ainsworth and J. T. Oden. A Posteriori Error Estimation in Finite Element Analysis. John Wiley & Sons, Inc., USA, 2000.
 [39] T. Grätsch and K. Bathe. A posteriori error estimation techniques in practical finite element analysis. Comput. and Struct., 83:235–265, 2005.
 [40] O. C. Zienkiewicz, D. W. Kelly, J. Gago, and I. Babuska. Hierarchical finite element approaches error estimates and adaptive refinement.
 [41] J. van Schoubroeck, D. K. Gupta, and M. Langelaar. An investigation of the effect of initial design choices in densitybased topology optimization (submitted). Struct. Multidiscip. O., 2018.
 [42] S. RojasLabanda and M. Stolpe. Benchmarking optimization solvers for structural topology optimization. Struct. Multidiscip. O., 52(3):527–547, 2015.
 [43] D. K. Gupta, M. Langelaar, M. Barink, and F. van Keulen. Topology optimization of front metallization patterns for solar cells. Struct. Multidiscip. O., 51:941–955, 2015.
 [44] M. P. Bendsøe and O. Sigmund. Topology optimization: Theory, methods and applications. Springer, Berlin, 2003.
 [45] J. T. Oden and S. Prudhomme. Goaloriented error estimation and adaptivity for the finite element method. Comput Math Appl, 41:735–756, 2001.
 [46] W. Bangerth, R. Hartmann, and G. Kanschat. deal.II – a general purpose object oriented finite element library. ACM Transactions on Mathematical Software, 33(4):24/1–24/27, 2007.
 [47] M. Matsumoto and T. Nishimura. Mersenne twister: A 623dimensionally equidistributed uniform pseudorandom number generator. ACM Trans. Model. Comput. Simul., 8(1):3–30, 1998.