Design and analysis adaptivity in multi-resolution topology optimization

Multiresolution topology optimization (MTO) methods involve decoupling of the design and analysis discretizations, such that a high-resolution design can be obtained at relatively low analysis costs. Recent studies have shown that the MTO method can be approximately 3 and 30 times faster than the traditional topology optimization method for 2D and 3D problems, respectively. To further exploit the potential of decoupling analysis and design, we propose a dp-adaptive MTO method, which involves locally increasing/decreasing the shape function orders (p) and design resolution (d). The adaptive refinement/coarsening is performed using a composite refinement indicator which includes criteria based on analysis error, presence of intermediate densities as well as the occurrence of design artefacts referred to as QR-patterns. While standard MTO must rely on filtering to suppress QR-patterns, the proposed adaptive method ensures efficiently that these artefacts are suppressed in the final design, without sacrificing the design resolution. The applicability of the dp-adaptive MTO method is demonstrated on several 2D mechanical design problems. For all the cases, significant speed-ups in computational time are obtained. In particular for design problems involving low material volume fractions, speed-ups of up to a factor of 10 can be obtained over the conventional MTO method.


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 density-based 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 three-dimensional (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 large-scale 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 h-refinement and p-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 Zienkiewicz-Zhu error estimator and (b) the shortest distance to the material-void 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 material-void boundary. Bruggi and Verani [14] progressed in the direction of the work proposed by [16], and proposed a goal-based AMR strategy that properly guides the progression of refinement and coarsening in TO. For refinement/coarsening, a dual-weighted residual based FEA indicator as well as a heuristic density-gradient based indicator were used. While most of these methods helped to achieve the required h-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 h-refinement and very clear material descriptions with low gray regions were obtained. Other adaptive formulations involving h-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 p-refinement, as stated earlier, where the mesh topology remains the same. Additionally, for smooth problems, the accuracy of p-refinement is dramatically higher than that of h-refinement for the same computational costs [13]. Increasing the polynomial order of the shape functions gives an exponential rate of convergence. Other advantages of p-refinement are its robustness against locking effects and high aspect ratios [24]. However, due to the fact that the conventional TO approaches assume an elementwise-constant density distribution, using higher-order 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 FE-based 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 material-void 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 elementwise-constant 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 p-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 'multi-resolution' 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 well-posedness 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, FCM-based 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 FCM-based TO. Their observations show close resemblance with those in [30]. Also, the authors showed that using FCM-based TO, remarkable speed-ups of more than 3-and 60-folds for 2D and 3D problems, respectively, could be obtained over the traditional TO approach. However, for certain configurations of FCM-based TO, it is possible that the design consists of 'QR-patterns', 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 QR-pattern artifacts [27,30,31], but has the undesired consequence of reducing the design resolution. Applying p-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 FCM-based 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 QR-patterns, 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 QR-pattern artefacts. We coin the term 'dp-adaptivity', an adaptive multiresolution TO scheme where both the design resolution d and FE polynomial order p can be locally adapted based on certain refinement/coarsening criteria. Here, the symbol 'd' should not be confused with the one in hp-d adaptivity, where it refers to domain decomposition and mesh overlaying [33]. It is assumed that computational costs are the limiting factor, and that the manufacturing-imposed 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 non-adaptive 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 QR-patterns 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 dp-adaptivity approach. The applicability of this

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 non-uniform. 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 non-uniformly in the domain. The overlapping background element comprises 3 × 3 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 k-means clustering method is used [34,35]. This approach divides the design domain into k segments (clusters) with roughly equal areas. The design points are assumed to be located at the centroids of these clusters. For self-containment, 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 equally-sized 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, 3 × 3 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 P 1 (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 element-level implementation. While choosing the local projection radius P 1 , 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 2 × 2 MTO elements, each comprising a Q1 finite element and 3 × 3 density cells. As can be seen, the distribution of design points can be non-uniform. The four MTO elements from top-left to bottom-right consist of 4, 9, 3, and 7 design points, respectively. In the bottom-right 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 P 2 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 P 2 (Fig. 1). Fig.  3 illustrates how these density values are computed. It shows a mesh of 2 × 2 MTO elements, comprising Q1 finite elements and the corresponding background domain with 3 × 3 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 distance-weighted projection of design densities found in the neighborhood of a certain radius R over the background mesh. In this work, density filtering is used for the projection [36].
The use of the background mesh facilitates d-adaptivity, i.e. the use of different numbers of design points in adjacent elements. In the absence of the background mesh, the non-uniform 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.

Mathematical formulation
In this paper, the applicability of a dp-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 where, J (·) denotes the objective functional, and K, u and f denote the global stiffness matrix, displacement vector and load vector, respectively. The vector z 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 V 0 .
Next, the details related to various steps associated with the proposed multiresolution modeling approach are described. The matrix K in Eq. 1 is obtained from the global assembly of the element stiffness matrices K e , which can be expressed as where B and D denote the strain-displacement matrix and constitutive matrix, respectively, and N g is the number of integration points. More details related to the choice of numerical integration are discussed in Appendix B. The subscript i refers to the i th integration point and w i denotes the respective integration weight. The construction of the D 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 where E 0 is the Young's modulus of the solid material and E min is a very small value (typically 10 −9 E 0 ) used to avoid singularity of the system stiffness matrix. Also,ρ i denotes the density calculated at the i th integration point, q is the penalization power and D 0 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 P 2 based on the density filtering method which is widely used in TO [36]. Mathematically, it can be stated as whereρ refers to density values for the cells contained in the background mesh with their centers lying within a distance R from the corresponding integration point (Fig. 3), and their number is denoted by nρ. Here, terms H ij reduce linearly with distance from the integration point, i.e., where dist(·) denotes the Euclidean distance operator. As stated in Section 2.1, the background mesh densities are calculated using the P 2 projection from the design mesh to the background mesh. For the p th MTO element, the density of the q th density cell is given asρ where, ρ s refers to the density value associated with the s th design point in the design domain contained within the p th MTO element, and lying within a distance r p from the centroid of its q th density cell. The number of such design points is denoted by n ρ , and r p is the radius of the projection for the p th element (Fig. 2). Here, h qs is defined as h qs = r p − dist(q, s).
As stated earlier, the projection radius r p 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 where dim denotes problem dimension, and L p is the edge-length of the p th MTO element. The operator · denotes ceiling function which rounds the contained floating-point number to the nearest greater integer value. The term Lp d 1/dim refers to edge-length of the density cells. Next, to obtain a projection length slightly larger than the diagonal, we multiply by 1.04(dim) 0.5 . Note that Eq. 8 has been obtained empirically through observations based on various design distributions obtained using the k-means clustering approach. For other approaches of choosing the locations of design points, where for any value of d, the distance between the design points can be provided mathematically, it is possible that even lower values of r p work. Lower values of r p can help to further reduce the loss in design resolution caused due to the choice of localized projection P 1 , and this could be a potential direction for future research. Fig. 4a 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 8 × 4 finite elements. For each MTO element, 225 design points, distributed in a square grid of 15 × 15, 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 element-level 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 p equal to 10, the number of deformation modes is 239, which is greater than 225. With p and d equal to 10 and 225, respectively, the MTO elements are referred as Q10/d225 type elements. For this example, the projection radius R is set to 0.13 times the element-length, which is equivalent to the size of 2 density cells. Fig. 4b 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. 4b, there are parts of the domain where even lower-order 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.

General description of the method
We present here a dp-adaptive version of the MTO method which is capable of enhancing further the ratio between the design resolution and analysis cost compared to non-adaptive MTO. The proposed MTO method efficiently distributes the design variables and locally adapts (increases/decreases) the polynomial order of the shape functions. A three-part 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 p and d) 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 dp-adaptive MTO, a cycle of TO is performed, using a certain initial uniform designand FE-discretization. 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 error-based indicator, a density-based indicator, and a QR-based indicator.
Here, QR-error refers to the error due to the incapability of the chosen shape function in modeling the displacement field arising from a high-resolution 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 d and p values for the whole domain, constitute one cycle of dp-adaptivity. The general structure of a dp-adaptive MTO cycle is as follows: 1. Perform optimization of an MTO problem with fixed p and d values. With the new dp-adapted mesh, the next cycle of TO is performed. Section 3.3 below describes each of the above steps in detail.

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 dp-adaptive scheme in itself is not dependent on the choice of refinement indicator, and can be coupled with other appropriate indicators as well.

Analysis-based 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 u 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 (p). 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., goal-oriented 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  Figure 6: Bounds for the design refinement indicator as a function of the adaptive cycle [19].
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 QR-error 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 p-adaptivity based on the Kelly error estimator, the analysis residual error vector Γ a needs to be defined. For the i th MTO element, Γ a i can be computed as: where, F refers to a face (edge in 2D) of the element and operator [·] denotes the jump in the argument across face F . Also, ∂i denotes the set of all faces of the element. The constant term c F is set to h F 2p F , where h F is the element diagonal and p F denotes the maximum among the polynomial degrees of the adjacent elements [40]. The residual errors Γ a are ranked, and the top 10% and bottom 5% of the elements are selected for increasing and decreasing the p values, respectively.
For illustration purposes, we perform a partial adaptive MTO run on the problem shown in Fig. 4a. Fig. 5a shows the optimized cantilever beam design obtained for this problem after one TO cycle. The design has been obtained on a mesh of 40 × 20 Q2 finite elements with 4 × 4 design points per element. The optimized design clearly shows typical artefacts (QR-patterns) of disconnected structural features. Fig. 5b shows the distribution of polynomial shape function orders obtained from p-adaptivity controlled by only the analysis-based refinement indicator. It is observed that coarsening (reduction in p) has mainly occurred in the void cells which are far from material-void boundaries. This is because the jumps in displacement gradients across the edges for these elements are zero. For refinement (increase in p), the elements at the boundary have been preferred.

Density-based refinement indicator
The density-based 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 high-resolution 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)  The optimized design used for adaptive refinement is shown is shown in Fig. 5. 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 (r l ≤ ρ ≤ r u ) and coarsening (ρ < c l or ρ > c u ) 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 density-based refinement indicator Γ d k for the k th MTO element is set to 0. To update Γ d k , we iterate over all the density cells of the k th MTO element and consider the sum of individual refinement or coarsening contributions of these cells. Let n d,k denote the number of density cells contained within the background mesh associated with the k th MTO element. Then Γ d k is updated as follows: • Iterate over j from 1 to n d,k : 1. Let the density of the j th voxel be denoted by ρ kj .
Here, the average density ρ avg is defined using the expression ρ avg = (ρ max + ρ min )/2. The variables r l , r u , c l and c u are the bounds used to characterize the refinement and coarsening zones as shown in Fig. 6, and are defined as follows: Here,k 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 dp-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. 4a. The shape function field (Fig. 7a) denotes the polynomial order of the shape functions used in every finite element. The design field (Fig. 7b) 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 density-based refinement indicator. From Fig. 7, it is seen that the material-void boundaries where the intermediate densities are prominent, have primarily been refined. Coarsening occurs in void parts of the domain.

QR-error 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 QR-error. A closed-form condition to predict this QR-error 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 dp-adaptivity, QR-errors must be quantified at element level. We have proposed a method in [32], where an approximation to the QR-error can be obtained for any element through a comparison with a reference solution obtained by local p-refinement. In this work, we use this cost-effective local QR-error indicator proposed in [32]. Once a sufficiently converged design has been obtained from a TO cycle, the QR-error is determined by evaluating the effect of local p-refinement, as follows.
Let K k denote the element stiffness matrix, displacement solution and internal load vector for the k th MTO element. Here, p denotes the polynomial degree of the shape functions used in this element. Let u (p+1) k denote the displacement solution obtained for the k th element using shape functions of polynomial order p + 1. Note that u (p+1) k will be obtained by solving the element-level system K Here, nodal load f where J p k refers to element-level strain energy for the k th finite element using shape functions of order p.
k have been used. This strain-energy-based criterion (Eq. 14) has been found to work well for the cases shown in this paper. Fig. 8a and 8b show an optimized design obtained after first cycle of MTO run, and the corresponding error distribution obtained using the QR-error indicator for the problem shown in Fig. 4a. Since the elementlevel test for QR-error is very conservative, it predicts higher error values compared to the actual full-scale TO problem [32]. Thus, to avoid undesired excessive increase in the values of p, we restrict the increment of p by only 1 per adaptive cycle based on the QR-error 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. 8c. It is observed that the regions where the QR-patterns exist, have been flagged for refinement. Moreover, elements at the material boundaries, which are partially void or solid, also show high value of QR-error and are flagged. An interesting observation in Fig. 8b is that the elements which are completely void or solid also show QR-error values in the range 0.3-0.5. Although significant, the QR-error 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 QR-error values in these regions is the use of low order shape functions. For low values of p, the displacement solution for even a uniform density field may not be accurately modeled. When solving element-level FE problems with low shape function orders p and p + 1, it is observed that the modeling accuracy significantly improves when p is increased. Due to this, nonzero large values of Q k are recorded in solid and void parts as well.

dp-adaptivity algorithm
The different steps of dp-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.
1. Once a cycle of TO run is completed, get the optimized design for dp-adaptivity.

Perform p-adaptivity based on analysis error criterion.
(a) Update Γ a = {Γ a 1 , Γ a 2 , . . . , Γ a n el } values for the whole analysis mesh (discussed in Section 3.2.1), where Γ a i is the analysis error indicator value for the i th MTO element. (b) Sort Γ a in ascending order such that a corresponding ordered setΓ a is obtained.
(c) Set the refine/coarsen flag of the k th element Θ k to -1 for the first α d c fraction of the MTO elements inΓ a , and Θ k = 1, for the last α a r fraction of the elements. Here, −1 and 1 denote that the cell has been flagged for coarsening (decrease in p value) and refinement (increase in p value), respectively. For no refinement/coarsening, Θ k is set to 0. (c) Update p-values by iterating over k from 1 to n el : i. For the first α d c fraction of the elements inΓ d , do: A. if p k = 1, set Θ k = −2. This helps to identify that the current element has been checked for coarsening. Since p k cannot be lower than 1, no coarsening is performed. B. if p k > 1 and Θ k = 0, set p k = p k − 1.
ii. For the last α d r fraction of the elements inΓ d , do: A. if Θ k = 0 or Θ k = −1, set p k = p k + 1. This means that if the element has been coarsened or left untreated based on the analysis indicator above, then refine it. ii. Find first perfect square (cube in 3D) number (d) greater than max(d el ).
iii. Set the number of density cells per MTO element equal tod. iv. Update projection radius r for every MTO element (Eq. 7).

Update p values to reduce the QR-error in every MTO element.
(a) Iterate over k from 1 to n el , do: i. Calculate the QR-error for the k th cell (discussed in Section 3.2.3).
ii. Update p k = p k + 1 for the k th element, if QR-error is greater than a certain error tolerance α QR .
The dp-adaptive MTO cycle is complete once the domain has been adaptively refined based on the three indicators. With the new dp-refined mesh, the next cycle of TO is performed.

Numerical tests 4.1 Definition of test problems
To demonstrate the applicability and effectiveness of dp-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 E 0 is set to 1 Nm −2 , ν = 0.3, and the SIMP penalization factor q is set to 3. The domain in each case is discretized using an initial mesh of 40 × 20 MTO elements, comprising quadrilateral finite elements with shape functions of polynomial order 2 and 4 × 4 design points per element. The radius R is set to 0.3h, where h is the edge-length of any MTO element in the mesh. As a stopping criterion for all the test cases used in this paper, the optimization process for thek th cycle is terminated when the change in objective value between two consecutive iterations is less than ∆J 1 × γ (k−1) . Here, ∆J 1 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, ∆J 1 and γ are set to 0.04 and 0.6, respectively, and these values have been found to work well. Based on this, the first (k = 1) and second (k = 2) 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 high-resolution traditional TO mesh with elementwise constant densities. In this paper, the reference mesh comprises 320 × 160 finite elements and the polynomial order p 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 F (Fig. 4a). For this case, F = 1 N and L = 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. 9a), and it is ensured that the load is consistently distributed over the various cycles of adaptivity. Here, F = 0.5N L and L = 1 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. 1 with z = f . The third case is a compliant mechanism problem where a force inverter needs to be designed, such that for a point load f in at one end, the displacement u out at the other end is maximized (Fig. 9b). Here, spring stiffnesses k in and k out are set to 1 Nm −1 and 0.001 Nm −1 , respectively. For the force inverter, z in Eq. 1 is a vector of zeros with 1 contained at the DOF where u out needs to be maximized. Thus, z = [0 . . . 0 1 0 . . . 0] . The flexure hinges that are formed in this compliant mechanism problem will have sub-element resolution, and this aspect makes also this problem an interesting test for our method.

Results
Here, we discuss the results obtained for the three test problems using a dp-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 dp-adaptive MTO schemes for 2D problems, and for this purpose, we will compare its performance with the non-adaptive MTO scheme, implemented in the same framework and evaluated on the same computing hardware. Figure 10: Optimized cantilever designs for the point load case shown in Fig. 4a, obtained using (a) a uniform MTO mesh and (b) dp-adaptive MTO approach. The maximum permissible material volume fraction is set to 0.45. A 4.5-fold speed-up as well as a superior objective value are obtained using dp-adaptivity. Additional information related to this test case is listed in Table 1. Table 1: Numerical findings of several dp-adaptive MTO cases. For all the cases, the domain has been discretized using 40 × 20 MTO elements, and the initial polynomial order of the shape funtions is set to 2 for every element. Each MTO element initially consists of 16 design points and the projection radius R is set to 0.3h, where h denotes element size. The maximum permissible values for shape function order p max and number of designs points d max are set to 5 and 64, respectively. For the reference solution, a globally uniform mesh comprising 320 × 160 finite elements with p = 3 is used. Below, V 0 denotes maximum allowed volume fraction of material, J and J 0 are the objective values for dp-adaptive MTO run and the non-adaptive MTO run, and J * denotes the reference solution. The N d and DOFs denote number of design points and free degrees of freedom employed in the last cycle of dp-adaptive MTO run.

Problem
Definition   Fig. 4a. The first design (Fig. 10a) has been obtained using the traditional non-adaptive MTO scheme, and the other (Fig. 10b) by our dp-adaptive approach. For the two cases, the maximum allowed material volume fraction V 0 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 dp-adaptive MTO result is the reduced computational cost. Adding the dp-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 N d and free DOFs used in the dp-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 th ) cycle of the dp-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 speed-up is the reduced number of iterations required in the final cycle of the dp-adaptive method under the same stopping criterion as used for the non-adaptive 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 non-adaptive MTO scheme uses a total of 56 iterations. Table 2 provides the details related to the dp-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.
In terms of performance, the cantilever design obtained from the dp-adaptive approach slightly outperforms the design obtained using non-adaptive MTO. The obtained performance ratio J /J 0 is equal to 0.98, where J and J 0 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 J /J * = 0.98, where J and J * 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 dp-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 dp-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 QR-patterns arising from the limitations of low order polynomial shape functions in those parts of the design [32]. Over the next cycles, p-refinement occurs in those regions and the QR-patterns are eliminated. Since the design points are distributed in the domain using k-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, (d) Cycle: 4 Figure 11: Optimized designs (right), and the respective shape function orders (middle) and design field (left) obtained for 4 cycles of dp-adaptive MTO run for a cantilever beam subjected to point load (Fig. 4a). The initial mesh is uniform and each element has shape functions of polynomial order 2 and 16 design points per element. The maximum allowed shape function order and number of design points are restricted to 5 and 64 per element, respectively. dp-adaptivity could be potentially beneficial. To investigate this, we study two additional cases of the point load cantilever beam involving lower values of V 0 . Fig. 12 shows the optimized designs for V 0 = 0.20 using conventional MTO (Fig. 12a) and dp-adaptive method (Fig. 12b), respectively. For V 0 = 0.20, the computational time advantage has increased to a factor of 8.3. Also, it is seen that the design obtained using the non-adaptive MTO method differs significantly from the result of dp-adaptivity. Moreover, in terms of performance, the design obtained using dp-adaptivity is relatively less compliant. The ratio J /J 0 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 V 0 = 0.10. Fig. 13 shows the optimized designs for this volume fraction obtained using the conventional MTO method and dp-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 speed-up of 10 times is observed. The increase in speed-up 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 non-adaptive MTO. For this case, it is observed that J /J 0 is 1.03, which implies that the design obtained using dp-adaptivity is slightly inferior to that obtained using the non-adaptive version. The analysis accuracy is also slightly lower than in the previous cases, with J /J * = 0.96.
An understanding on the convergence of the dp-adaptive MTO process for V 0 = 0.10 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 V 0 = 0.45, it is observed that QR-patterns are formed here as well, which are removed Figure 12: Optimized cantilever designs for the point load case shown in Fig. 4a, obtained using a uniform MTO mesh (left) and dp-adaptive MTO approach (right). The maximum permissible material volume fraction is set to 0.20. A speed-up of 8.3 times is obtained using dp-adaptivity. Additional information related to this test case is listed in Table 1.  Fig. 4a, obtained using a uniform MTO mesh (left) and dp-adaptive MTO approach (right). The maximum permissible material volume fraction is set to 0.10. A 10-fold speed-up is obtained using dp-adaptivity. Additional information related to this test case is listed in Table 1.
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 non-adaptive as well as the dp-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 dp-adaptive method for V 0 = 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.

Compliance minimization for distributed load
For the cantilever beam subjected to a distributed load (Fig. 9a), V 0 is set to 0.45. Fig. 15 shows the optimized designs obtained using a uniform MTO mesh (Fig. 15a) and the dp-adaptive approach (Fig. 15b). 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 non-adaptive and adaptive variants of MTO are very similar. In terms of performance, a speed-up of 4.6 times is observed, and the accuracy of the obtained solution is close to 1. The obtained J /J 0 value is 0.98, which implies that the dp-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 dp-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 Optimized designs (right) and the respective shape function orders (middle) and design field (left) obtained for 4 cycles of dp-adaptive MTO run for a cantilever beam subjected to point load (Fig. 4a). The initial mesh is uniform and each element has shape functions of polynomial order 2 and 16 design points per element. The maximum allowed shape function order and number of design points are restricted to 5 and 64 per element, respectively. 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 non-design 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. 16a), the initial value of q is 3 and it is increased by 1 at every cycle. For the second case (Fig. 16b), 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 dp-adaptive refinement, see Fig.  17. Due to the low order of the shape function used in Cycle 1, QR-patterns 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 dp-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 dp-adaptive MTO are lower than in the conventional MTO method. Together, these two factors make the dp-adaptive MTO method 4.6 times faster in this case. Figure 15: Optimized cantilever designs for the distributed-load case shown in Fig.9a, obtained using a uniform MTO mesh (left) and dp-adaptive MTO approach (right). A 4.6-fold speed-up is obtained using dp-adaptivity.

Force inverter compliant mechanism
(a) dp-adaptive MTO (q =3, 4, 5 and 6) (b) dp-adaptive MTO (q =3, 5, 7 and 9) Figure 16: Optimized cantilever designs for the distributed-load case shown in Fig.9a, obtained using dpadaptive MTO approach. For both the cases, adaptive penalization has been used. For the 4 cycles of the dp-adaptive MTO run, the values of q used have been reported in the sub-captions.
To demonstrate the applicability of dp-adaptivity on topology optimization of compliant mechanisms, it is applied to the force inverter problem shown in Fig. 9b. The allowed volume fraction V 0 is set to 0.45 and the goal of the problem is to distribute the material in a way that the displacement u out is maximized. Fig. 18 shows the optimized designs obtained using conventional MTO (Fig. 18a) and the dp-adaptive method (Fig. 18b). 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 J /J 0 is 1.01. Since this is a maximization problem, a value of J /J 0 higher than 1 denotes that the design obtained using dp-adaptive MTO performs better. J /J * 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 dp-adaptivity. Similar to the other cases discussed in this paper, QR-patterns are observed in the results of the first cycle. Nevertheless, the overall material distribution after Cycle 1 already corresponds to the final solution. The QR-patterns 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.

Discussions
The primary goal of using an MTO scheme is to obtain a high-resolution 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 high-resolution 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, QR-patterns) where scope of improvement existed. The dp-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 dp-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 , and the respective shape function orders (middle) and design field (left) obtained for 4 cycles of dp-adaptive MTO run for a cantilever beam subjected to distributed load (Fig.  9a). The domain has been discretized using 40 × 20 quadrilateral finite elements (r = 0.3h). The initial mesh is uniform and each element comprises shape functions of polynomial order 2 and 16 design points per element. The maximum allowed shape function order and number of design points are restricted to 5 and 64 per element, respectively.
(a) MTO (J0 = 2.224m) (b) dp-adaptive MTO (J = 2.258m) Figure 18: Optimized cantilever designs for the force inverter problem shown in Fig.9b, obtained using a uniform MTO mesh (left) and dp-adaptive MTO approach (right). A speed-up of 6.2 folds is obtained using dp-adaptivity. : Optimized designs (right), and the respective shape function orders (middle) and design field (left) obtained for 4 cycles of a dp-adaptive MTO run for the force inverter problem shown in Fig. 9b. The initial mesh is uniform and each element comprises shape functions of polynomial order 2 and 16 design points per element. The maximum allowed order of the shape functions and number of design points are restricted to 5 and 64 per element, respectively.
involving other physics as well as in 3D settings. In [31], it has been shown that MTO can bring a speed-up of up to 32 folds over the traditional TO scheme. The improvement in 3D is significantly higher than that observed in 2D. As dp-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 dp-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 dp-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. 20a shows an optimized design obtained for a linear thermal conduction problem using the traditional TO approach. A mesh of 400 × 400 elements was used and R 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 dp-adaptive MTO method over its non-adaptive variant. Fig. 20b shows an optimized solar cell front metallization design obtained using the traditional TO approach on a mesh of 400 × 400 finite elements and R set to 1.5 times the element edge length [43]. This design has been obtained by solving a nonlinear electrical conduction problem, and only 4-5% of the domain is filled with material. For this case, it is seen that significant parts of the domain consists of void regions, (a) Optimized design for a linear heat conduction problem (b) Optimized metallization design for the front surface of a solar cell [43] Figure 20: Optimized designs obtained using the traditional TO approach on a mesh of 400 × 400 finite elements, with R set to 1.5 elements. The two cases refer to (a) linear heat conduction problem with V 0 set to 0.3, and (b) nonlinear electrical conduction problem [43].
which can be easily modeled with low values of d and p. Clearly, for such cases, the dp-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 dp-adaptivity could certainly have a potential value for problems where designs feature extended void regions.
To demonstrate the concept of dp-adaptivity, a composite indicator has been formulated in this paper. This indicator consists of an analysis error indicator, a density-based indicator and a QR-indicator. 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., goal-oriented 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 Γ a are chosen for refinement and coarsening, respectively. There is no particular motivation to choose these cut-offs. 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 density-based indicator, both α d r and α d c are set to 1.0 for the current study. This ensures that all the elements with Γ d > 0 are refined and all elements with Γ d < 0 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 α d r and α d c to values less than 1, in order to avoid excessive refinement and coarsening. The tuning of all these meta-parameters forms an optimization problem in itself, and as adaptive design approaches become more sophisticated, setting such parameters can become highly nontrivial and time-consuming. 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 meta-parameters during the adaptive optimization itself, to take this burden away from the user.
For the MTO method, dp-adaptivity serves as an add-on 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 R could lead to further improvements. In the context of adaptive h-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 dp-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 element-level 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 system-level bounds (also derived in [30]) did not have any detrimental impacts. Hence, we decided to not incorporate the system-level bounds in the method. However, for more complex problems, where the objective is very sensitive to small design changes, the system-level 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. Appendices A k-means clustering k-means clustering is a cluster analysis technique popularly used in data mining [34]. It aims to partition ψ observations into k 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 k-means clustering, a large number of uniformly distributed random points ψ are chosen within the FE using Mersenne twister pseudorandom number generator [47]. Given that k design points' locations need be to chosen in the FE, we choose ψ = 1000k. Next, an initial set of k points is chosen in the FE using k-means++ cluster center initialization algorithm [35]. These points serve as the initial k means for the ψ observations. Let m 1 , m 1 , . . . , m (k) 1 denote the initial locations of k design points, then the following two steps are iteratively performed to optimize these locations: 1. Assignment step: Each observation x p is assigned to exactly one out of k clusters based on the shortest Euclidean distance. Thus, during the t th iteration, x p is assigned to the i th cluster, if 2. Update step: The new centroids of each of the k clusters then become the new locations of the design points. The centroids are calculated as follows: The two steps are repeated until locally optimal cluster-means 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 Q-type FE. The optimized design distribution has been obtained using the k-means clustering algorithm. Clearly, in the optimized design field, the design points are more uniformly distributed and away from the boundaries of the element.

B Numerical integration scheme
The element stiffness matrix K e needs to be accurately integrated for every finite element. For the traditional TO using Q1 elements with elementwise constant densities, a 2 × 2 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 K e . One of the possibilities is to use higher order integration schemes. A drawback of this approach is that a solid-void 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 voxel-contributions 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 solid-void 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. Table 3: Choice of integration scheme for different combinations of design fields and polynomial shape functions for Q-type finite elements. Here, d and p denote the number of design points and polynomial order of the shape functions, respectively, n sup refers to the number of support points, and P(d) and P(K) denote the maximum possible polynomial order of the design field and stiffness matrix, respectively.  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.