Topology optimization with supershapes

This work presents a method for the continuum-based topology optimization of structures whereby the structure is represented by the union of supershapes. Supershapes are an extension of superellipses that can exhibit variable symmetry as well as asymmetry and that can describe through a single equation, the so-called superformula, a wide variety of shapes, including geometric primitives. As demonstrated by the author and his collaborators and by others in previous work, the availability of a feature-based description of the geometry opens the possibility to impose geometric constraints that are otherwise difficult to impose in density-based or level set-based approaches. Moreover, such description lends itself to direct translation to computer aided design systems. This work is an extension of the author’s group previous work, where it was desired for the discrete geometric elements that describe the structure to have a fixed shape (but variable dimensions) in order to design structures made of stock material, such as bars and plates. The use of supershapes provides a more general geometry description that, using a single formulation, can render a structure made exclusively of the union of geometric primitives. It is also desirable to retain hallmark advantages of existing methods, namely the ability to employ a fixed grid for the analysis to circumvent re-meshing and the availability of sensitivities to use robust and efficient gradient-based optimization methods. The conduit between the geometric representation of the supershapes and the fixed analysis discretization is, as in previous work, a differentiable geometry projection that maps the supershapes parameters onto a density field. The proposed approach is demonstrated on classical problems of 2-dimensional compliance-based topology optimization.


Introduction
Density-based and level set topology optimization methods have been employed with wide success to produce organic, free-form designs across multiple application realms. Despite this success, a withstanding challenge in topology optimization is the ability to impose geometric considerations to render designs that are amenable to production. As most industry practitioners of topology optimization can attest to, a design obtained via topology optimization is rarely produced as-is. Consequently, following the topology optimization, design engineers create designs that try to capture as much as possible the material distribution of the optimal topology, but incorporate necessary aspects of the production process. This 'interpreted' design often incurs in a significant detriment of the structural performance. As a result, the engineer must spend a significant amount of time performing design changes by trial and error, which is not only a resource-intensive task, but one that leads to a suboptimal design. This is arguably the primary reason why topology optimization has not yet earned a permanent place in engineering workflows, despite the wide availability of commercial software.
One aspect that has been investigated in topology optimization (and is a focus of ongoing work) with regards to the physical realization of topology optimization designs is the incorporation of various 'manufacturing constraints.' Their purpose is to render designs that conform to certain aspects of manufacturing processes. For example, one such widely used constraint is the so-called 'casting draw' constraint, which renders designs that do not exhibit undercuts and thus require no cores for mold casting. The reader is referred to the survey in Liu and Ma (2016) for a review of manufacturing considerations in topology optimization.
In addition to ease of manufacturing, there are other production considerations that necessitate consideration of geometric requirements. These include, for example, the availability of high-level geometric features that facilitate positioning in fixtures, assembly and inspection. In many existing workflows for design of mechanical components (that do not use topology optimization), a widely used geometric representation that accommodates these requirements is the combination of geometric primitives via Boolean operations. Primitives have convenient datum entities (faces, edges, and vertices) that facilitate positioning in manufacturing fixtures and assembly. They also make it easier to specify dimensional tolerances, which is not only useful to perform tolerance stackup analysis but also for dimensional quality control and inspection. Thus, even if an organic design obtained via topology optimization can be manufactured by a given process, it still may be more convenient to replace it with a design made of primitives that accommodates the other aforementioned requirements. It should be noted that, in fact, most computer aided design systems for solid modeling work primarily by representing designs via primitives.
Despite the effectiveness of the aforementioned manufacturing constraints to enforce certain geometric requirements, it remains difficult to introduce general geometric constraints in the topology optimization. At the root of this difficulty lies the way in which density-based and level set-based topology optimization techniques represent the structure, namely by means of a density field or the level set of a function, respectively. These representations naturally and easily accommodate topological and shape changes, leading to organic designs. However, they make it difficult to impose geometric considerations such as, for example, obtaining a design with convenient datum entities for ease of production. Therefore, it is desirable to have topology optimization techniques that are endowed with a featurebased geometric representation of the structure that allows to impose these kinds of considerations, while retaining the advantage of existing techniques of employing a fixed grid for the analysis and optimization to circumvent remeshing.
Methods to incorporate geometric features with analytical geometry descriptions in topology optimization generally fall within two categories. First, there is a group of methods that, through various means, combine free-form topology optimization with embedded components or holes shaped as geometric primitives (cf. Chen et al. 2007;Qian and Ananthasuresh 2004;Xia et al. 2012;Zhang et al. 2011Zhang et al. , 2015Zhou and Wang 2013;Zhu et al. 2008). The second group, to which this work belongs, includes methods that represent the structure solely through a feature-based representation. The precursor of these methods is, in fact, one of the first topology optimization techniques: the bubble method (Eschenauer et al. 1994), where splines are used to represent the structural boundaries. This method requires re-meshing upon design changes, which limits its use. More recently, novel topology optimization techniques have been advanced that, as their density-based and level set-based counterparts, employ a fixed grid for the primal and sensitivity analyses, thereby circumventing the need for remeshing. Table 1 provides a comparison of these methods. The table includes only methods where the geometry is entirely composed of geometric primitives (to describe either the solid or the void regions).The geometric shapes listed for each method are the ones reported in the cited works; this does not imply these methods cannot potentially accommodate other geometric representations. Moreover, the analysis technique reported in the table is the one demonstrated in the references cited, however it can be argued that sharp interface methods (such as the extended finite element method) and ersatz material methods can be used interchangeably.
Among the methods listed in Table 1, two methods represent the structure as the union of primitive-shaped geometric components with fixed shape but variable dimensions and position: the moving morphable components method (Guo et al. 2014a(Guo et al. , 2016Zhang et al. 2016b;2017b), and the geometry projection method (Bell et al. 2012;Norato et al. 2015;Zhang et al. 2016a). The former method performs the topology optimization of a) 2d-structures, by employing as geometric components rectangles, approximated via superellipses, and quadrilateral shapes with two opposite sides described by polynomial and trigonometric curves; and b) 3d-structures, by using cuboids. The latter method uses rectangles with straight or semicircular ends for the design of 2d-structures, and plates with semicylindrical edges for the design of 3d-structures. The use of fixedshape components facilitates the design of structures made of stock material.
This work extends the geometry projection method to incorporate supershapes as a more flexible geometric representation to accommodate via a single formulation both designs of similar complexity to those obtained with free-form topology optimization, as well as designs made exclusively of primitives. While the designs presented in this work are still far from manufacturable as they do not incorporate many important geometric considerations, ease of production is the driving motivation, and our method aims to be a step in that direction by advancing the ability to represent primitive-only designs with a unified formulation.
The remainder of the paper is organized as follows. In Section 2, we introduce the supershapes. Section 3 describes the projection of the supershapes onto a density field for the analysis. The computation of the signed distance to Table 1 Existing methods for topology optimization with feature-based geometries. Abbreviations: topological derivative (TD), finite element analysis (FEA), extended finite element method (XFEM), topology optimization (TO), radial basis function (RBF) and non-uniform rational B-splines (NURBS) Method the boundary of a supershape, which is required for the geometry projection, is described in Section 2. Sections 5 and 6 describe the optimization problem and details of the computer implementation respectively. Several examples that demonstrate the proposed method are presented in Section 7. Finally, Section 8 draws conclusions of this work.

Supershapes
Supershapes (also called Gielis curves, Gielis 2017) are a generalization of superellipses that can exhibit variable symmetry as well as asymmetry and can describe through a single equation, the so-called superformula, a wide variety of shapes, including geometric primitives (Gielis 2003). The appeal of supershapes is their ability to express through a single equation a wide range of primitives and organic shapes. The superformula is given by (1) Figure 1 shows some examples of supershapes. When m = 4 and n 1 = n 2 = n 3 = n, this expression reduces to the equation for a superellipse (also known as a Lamé oval, Gridgeman 1970). Unlike superellipses, supershapes need not be symmetric; the parameter m controls the rotational symmetry. The values of a and b control the size, and the exponents n 1 , n 2 and n 3 control the curvature of the sides. The superformula can produce a wide range of shapes, including many shapes found in nature; the reader is referred to Gielis (2003;2017) for more examples.
For the purpose of this work, two modifications are made to the superformula of (1). First, the absolute value function is not differentiable at zero, which would preclude the use of efficient gradient-based optimization methods. Therefore, we replace it with the approximation Second, use of the power reduction formulas 2 cos 2 x = 1 + cos 2x and 2 sin 2 x = 1 − cos 2x provides a simplification both for computation and for the sensitivity analysis of the square trigonometric functions that result after replacing the absolute value with the above approximation. The modified superformula is given bỹ A Mathematica script to interactively plot a supershape with the foregoing expressions while modifying its parameters is presented in Appendix A. An important note about (4) and (5) is that the angle θ must be in the range [−π, π], as otherwise the supershape will not be closed. Therefore, if an angle is passed to these formulas, it must be appropriately wrapped around this range.
In addition to the foregoing parameters that define the size and shape of the supershape, we also wish to control its scale, location and orientation in space. Therefore, we introduce a scale factor s, the position x c of the origin of the coordinate system e 1 − e 2 on which (3) is applied, and the rotation φ = e 1 · e 1 of this coordinate system with respect to a global coordinate system e 1 − e 2 as additional design parameters, as shown in Fig. 2. We note that θ in (3) is defined with respect to e 1 . With this notation, a point on the boundary of the shape is given by where

Geometry projection
To perform the analysis on a fixed grid, we use an ersatz material, whereby the elasticity tensor at a point is modified as a function of a density variable to reflect a fully or partially solid material, or void. This strategy is the same used by all density-based and many level set-based topology optimization methods. The design variables in our method, however, are the supershape parameters described in the preceding section. Therefore, we need a mapping between the supershape parameters and the aforementioned density. This mapping must be differentiable so that we can obtain design sensitivities with respect to the supershape parameters and employ efficient gradient-based optimization methods. To this end, we use the geometry projection method (Norato et al. 2004;Bell et al. 2012;Norato et al. 2015;Zhang et al. 2016a). The idea of the geometry projection is simple: the density at a point p in space is the fraction of the volume of the ball B r p of radius r centered at p that intersects the solid geometry, i.e.: where ω denotes the region of space occupied by the structure and | · | is a measure of the volume (or area). If we make the assumption that r is small enough that the portion of the solid boundary that intersects the ball, ∂ω ∩ B r p , can be approximated by a straight line, then the density of (9) can be readily computed as a function of the signed distance d between p and ∂ω. As depicted in Fig. 3, in 2-d the density under this simplifying assumption corresponds to the fraction of the solid circular segment area to the area of the ball B r p , giving where the arguments of φ have been omitted on the righthand side of the equation for conciseness. Since r is fixed in the proposed method, it will be omitted as an argument hereafter. The vector of design parameters z for the entire structure is the aggregation of the vectors z i of design parameters for each of the supershapes, i.e. z = where N s is the number of supershapes. The signed distance d is positive if p is outside the shape, zero if p lies on its boundary, and negative if it lies inside. The computation of the signed distance is discussed in Section 4.
In addition to the design parameters that dictate the shape, scale, position and orientation of a supershape, we also ascribe a size variable α to each supershape. A unique feature of the method presented in Bell et al. (2012), Norato et al. (2015), and Zhang et al. (2016a) is that in the ersatz material this variable is penalized as in density-based topology optimization. Therefore, the elasticity tensor at a point p is given by where C is the ersatz elasticity tensor, C 0 is the elasticity tensor of the fully solid material,ρ is an effective density and q is a penalization power. The size variable α i of supershape i is part of its vector of design variables, i.e., This penalized size variable is an important ingredient of our formulation, as it allows the optimizer to entirely remove a geometric component from the design when the size variable becomes zero. Equation (10) is used for the calculation of the density with respect to a single geometric component. More generally, we consider a structure made by the union of multiple shapes. Here, this union is made through the maximum function, which corresponds to the Boolean union of implicit functions (Shapiro 2007). Specifically, at a point p, a composite density is defined as whereρ i and z i denote the effective density and design parameters for supershape i respectively. Since the maximum function is not differentiable, we replace it with a lower-bound Kreisselmeier-Steinhauser (KS) approximation (as in Zhang et al. 2017a): where k is a specified parameter and N s is the total number of shapes. Besides other benefits (described in Zhang et al. 2017a), this approximation has the advantage that it approximates the true maximum from below. This helps the method produce designs with less gaps between shapes in the optimal designs. These gaps occur because points in the gap regions obtain some stiffness from the nearby shapes through the geometry projection and thus have a stiffness that is appreciably larger than that of void regions away from the shapes. This phenomenon is exacerbated by the use of other aggregation functions (such as the pnorm) that approximate the maximum from above, as they consequently increase the effective stiffness of these gap regions. In our numerical experiments, the lower-bound KS function does a much better job of rendering designs without gaps. We note that if everyρ i for every supershape i is zero at a point p, the composite density will be consequently zero and the analysis will be ill-posed. Therefore, as is customary in ersatz material methods, a small lower bound ρ min must be imposed on the composite density. We achieve this viã where 0 < ρ min 1. Finally, the ersatz material is modeled using the composite density as In our method, we use the finite element method to solve the analysis problem. A composite density is computed at each finite element, with p above corresponding to the element centroid.

Signed distance calculation
Even though supershapes are given in polar coordinates, the computation of the signed distance is not as straightforward as in the case of the geometric representations used in our previous works, whereby the discrete geometric components are determined by a medial axis (in the case of approximately cylindrical bars, Norato et al. 2015) or a medial surface (in the case of plates, Zhang et al. 2016a). In these geometric representations, which in effect are implicit, distance-constrained offset surfaces (Bloomenthal 1990), the signed distance can be readily computed in closed form as the distance to the medial axis (or surface) minus the offset (e.g., the half-width of bars or half-thickness of plates). Due to the form of (1), however, it is not possible to obtain a closed-form expression for the signed distance from any point to the supershape boundary, and therefore it must be obtained numerically. It is worth noting that (6) constitutes an explicit representation, since this expression constitutes a rule to generate points on the boundary of the supershape, as opposed to implicit representations, whereby a rule determines whether a point is outside or inside the solid (Shapiro 2002).
Given a supershape with parameters z i , the distance from a point p to its boundary is given by where c, the point on the boundary of the supershape closest to p (cf. Fig. 2), is given from (6) by where θ * is the solution to the problem The foregoing problem can be solved numerically using root-finding methods to solve the first-order optimality condition f (θ) = 0. We employ a damped Newton method to solve this problem, with θ at iteration I given by: where the damping parameter is given by The Newton update is repeated until |f | falls under a specified tolerance restol. The derivative f can be analytically computed from 3) and (6) (cf. (48)). To compute f , we use forward finite differences. The Newton update can fail if f (θ) = 0, i.e. at an inflection point of the squared distance. However, we did not observe this situation in our numerical experiments. While a more sophisticated, robust and efficient solution strategy could be used to solve the shortest distance calculation (for example, a trustregion method), the above iteration is easy to vectorize for multithreading and works well in our examples.
A good initial guess θ 0 for the Newton iteration corresponds to the point c 0 where the line betwen p and x c intersects the boundary (shown in Fig. 2), which can be readily computed by A four-quadrant version of the arctan function must be employed and, as discussed in Section 2, the resulting angle must be wrapped to the interval [−π, π].
The sign of the signed distance can be computed, for instance, as Since we compute one composite density per element, a signed distance calculation to the element centroid is required for every element in the mesh. This calculation is evidently more expensive than the one used in our previous works, where a closed-form expression in terms of the design parameters is available. This cost is alleviated by the fact that these calculations are embarrassingly parallel, hence we can use parallel computing to speed up the computation. The computational cost of one signed distance calculation for the entire mesh is O(N el × N s ), where N el is the number of elements in the mesh. We note that this cost is linear in N el , whereas the cost of the finite element analysis is O(N q eq ), where N eq > N el is the number of equations and q > 1. Therefore, as the number of elements grows, the finite element analysis dominates the cost of the optimization as usual. Finally, it is possible (but not done in this work) to use efficient spatial partitioning techniques such as kD-trees to further reduce the cost of the signed distance calculation.

Optimization problem
In this work, we consider the classical topology optimization problem of minimizing compliance subject to a volume fraction constraint: In the above problem, C denotes the structural compliance. While we consider compliance minimization for simplicity in this work, we note that stress considerations can be incorporated in geometry projection methods (Zhang et al. 2017a). The set denotes the design envelope, i.e. the region of space where material can be distributed. We assume the tractions t are design-independent, and they are applied on a fixed portion of ∂ , which we denote by t . The volumes of the structure and the design envelope are denoted by V and | | respectively, and v * f is the constraint limit on the volume fraction. The volume of the structure is computed via Equilibrium of the structure is enforced for any design produced by the optimizer via the second constraint in (25), where a and l are the energy bilinear and load linear forms respectively, given by and U := {u ∈ H 1 ( ) : u = 0 on u } is the set of admissible displacements. The portion of ∂ on which homogeneous displacement boundary conditions are applied is denoted by u , and we assume this boundary to also be design-independent. The effective elasticity tensor C in (27) is computed from (17). Finally, in the above we have also assumed that there are no body loads applied on the structure. An important note is due here: for the SIMP penalization on the supershape size variables to be effective, the elasticity tensor C in (27) must be computed using a value of q in (12) larger than the value used to compute the volume in (26). Here, we use q = 3 for the former and q = 1 for the latter.
The last line in the problem of (25) corresponds to bounds on the design variables. The purpose of these bounds is discussed in detail in the following section. Since the problem is highly nonlinear, we impose move limits on the design variables so that only limited design steps are taken by the optimizer at each iteration, thus avoiding early overshooting to a poor design from which the optimization cannot recover. Since different variables represent different quantities and have different bound ranges, it is easier to impose move limits if we first scale all of the design variables to the range [0, 1], for example viaẑ i := (z i − z low )/(z upp − z low ), and then impose the same move limit on all of the variables. For a move limit 0 < M ≤ 1, we then replace the last line in (25) at iteration I with max(0,ẑ The sensitivities are accordingly scaled by ∂ẑ i /∂z i = 1/(z upp − z low ). Appendix B sumarizes the sensitivities of the functions in (25).

Computer implementation
We start the description of our computer implementation by providing an algorithm for our method, shown in Fig. 4. Our method is implemented in MATLAB. The finite element analysis is performed on a regular grid of N el square bilinear elements of size h. Each of the outer repeat iterations corresponds to a design update in the optimization, which is executed until the norm of the Karush-Kuhn-Tucker (KKT) condition or the relative change in the compliance between consecutive iterations fall below specified tolerances.

Signed distance and projected density calculation
We compute one projected density (10) per element, for which we compute the signed distance from the centroid of the element to each of the supershapes. This calculation corresponds to the first double for loop inside the repeat loop. This distance is computed by using the damped Newton iteration described in Section 4. The step size of the finite difference approximation of f (20) is 1E-6. The initial damping parameter is β (0) =1E-3, and β min =1E-6. The stopping criterion is restol =1E-8.
By using N el -vector representations of the angles in (20) and (23) and of the distance in (18) and (20), we take advantage of MATLAB's automatic multithreading of vector operations to parallelize the inner for loop. While this incurs in some redundant computation for elements for which the iteration converges faster than others, this strategy is still much more efficient than the alternative sequential solution. The computation of the sign of (24) is similarly Fig. 4 Algorithm for topology optimization with supershapes vectorized. We note that while in our implementation we use multithreading, the signed distance calculation is embarrassingly parallel and can be readily implemented in distributed memory programming.
For the calculation of the projected density, we employ a sample window radius of r = 2.5h in (10). In previous implementations of the geometry projection method we have used r = √ 2h/2, i.e., the sampling window circumscribes the element. However, in the numerical experiments performed in this work we found that a larger radius promotes faster convergence, particularly in early iterations. The calculation of the projected density of (10) is also done using efficient vector operations. An array of projected densities per element and per supershape is stored for subsequent calculations.

Finite element primal and sensitivity analyses
With the projected density for each element/supershape available, we proceed to solve the finite element analysis (corresponding to the second constraint in (25). To compute the stiffness matrix K j of element j , we compute the ersatz elasticity tensor of (17), which requires the calculation of the composite densityρ of (16). As mentioned before, we use q = 3 in the calculation of the effective density of (12) to penalize intermediate values of the size variable α i of supershape i.
The lower bound ρ min in (16) is set to 1E-4 in all presented results. Lower values increase the ill-conditioning of the stiffness matrix, and thus decrease the accuracy of the primal solution and, consequently, of the sensitivities. As ρ min decreases, the disagreement between the analytical sensitivities of the compliance and those approximated via finite differences increases. Inaccurate sensitivities hinder the convergence and efficacy of the gradient-based optimization. On the other hand, too high a ρ min leads to unrealistically stiff void regions. The foregoing value seems to consistently work well in our experiments, and is used throughout this work.
Finally, we use a value of k = 32 for the calculation of the lower KS function of (16).
The solution u to the primal analysis is used to compute the compliance C of (25). The computation of the volume of (26), on the other hand, is performed with q = 1 in (12). Finally, sensitivities of the compliance and volume fraction are computed following the expressions in Appendix B.

Optimization
Before taking the design step, we impose move limits on the design by modifying the lower and upper bounds on the scaled design variables as per (29). To perform the design update, we employ the method of moving asymptotes (MMA) (Svanberg 1987). Specifically, we use the MATLAB implementation corresponding to the version of MMA presented in Svanberg (2007). All MMA parameters are set to the default values suggested in Svanberg (2007) for the standard problem formulation, namely a 0 = 1, and a l = 0, c l = 1, 000 and d l = 1 for every constraint l in the optimization. The asymptote increase and decrease factors are set to the default values of 1.2 and 0.7 respectively. The reader is referred to Svanberg (2007) for a detailed explanation of these terms, and here we provide them solely for the purpose of reproducibility of our results.

Design variables bounds
An important aspect of our method that we discuss in this section is the need for appropriate bounds on the various design parameters. This need arises from considerations on 1) accuracy of the geometry projection, 2) robustness of the sensitivities calculation, and 3) physical realization.
With regards to accuracy of the geometry projection, we recall that the definition of the projected density of (10) rests on the assumption that the portion of the supershape boundary intersecting the sample window, ∂ω ∩ B r p , can be approximated as a single straight line. Therefore, we want to avoid shape boundaries that are too curvy in relation to the element size, as well as thin slivers that would lead to disjoint intersections. Loosely speaking, we want the element size h to be small in relation to the supershape. This is much easier to control in the offset surfaces reported in our previous work, as they accommodate (by design) very little shape variation.
To control the supershape proportions in relation to the element size h, we impose various bounds on its parameters. We set lower bounds on a and b to avoid shapes that resemble thin slivers. We also impose upper bounds on these parameters, since slivers can also develop for extreme a/b ratios if, for example, the exponents n 1 , n 2 and n 3 are such that convex sides form. With the same motivation, we impose lower bounds on these exponents to avoid highly convex sides, and, in some examples, we also set n 1 = n 2 = n 3 . In general, we aim to attain an adequate range for the aspect ratio of the shape. We note that it is difficult to strictly enforce the requirement that the intersection between the boundary of a shape and an element is approximately straight, since supershapes can develop corners. Nevertheless, this did not seem to be a problem in our numerical experiments given the applied bounds.
The second consideration as to why parameter bounds are needed pertains the calculation of sensitivities. As discussed in our previous work (cf., Norato et al. 2015;Zhang et al. 2016a), when the closest point on the boundary of the supershape is not unique, the sensitivities are not defined. In the case of offset surfaces, this occurs when the signed distance is evaluated at a point lying on the medial axis (in 2-d) or medial surface (in 3-d), and the sample window size is equal or greater to the bar width or plate thickness. In that case, it is easy to prevent this situation by requiring that the sample window size (and hence the element size) are sufficiently small with regards to these dimensions. In the case of supershapes, on the other hand, it is more difficult to prevent this situation, given the wide variation of shapes that can be accommodated. As before, controlling the aspect ratio of the supershape via bounds on the aforementioned parameters seems to effectively avoid this situation in practice.
Non-uniqueness of the signed distance can also arise from sharp corners in the supershapes. However, the substitution of the absolute value in the supershape formula with the approximation of (2) rounds the corners. We use = 1E-3 in (2) for all of our examples. We also place upper bounds on the exponents n 1 , n 2 and n 3 to avoid increasing the non-linearity of the signed distance in the supershape parameters, and consequently of the optimization problem functions.
The third and final consideration to put bounds on the supershape parameters is simply physical realization. In particular, we want to avoid shapes so small that they could not be fabricated. The most natural parameters to control the shape size are a and b along with the scaling factor s. In general, however, enforcement of a strict minimum size in the structure is not straightforward. 1 Convex sides can render a supershape size smaller than a and b. Moreover, Fig. 5 Short cantilever beam problem the intersection of supershapes can be smaller than the minimum size of the intersecting supershapes. A strategy to enforce a minimum size in geometry projection methods is described in Hoang and Jang (2016). We defer investigation of this issue to future work.
In summary, an important aspect of using supershapes, given the wide freedom of shapes they can represent, is the imposition of bounds on the parameters to ensure robust behavior of the proposed method. The values of the bounds presented in our examples seemed to consistently work well. However, a more in-depth understanding of the various interplays between supershape parameters would render a more systematic strategy to control the resulting supershapes. Such a study is deferred to future work.

Results
In this section we present numerical examples to demonstrate the proposed method. In all examples, we consider an isotropic elastic material with Young's modulus E = 1E5 and Poisson's ratio of ν = 0.3 (since minimum-compliance designs depend only on the volume fraction and not on the magnitude of the applied load or the material's modulus, we henceforth omit units for brevity). The magnitude of the applied load in all examples is F = 10. All supershapes have a thickness of t = 1, and in all of the initial designs we assign α q = 0.5 to all the shapes. Unless noted, in all examples, we use a move limit of m = 0.05.

Short cantilever beam
We first design a short, cantilever beam (cf. Fig. 5) with a volume fraction constraint of v * f = 0.5. A 40 × 40 grid is used for the analysis, i.e., h = 0.05. In this example, we use a single component, with the intention of showing the flexibility of supershapes. Therefore, this constitutes a shape optimization example. Geometry projection methods have been used for shape optimization in Norato et al. (2004;Wein and Stingl 2018).
The supershape plot is clipped to the design region , since this is the only portion of the supershape that has an effect in the compliance and the volume (correspondingly, the fabrication of the structure should cut material outside  Fig. 7 Composite densityρ for optimal short cantilever beam the design region). The composite densityρ for the optimal design is shown in Fig. 7; note that in this case the composite density is the same as the effective densityρ, since there is only one shape in the design. The optimal design is also shown in Fig. 8 without clipping to . In the same figure, we compare the optimal supershape design to the shape corresponding to uniform energy density along the boundary, and with volume fraction of 0.5. For an elastic structure with a single designindependent loading, and as long as the design region allows it, the minimum compliance optimal shape has uniform strain energy density along the boundary (Pedersen 2000). Moreover, under the same conditions, the stiffest design is also the strongest design. Thus, the uniform-energy shape is the same as the uniform-strength shape, which Galileo found to be parabolic for a cantilever beam (Timoshenko 1953). Using Euler beam theory, it can be readily shown that Fig. 8 Optimal supershape for short cantilever beam design (continuous line), without clipping to design region . Dashed line corresponds to the uniform-energy boundary for 0.5 volume fraction the uniform-energy shape for this design region and for a volume fraction of 0.5 is given by The above expression is based on an origin (x, y) = (0, 0) at the midpoint of the left side of the design region. As seen in Fig. 8, there is good agreement with the uniformenergy shape.

MBB beam
Our next example corresponds to the Messerschmitt-Bölkow-Blohm (MBB) beam, extensively studied in topology optimization. The problem setup is shown in Fig. 9. The volume fraction constraint for this problem v * f = 0.4. A 160 × 40 grid is used for the analysis, i.e., h = 0.125. Taking advantage of the symmetry of the problem, we perform the analysis and design on half of the design region.
The bounds for this problem are a, b ∈ [0.4, 5], m ∈ [2, 6], s ∈ [1, 2], x c x ∈ [0, 20], x c y ∈ [0, 5] and φ ∈ [−π, π]. In these examples, we assign n 1 = n 2 = n 3 = 10 and fix these parameters in the optimization to avoid unfavorable, high aspect ratio and/or highly concave shapes that could cause the geometry projection to fail. The convergence tolerances are kktol = 1E-3 and obj tol = 1E-5. We solve the topology optimization using different initial designs with supershapes arranged in an N x ×N y grid, with positions x c x and x c y evenly spaced by 20/(N x + 1) and 5/(N y + 1) respectively. The remaining parameters for the shapes in the initial design are a = b = 0.5, m = 4, s = 1.0, and φ = 0. Figure 10 shows the initial and optimal designs for beams with 5 × 2, 10 × 2 and 10 × 4 supershapes. The optimization converges in 148, 131 and 152 iterations, and the compliance values of the optimal designs are 0.5938, 0.5439, and 0.5006 respectively. The volume fraction of the optimal design in all cases is less than 1E-3 from the imposed limit of 0.4. We first observe that the optimal designs are reminiscent of known density-based optimal topologies for the MBB beam. Although there are differences in the designs obtained, the bounds on the supershape parameters clearly impose a minimum size on the design because we do not see thinner members in one design than in the others. We note that the optimizer removes some shapes in the 10 × 2 and 10 × 4 designs.
This example illustrates one aspect of our formulation, and that is that in some cases there are gaps between the shapes. As discussed in Section 3, this owes primarily to the fact that these gaps are much stiffer than void regions that are away from the shapes. Since the composite density in these gaps is in general closer to unity than it is to zero as a result of the geometry projection with respect to the nearby shapes, and since this composite density is reflected in the analysis via (17), it is not unreasonable to close these gaps in the actual design. This can be easily done by modifying the primitives around the gap.
We also use this example to discuss the effect of mesh size on the results. As discussed in Section 6.4, we require the element size to be smaller than the minimum size of the shape to ensure well-defined sensitivities. Beyond this requirement, the mesh can be subsequently refined.
To illustrate the effect of mesh refinement, we perform the optimization for the MBB beam with 5 × 2 shapes (with initial design shown in Fig. 10, top/left) and with finite element meshes of 80 × 20, 160 × 40 and 200 ×    Fig. 10, top). The results are shown in Fig. 11. The optimization converges in 250, 148 and 201 iterations, and the compliance values of the optimal designs are 0.5283, 0.5439, and 0.5714 respectively.
Two observations are worth making. First, a minimum length scale is in general satisfied, and this is dictated by the fixed number of shapes, and by the lower bounds on their dimensions. That is, regardless of the mesh size, we expect to obtain designs of similar complexity if we use the same design representation (number of shapes and bounds on shape parameters). Second, as we expect, a finer mesh provides a better resolution of the shapes upon projection, and therefore it provides better control over the shape boundaries than coarser meshes. The compliance values actually increase when refining the mesh, but part (if not all) of that increase is due to the fact that the finer mesh is more compliant, since it has a larger number of analysis degrees of freedom.

Long cantilever beam
The designs obtained in the preceding section are reminiscent of MBB beam designs obtained with free-form topology optimization methods. This is an indication that the proposed method obtains expected results. However, these designs can be more easily and more efficiently obtained using, for example, density-based topology optimization. Furthermore, one could post-process the optimal free-form design and 'fit' a set of supershapes to it, as it has been done in the past using, e.g., splines (cf. Bendsøe and Sigmund 2003 and the extensive references therein).
Where supershapes provide interesting possibilities is in producing designs made of geometric primitives such as, for example, rectangles, ellipses and triangles. As explained in Section 1, the availability of datum geometric entities helps designs with primitives facilitate certain aspects of the production of mechanical components that are not available in free-form designs. To demonstrate the design of structures made of primitives with our method, in this section we design a long cantilever beam with the same dimensions and analysis mesh as the beam of the preceding section, but with the load and boundary conditions shown in Fig. 12. The volume fraction constraint for this problem is again v * f = 0.4. All examples employ a 5 × 2 grid of supershapes in the initial design, and the convergence tolerances are kktol = 1E-3 and obj tol = 1E-4. In all examples, a = b = 0.5, φ = 0, s = 1 and α = 0.5 for all supershapes in the initial design. Likewise, the bounds φ ∈ [−π, π], s ∈ [1, 3], x c x ∈ [0, 20] and x c y ∈ [0, 5] are applied to all examples. Figure 13 shows the initial and optimal designs for a beam whose initial design is made of triangles, ellipses, rectangles, and a combination of ellipses and rectangles. For triangles, we fix m = 3, n 1 = 4, and n 2 = n 3 = 8; for ellipses, we fix m = 4, n 1 = n 2 = n 3 = 2; and for rectangles, we fix m = 4, n 1 = n 2 = n 3 = 8 (note that rectangles are effectively modeled using superellipses). We note that in the case of combined ellipses and rectangles, the foregoing parameters for each supershape are fixed throughout the optimization, i.e. a supershape that is an ellipse or a rectangle in the initial design remains an ellipse or rectangle throughout the optimization. In all cases except Optimal designs for cantilever beam problem using shapes that morph into ellipses or rectangles with different initial values of the parameter n = n 1 = n 2 = n 3 . The color of the shapes (corresponding to the horizontal color bar) denotes their size variable α. Supershapes are clipped to the design region the triangles, we impose the bounds a, b ∈ [0.4, 5]. In the case of triangles, it is more difficult to capture a nonequilateral triangle with the supershapes, and larger a/b aspect ratios tend to render shapes that look more like water drops, and so we impose a smaller range of a, b ∈ [0.4, 4]. Also in the case of triangles, we fix the value of b = 0.4, because for b < a one obtains 'bowtie' shapes. Note that we can still get triangles with size larger than 0.4 thanks to the scale factor s. From top to bottom in Fig. 13, the compliance values for the optimal designs are C = 0.5759, 0.5263, 0.566 and 0.5348. These designs are attained in 150, 127, 129 and 141 iterations respectively. All designs satisfy the volume fraction constraint within 1E-3. Clearly, the design that performs the worst is the one with triangles. We posit this is because the geometry of the better designs cannot easily be captured with ten triangles.
An even more interesting possibility is to let shapes morph into a selected set of primitives. 2 We illustrate this idea by repeating the cantilever beam design, while forcing shapes to morph into either rectangles or ellipses. Both shapes have the same symmetry (m = 4), but different exponents (n = 2 for ellipses, n = 8 for rectangles, n 1 = n 2 = n 3 = n). In the optimization, we fix m, but we let n vary within [2,8]. To steer shapes towards becoming ellipses or rectangles, we add the following 'shape' constraint to the optimization problem of (25): In the above expression, n i is the value of n for shape i. Clearly, if n i = 2 or n i = 8, the contribution of shape i to the sum is zero. The number 81 in the denominator ensures that the value in the sum is at most 1.0 (which occurs when n i = 5). The factor outside the sum helps ensure that g s ∈ [0, 1], hence the choice of the constraint limit g * s can be made independent of the number of shapes. In the following examples, we use g * s =1E-3. To avoid the shapes from 'locking' prematurely into the desired primitives, we use a continuation strategy, whereby we set g * s = 1 in the first iteration, and then decrease it by 2E-3 in each iteration until we reach the desired value of g * s =1E-3 (that is, it takes at least 50 iterations to attain the desired constraint limit). In addition to this continuation strategy, for this example we employ a tighter move limit of M = 0.025. Figure 14 shows several design iterates using a starting design with n = n 1 = n 2 = n 3 = 4. The optimal design is made of two rectangles (n ≥ 7.99) and six ellipses (n ≤ 2.124). The figure shows how the shapes morph into one of the two primitives as the optimization progresses.
This example also serves to illustrate a characteristic of the proposed method, namely that it is more prone to converging to unfavorable local minima than free-form topology optimization methods. This is a trait of geometry projection methods we have observed before (Norato et al. 2015;Zhang et al. 2017a), and it is likely due to the more restrictive representation of possible geometries imposed by the discrete geometric components. To illustrate this dependency, we show in Fig. 15 the results of the same optimization problem with different initial values of n. Table 2 lists results for the corresponding optimal designs, including compliance values, number of iterations to convergence, and number of rectangles (n ≈ 8) and ellipses (n ≈ 2) with their respective ranges for n i . All designs satisfy the volume fraction and shape constraints tightly.

Conclusions
This paper demonstrates a method to perform topology optimization using supershapes. The numerical examples show the method effectively produces designs similar to those obtained with free-form topology optimization techniques. The appealing feature of the proposed method is the ability to produce designs made exclusively of various geometric primitives that are all modeled with a single representation, the superformula. The fact that different primitives can be represented via a single equation also allows to produce designs whereby shapes morph into specified types of primitives. We demonstrated this capability by introducing a shape constraint that steers the shapes into becoming one of two types of primitives (rectangles and ellipses).
To fulfill its final purpose, which is to produce designs made of primitives that are amenable to production, the proposed method requires further advancement. The extension to 3-dimensional problems is the most immediate one. The strict imposition of minimum size constraints, which in this work are loosely enforced by imposing bounds on the supershape parameters, is also an important need. A deeper understanding of the interactions between supershape parameters is needed in order to facilitate the imposition of manufacturing-driven geometric constraints. Finally, strategies to circumvent convergence to an undesired local minimum are necessary. These strategies are not only important from the point of view of how good the optimal design is, but also to free the designer from having to try different initial designs. These research directions are a matter of ongoing and future work. of φ. To cirumvent this, one would have to make the plotting range a function of φ. An easier solution is to use the ParametricPlot function. The entire script is as follows, which produces an interactive plot like the one shown in Fig. 16. r[t_, m_, a_, b_, n1_, n2_, n3_, s_, eps_]:= s ((eps+(1/(2 aˆ2)) (1+Cos[m t/2]))ˆ(n2/2)+ (eps+(1/(2 bˆ2))(1-os[m t/2]))ˆ(n3/2))( -1/n1) x [t_,m_,a_,b_,n1_,n2_,n3_, eps_, phi_, s_] Figure 16 shows the parameter values for the optimal design of the example in Section 7.1. The plotting range can be adjusted for each one of the parameters in the last statement of the script. If, say, n = n 1 , n 2 , n 3 , then the last statement should change to Manipulate[ ParametricPlot[{x[t,m,a,b,n,n,n,e,phi,s],y[t,m,a,b,n1,n2,n3,e,phi,s]

Appendix B: Sensitivity analysis
In this Appendix, we list the expressions necessary to obtain design sensitivities necessary for the gradientbased optimization. For brevity, we do not derive these expressions; however, they can be readily obtained. We start by stating the sensitivities for the geometry projection and the signed distance, and then state sensitivities for the optimization functions. where D zρj , U j and K j 0 are the sensitivity of the composite density, the vector of nodal displacements, and the fully solid stiffness matrix of element j respectively.
The sensitivity of the volume of (26) is similarly given by and can be simply computed in the finite element discretization as D zρj v j , where v j is the volume of element j . Finally, the sensitivity of the shape constraint of (31) is simply