A MATLAB code for topology optimization using the geometry projection method

This work introduces a MATLAB code to perform the topology optimization of structures made of bars using the geometry projection method. The primary objective of this code is to make available to the structural optimization community a simple implementation of the geometry projection method that illustrates the formulation and makes it possible to easily and efficiently reproduce results. A guiding principle in writing the code is modularity, so that researchers can easily modify the program for their own purposes. Another goal is efficiency, for which extensive use of vectorization is made. This paper details the formulation of the geometry projection, discusses implementation aspects of the code, and demonstrates some of its capabilities by presenting several 2D and 3D compliance minimization examples.


Introduction
The prevalent techniques to perform topology optimization of continua are the density-based and the level-set methods (Sigmund and Maute 2013).These techniques produce organic designs that are highly efficient.In some cases, a design that closely follows the optimal topology can be manufactured using, for example, additive manufacturing or casting techniques.However, when the most economical fabrication process consists of joining stock material such as bars or plates, it can be very difficult to translate the optimal topology into a design that conforms to that process.This difficulty has motivated the development of topology optimization techniques that produce designs exclusively made of geometric primitives.
One of the topology optimization techniques that render designs made of geometric components is the geometry projection method (GPM) (Bell et al. 2012;Norato et al. 2004Norato et al. , 2015)).There exist other techniques to perform topology optimization using geometric components, such as the moving morphable components method (Guo et al. 2014;Zhang et al. 2016b); a review of these techniques is outside of the scope of this paper; and we refer the reader to the recent review by Wein et al. (2019).The purpose of this work is to introduce a MATLAB code to illustrate the GPM for the topology optimization of 2D and 3D structures made of bars.The GPM is described in detail in Section 2. In particular, we aim to demonstrate how the geometry mapping can be performed in an efficient manner using vectorized operations.Unlike other educational codes published in this journal, we do not attempt to fit the code into a relatively low number of lines and include it in the manuscript.Although this would be in principle possible, we believe in our case it would sacrifice clarity and therefore hinder the objective of explaining the inner workings of the geometry projection.Instead, this manuscript provides an overview of the program, and the code is released for free and made available through GitHub at https:// github.com/jnorato/GPTO.This approach allows us to better modularize the code but also to add various functionalities and options that users can experiment with, which we believe will be beneficial to the research community.The code is released under a Creative Commons CC-BY-NC 4.0 license, which means it is free for non-commercial use and that appropriate credit must be given.The program is not an open source program in the sense that a repository to which users can contribute to further development of the code is not available.Nevertheless, the authors welcome any suggestions that users may have for future improvements.
The rest of the manuscript is structured as follows.Section 2 presents formulation of the geometry projection method.In Section 3, details of the implementation are provided.Some usage examples are presented in Section 4, and we draw conclusions of this work in Section 5.

Geometry projection formulation
The basic idea of the geometry projection is to take a highlevel parametric description of a geometric component ω b and map it onto a pseudo-density field over a design region Ω ⊇ ω b .This field is subsequently discretized via a fixed mesh for analysis.The mapping must be differentiable so that design sensitivities of the optimization functions are well defined, and thus efficient gradient-based optimizers can be employed.This section describes the formulation of the geometry projection and of the sensitivities.

Projected, penalized, and combined densities
We first consider the projection of a single geometric component.The projected density at a point x is the volume fraction of the intersection between the ball B r x of radius r centered at x and the component ω b (Norato et al. 2004), i.e., where z b is the vector of geometric parameters that describe ω b .For arbitrary shapes of the component ω b , computing exactly this intersection is not straightforward.Therefore, we seek an approximation of (1) that is easy to compute and is differentiable.If we assume ω b to be a smooth closed manifold and r to be sufficiently small, then ∂ω b ∩B r x can be approximated as a line segment in 2D and a circle in 3D (since B r x is a circle and a sphere, respectively).With that assumption, the projected density can be approximated as the area fraction of the circular segment of height h = r − ϕ b in 2D (or volume fraction of the spherical cap of the same height in 3D), where ϕ b (x, z b ) is the signed distance from x to ∂ω b (see Fig. 1). 1 The projected density is thus given by where The notation e H is used here to bring attention to the fact that this is a regularized Heaviside; however, we emphasize that it corresponds to the area (volume) fraction calculation of the circular (spherical) cap.The size r of the sample window is fixed throughout the optimization; hence, it is hereafter omitted as an argument.We also omit arguments on the right-hand side of equations for brevity.The calculation of ϕ b depends on the particular representation of the geometric components.In the case of the offset bars considered in this work, it will be detailed in Section 2.2.
A key ingredient of the geometry projection method is that, in addition to the parameters that describe the shape of each component, a size variable α b ∈ [0, 1] is ascribed to each component b.This size variable is penalized in the same spirit of penalization schemes employed in density-based topology optimization, so that a value of α b = 1 indicates the geometric component must be part of the structure, while a value α b = 0 indicates the component must be removed from the design.This feature makes it easier for the optimizer to remove geometric components and modify the topology.The penalization is achieved via the definition of a penalized density b ρ b that incorporates the size variable α b as where μ is a penalization function, q is the penalization parameter, and we note that α b ∈ z b .For example, for a power law penalization, as in the solid isotropic material penalization (SIMP) commonly used in density-based topology optimization, we have that μ(α b ρ b , q) = (α b ρ b ) q .Importantly, for the penalization to be effective, it is required that the volume of the structure be computed with the unpenalized density (e.g., q = 1 for SIMP).
Here we note that, unlike our previous works, the projected density ρ b is also penalized, as otherwise the material interpolation is linear in ρ b when α b = 1, which renders a physically unrealistic material wherever intermediate-density material appears (i.e., along the boundaries of geometric components), as demonstrated in Bendsøe and Sigmund (1999).We also observe slightly better convergence with this modification, and this strategy produces less gaps in between components in the final design that are likely due to unrealistically stiff gray regions.
When multiple bars overlap, we combine the penalized densities for all bars into a combined density given by where n b is the number of geometric components, g max denotes a smooth approximation of the maximum function, Z≔ is the vector of design variables for all components, and 0 < ρ min ≪ 1 is a positive lower bound to prevent an ill-posed analysis.An example of a function that embodies (5) is the modified p-norm described in Zhang et al. (2016a):1 regardless of the number of geometric components.Finally, the combined density is reflected in the analysis by using an ersatz material, with the elasticity tensor modified as where ℂ 0 is the elasticity tensor for the material the geometric components are made of.
The foregoing formulation can be used for geometric components of any shape that are made of a single, isotropic material, so long as a signed distance and its derivatives with respect to the geometric parameters can be computed.For instance, this scheme has been employed to design structures made of bars (Norato et al. 2015), flat plates (Zhang et al. 2016a), curved plates bent along a circular arc (Zhang et al. 2018), and supershapes (Norato 2018), which are a generalization of hyperellipses with variable symmetry.The case of multi-material structures, where components can be made of one out of a given set of isotropic materials, and where the optimizer simultaneously determines the optimal components layout and the best material for each component, requires a different strategy to combine the geometric primitives and to interpolate the properties of the different materials (cf.Watts and Tortorelli (2017) and Kazemi et al. (2018)).

Distance
We now describe the computation of the signed distance ϕ b in 1.Here, as in some of our previous works, we represent bars as offset surfaces (Norato et al. 2015).Specifically, the boundary of the bar is given by the set of all points at a distance r b of the line segment with endpoints x 1b and x 2b (cf. Figure 1).This definition represents bars as rectangles with semicircular caps in 2D and cylinders with semispherical caps in 3D.The user must define an initial design made of a set of bars; bars can be removed from the design during the optimization (e.g., by setting its size variable to zero) or reintroduced (by increasing a zero size variable), but in this code, bars cannot be introduced during the optimization.
The vector of design parameters for bar b in the offset surface representation is therefore given by z b = {x 1b , x 2b , r b , α b }.Our formulation allows for the case where endpoints are shared by two or more bars, thus they remain "connected" at these endpoints throughout the optimization (much like ground structure approaches).The total number of design variables is 2(n d n p + n b ), where n d = 2, 3 in 2D and 3D, respectively, and n p is the total number of endpoints.Alternatively, a bar can be defined to have its own medial axis endpoints, thus it can "float" within the design region independently of other bars.If all bars are "floating," then n p = 2n b .
The advantage of using the offset surface representation is that the signed distance to the boundary of the bar can simply be computed as the distance to the medial axis (d b ) minus the bar radius as To compute the distance to the bar's boundary, it is therefore only necessary to compute the distance to the bar's medial axis.Moreover, as we show in the sequel, the distance d b and its design sensitivities can be computed in closed form as a function of the bar's design parameters z b .
Though not strictly necessary, an element-uniform combined density (5) is employed for the computation of the ersatz material of (7).Therefore, the combined density and consequently the signed distance are computed at the element centroid, which we denote as x e (cf. Figure 2).We also use the notation The length of bar b is given by and the unit vector along the bar's medial axis is given by Consider a Cartesian coordinate system for each bar in which x 1b is the origin and a b is the first axis.On this coordinate system, ℓ be and r be are, respectively, the parallel and perpendicular components of x e/1b given by The distance from the medial segment of bar b to the centroid of element e is given by

Optimization problem
The optimization problem solved by default in our code is the compliance minimization problem given by where c is the compliance, v f is the volume fraction and v * f its limit, and ρ is the combined density of (5).The domain Ω denotes the region occupied by the design envelope, and Γ t and Γ u denote the portions of ∂Ω where traction t and displacement boundary conditions u are imposed, respectively, with Γ t ∪ Γ u = ∂Ω, Γ t ∩ Γ u = ∅ as usual.We assume there are no body loads and the applied traction is design-independent.The admissible sets for the trial displacement functions u and the test functions v are given by U≔ u u∈H , respectively.The energy bilinear form a and the load linear form l in ( 17) are computed as with ℂ the ersatz elasticity tensor of ( 7), ∇ S the symmetric gradient operator and

Sensitivities
In this section, we present the expressions required to compute analytical sensitivities of the compliance and volume fraction.
For an optimization function θ that depends on the design exclusively and implicitly through the displacements, its sensitivity with respect to a design variable z i is given by where n e is the number of elements.Using adjoint analysis, and as customary in density-based topology optimization, for the compliance (i.e., θ ≡ c) we have where u e is the vector of nodal displacements for element e, K e is the element stiffness matrix computed using the ersatz material of (7), and K e 0 is the "fully solid" element stiffnes matrix corresponding to ρ e = 1.It should be noted that K e 0 only needs to be computed once and for all in the optimization and stored in memory.
The sensitivity of the combined density, ∂ρ e /∂z i in (23), is obtained from (5) as where b ρ be is the penalized density for element e and bar b of (4).It should be noted that ∂b ρ be =∂z i is zero when z i does not correspond to one of the design variables of bar b.The difference between "floating" and "connected" bars is that in the former case ∂b ρ be =∂z i is non-zero for (at most) one bar and only the summand for the corresponding bar needs to be computed.The actual expression for the derivative of the smooth maximum depends of course on the particular function chosen.
The sensitivities of the penalized density b ρ be are obtained from (4) as If z i corresponds to one of the components of x 1b or x 2b , then ∂α b /∂z i = 0 and the second term on the right-hand side of (26) vanishes.Conversely, if z i ≡ α b , then ∂ρ be /∂z i = 0 and ∂α b / ∂z i = 1, thus the first term on the right-hand side of ( 26) vanishes.Finally, if z i does not correspond to one of the design parameters of bar b, then the entire expression is zero.Consequently, each of the terms in the right-hand side of ( 26) is computed separately for each bar, and the sensitivities are subsequently "assembled" into a matrix.
From ( 2), the sensitivity of the projected density ρ be is given by This expression highlights the fact that the sensitivities of the projected density are non-zero only in a band of width 2r along the bar boundary.From (3), we find with x = ϕ be /r.The sensitivities of the signed distance of (8) are given by Finally, the sensitivities of the distance d be with respect to the bar endpoints x sb , s ∈ {1, 2} are given by where δ s k ≔ 1 if s ¼ k; 0 otherwise f g is the Kronecker delta.It can be noted from (30) that ∂d be /∂x sb is undefined when d be = 0 (Norato et al. 2015;Zhang et al. 2016a).In this case, an element centroid x e lies exactly on the medial axis of bar b.This situation can be circumvented by making sure the sample window size is smaller than the bar's width, i.e., r < r b , since an infinitesimal perturbation of the bar's medial axis endpoints will leave the projected density of any point on the medial axis unchanged (ρ be = 1), and thus ∂ρ be /∂x sb in (27) must necessarily be zero.
When the bar's axis collapses to a point (so that ℓ b = 0), the sensitivities are still well defined, since the first branch of (30) applies.However, since the code still computes a b /ℓ b , we check ℓ b against a small tolerance to prevent a division by zero.
The sensitivities of the volume fraction of (16) (i.e., θ ≡ v f ) are computed as where v e is the volume of element e and V ¼ ∑ n e e¼1 v e is the volume of the design region.The term ∂ρ e /∂z i is computed as before using (25-30), however with the clarification that the volume fraction must be computed with the unpenalized density.Consequently, ∂μ/∂(α b ρ be ) = 1 in (26).
Most of the equation numbers for this section have been added to comments in the code so that the user can readily find them (for instance, as shown in Appendix A).

Algorithm
We complete the presentation of the geometry projection formulation with the pseudo-code shown in Fig. 3.This algorithm presents a conceptual (but logically correct) flow of the steps taken to perform the optimization.However, as will be mentioned in the following section, an efficient implementation of the code does not exactly follow this pseudo-code.In Fig. 3, b z denotes the vector of design variables after they have been scaled so that they lie in the range [0, 1].This scaling and the imposition of move limits on each design update are discussed in detail in Section 3.5.

Implementation
In this section, we discuss general aspects of the implementation, particularly in terms of functionality of the code.We also discuss in some detail the calculation of the geometry projection described in the previous section.The code was developed and tested using MATLAB, version R2018b.It was tested in the Red Hat Enterprise Linux 7.4, macOS Mojave, and Windows 10 operative systems.Because we make extensive use of vectorization, the code performs better in multi-core machines.The code is executed by running the script GPTO.m (without any arguments).To check that the program is running, we suggest the user simply runs this script, which executes the optimization of a 2D cantilever beam (not described in this work).
It should be noted that performance will suffer drastically on non-Intel CPUs because Intel's math kernel library (MKL) disables by default all but the most basic vectorization on non-Intel CPUs.However, this can be circumvented by setting an environment variable: MKL_DEBUG_CPU_TYPE = 5.

Organization
Since we do not attempt to write the code in as concise a manner as possible, we are free to modularize it as much as possible, which we have attempted to do.Similar functions are placed under the same subfolder-for instance, those related to the finite element analysis, geometry projection, optimization, plotting, etc.This allows the user to more easily navigate the code, examine and focus on particular portions of the implementation, and make changes.We have adhered to the rule that one routine corresponds to one MATLAB script (i.e., no script has multiple functions declared in it).The modular structure also facilitates having multiple options for, e.g., the penalization function μ of (4), the aggregation function of (5), the mesh input or generation, etc.The root folder only contains the main script GPTO.m and another script to invoke the input functions (get_inputs.m);all other scripts are inside subfolders.
Another important aspect of our implementation is that we use three MATLAB structures, named FE, GEOM, and OPT, to store all information related to the finite element analysis, geometric components, and method parameters, respectively.These structures are declared as global variables in any script where they are needed.Having only a few structures and declaring them as global variables avoids having to pass long lists of arguments to the functions that need the information stored in them.If a new field is added to one of these structures, it becomes immediately available to any routine invoking the structure as a global variable.Using global variables is in general discouraged in modern compiled languages because they may not be thread-safe; however, since we are only using MATLAB's own vectorized operations, we assume that potential problems such as race conditions are managed by MATLAB (and we have to date not observed any issue related to this in our numerical experiments).Also, the use of global variables makes the code slightly more efficient, since local copies of the structures (some of which may be relatively large, such as FE) need not be created in each routine that uses them.

Inputs
The user must provide all of the inputs using MATLAB scripts.This strategy circumvents having to parse text files, allows in providing inputs in any order, facilitates the incorporation of comments, and makes it easier for the user to customize the code.Although not strictly necessary, we have placed all the input files for the sample problems into the subfolder \texttt (in fact, we created additional subfolders inside this subfolder for each one of the examples).To switch from running one example to another, the user needs only to update the single line with the run command in the inputs.mscript located in the root folder to indicate the location of the master input file for the run.Having a script that calls the master input files makes it easy to keep different inputs for different runs.Note that the \input_files subfolder is not added to the search path out of precaution to avoid potential name conflicts and because it is not necessary for the code to run.
Since there are quite a few options in the program, the easiest way to create input files for a new run is to copy and modify the files for one of the sample problems.The input files for all problems are well commented to make it easy to modify them.All of the inputs are used to initialize the aforementioned three data structures that pertain to the finite element analysis, the bars that make up the initial design, and all parameters related to the geometry projection, the optimization problem, the optimizer, and the finite difference check of sensitivities (if requested).

Finite element analysis
The parameters relating to the finite element mesh are placed in the FE.mesh_input structure.Three options are available to create or read a mesh, which are indicated in the type field.The 'generate' option creates a mesh on the fly for rectangular Fig. 3 Algorithm for topology optimization using geometry projection and cuboid design regions in 2D and 3D, respectively.The user must specify the dimensions of the region and the number of elements along each dimension using the box_dimensions and elements_per_side vectors, respectively.The second option, 'read-home-made', can be used to read a mesh that has been previously created (e.g., using the makemesh script provided in the code) and saved to a MATLAB .matfile, whose path must be provided in the field mesh_filename.The option 'read-gmsh' allows the user to read a mesh created with the open source software Gmsh (Geuzaine and Remacle 2009).The user must create a mesh made of quadrilateral or hexahedral elements in 2D or 3D, respectively (a transfinite mesh in Gmsh parlance); export it to MATLAB format; and indicate the path of this file in the gmsh_filename field.This third option is very convenient for design regions that are not cuboid-shaped.
The mesh only contains the node locations and the element connectivity.Displacement boundary conditions and loads must be specified in a separate MATLAB script that the user must create and whose path must be specified in the bcs_file field.This separation facilitates using the same mesh for different problems with different boundary conditions.As before, it is easier to copy the file for one of the sample problems a n d m o d i f y i t .F o r c u b o i d -s h a p e d m e s h e s , t h e compute_predefined_node_sets function provides a very convenient utility to retrieve node sets corresponding to prescribed points, edges, and faces of the domain, which can be subsequently used to impose displacement boundary conditions or to apply loads.These node sets are stored in the FE.node_set structure.For example, for a 2D rectangular domain, the set BR_PT contains the node in the bottom-right corner of the domain, and the set L_edge contains all nodes on the left-hand side edge of the domain.Needless to say, the user must ensure that the problem is sufficiently restrained to get a well-posed analysis.
The code for the finite element analysis closely follows the sparse data structures and vectorized operations presented in Andreassen et al. (2011).It provides the option to use a direct o r a n i t e r a t i v e s o l v e r t h r o u g h t h e p a r a m e t e r FE.analysis.solver.type.In the former case, we simply use MATLAB's matrix left division operator "\", which uses Cholesky factorization.In the latter case, which may be useful for larger systems, we use the preconditioned conjugate gradient solver with an incomplete Cholesky preconditioner (using MATLAB's pcg and ichol functions, respectively).For the iterative solver, the user must specify the convergence tolerance and the maximum number of iterations in the tol and maxit fields, respectively.
The code also has an option to use a GPU card to solve the system of linear equations, by setting the parameter FE.analysis.solver.use_gpu to true, which can be used if the system has an NVIDIA GPU card.In this case, only the iterative solver can be used, and a Jacobi preconditioner is used instead of the incomplete Cholesky preconditioner, since using the latter incurs a high data transfer cost to the GPU memory and is inefficient.The trade-off is that the iterative solution requires more iterations with the Jacobi preconditioner; nevertheless, for large meshes, it is still faster than using the iterative solver with CPUs.
Table 1 shows the time it takes to perform one optimization iteration using CPUs and GPUs in a system with an NVIDIA GTX 1070 Max-Q GPU card and an Intel core i7 8750H (with 6 cores) running on Ubuntu 19.10.The times correspond to the average iteration time of the first two iterations for the optimization of the default problem in the program (a 2D cantilever beam with eight bars).These times thus include not only the finite element analysis but also the geometry projection and the design update by the optimizer.Although the Jacobi preconditioner requires far more PCG iterations to converge, it still outperforms the use of CPUs for these mesh sizes (in this system).Also, it should be noted that the time and number of iterations will vary throughout the optimization, since different designs will lead to different condition numbers for the system matrix.As the table shows, it can be very convenient to use GPUs in this code.For instance, a mesh with half a million elements requires about one and a half minutes per iteration, therefore, 100 iterations of the optimization can be completed in approximately two and a half hours, which is-presently-very good for performing topology optimization in MATLAB on an individual workstation (recalling, however, that this is a 2D model).

Initial design
Through the GEOM.initial_design.pathfield in the master input file, the user must specify the path of the script containing the specification of the bars that make up the initial design.This script contains two arrays that describe the initial layout of bars.The point_matrix array has as many rows as points, with the first column containing an integer identifier for the point, and the remaining columns containing the spatial coordinates of the point.The bar_matrix array has as many rows as bars.Its first column corresponds to an integer identifier for the bar; the next two columns have the IDs of the points in the point_matrix array that correspond to the endpoints of the bar's medial axis; and the fourth and fifth columns contain the initial values of the bar's size variable and its radius, respectively.
It should be noted that this specification of points and bars allows for two possibilities: bars can be "floating" or "connected," as detailed in Section 2.2.The computation of sensitivities (cf.( 25)) in the code accounts for both situations.
As noted in Section 2.4, for sensitivities of the projected density of (2) to be well defined, it is necessary that r < r b .Since we typically consider a sample window that at least covers each finite element, as detailed in Section 3.9, this requirement imposes a minimum element size.For example, if r b min ¼ 1, and the sample window radius r is chosen to be twice the size of the element radius r e , then r e < 1/2 and therefore the element size h e should be at most ffiffi ffi 2 p =2 in 2D and ffiffi ffi 3 p =2 in 3D.To avoid undefined sensitivities in the initial design, it is also necessary that the length of each bar is not zero, i.e., x 1b ≠ x 2b .Therefore, the user should not create initial designs with zero-length bars.
When running the optimization, the code saves the current design (i.e., the arrays of points and bars) to a MATLAB .matfile.If the flag GEOM.initial_design.restart is set to true, the code reads the initial design from this file.This is a useful feature to start an optimization from the design obtained by another optimization run or to perform a finite difference check of the sensitivities (as detailed in Section 3.7) on the final design of a run.We note, however, that this does not constitute a clean continuation of a previous optimization run, since gradientbased optimizers (certainly the ones used in this code) use information from two or more consecutive iterations to construct approximations of the optimization functions.

Optimization
Although the code in its current form is written to solve the optimization problem of ( 15)-( 20), it is also structured so that (a) multiple constraints can be imposed and (b) any function can be chosen as the objective.This is achieved by indicating in the master input file the name of the objective function (e.g., OPT.functions.objective= 'compliance';) and the names and limits of the constraint functions (e.g., OPT.functions.constraints={'volumefraction', 'my constraint'}; and OPT.functions.constraint_limit= [0.31.0];).In this example, the user has to implement the function named 'my constraint' and add it to the list of available functions in the script init_optimization.m.This should facilitate the modification of our code to address other optimization problems.
As opposed to density-based topology optimization, the magnitudes of the design variables in the geometry projection scheme can differ greatly.Moreover, the optimization functions can be highly nonlinear with respect to the design variables, and therefore the optimizer can take too large a design step and produce poor iterates.To improve the performance of the optimization, we thus impose move limits on the design variables.In order to impose the same relative move limit on all variables, we first scale the design variables as where b z i is the scaled variable i and b z i and b z i are lower and upper bounds on the variables, respectively.For the components of the medial axis endpoints x 1b and x 2b , these bounds are given by the bounding box of the design region; for the bar radius, they correspond to the values specified by the user in GEOM.min_bar_radius and GEOM.max_bar_radius; and the size variables do not need any scaling.Variable scaling can be turned on or off using the parameter OPT.options.dv_scaling in the master input file.A move limit m is imposed on the design variables at iteration I as max 0;b z The move limit value is assigned with the parameters OPT.options.move_limit in the master input file.
Two optimizers can be used with our code: MATLAB's fmincon and the method of moving asymptotes (MMA) (Svanberg 1987(Svanberg , 2002(Svanberg , 2007)).The choice of optimizer is made with the parameter OPT.options.optimizer in the master input file.In the case of fmincon, we employ the activeset algorithm, since it allows us to impose move limits (through the 'RelLineSrchBnd' option as a relative bound on the line search step length), and it performs well in our numerical experiments.fmincon handles well the situation where the lower bounds on a variable equal the upper bounds, and therefore a fixed bar radius can be imposed if desired by setting GEOM.min_bar_radius = GEOM.max_bar_radius.
In addition to scaling the design variables, it is important to remember that in the case of MMA, it is recommended that the constraint limits and the objective function are between 1 and 100 for reasonable values of the design variables.The volume fraction automatically satisfies this requirement, but the compliance does not.If the magnitude of the applied loading is too large and/or the elasticity modulus of the material is relatively low, the compliance may easily exceed this range.The minimum compliance design for a given volume fraction does not depend on the magnitude of the load or the elastic modulus.Therefore, the user can modify, for example, the load magnitude to make sure the compliance is within the recommended range.Another, more general strategy is to simply scale the compliance and its sensitivities some factor, which requires modifying the compute_compliance.mscript accordingly.Ensuring proper scaling is important, as MMA may fail altogether if the magnitude of the compliance is very large.
In the case of MMA, we provide the code to call it in the script runmma.m,but the user must obtain the MATLAB version of MMA from its author and copy the files mmasub.m,subsolv.m,and kktcheck.m in the optimization folder of our code.We employed the 2007 version of MMA (Svanberg (2007)) for numerical tests with our code.Unlike fmincon, MMA does not handle well the situation where the lower and upper bounds of a variable are the same.Therefore, to approximately impose a fixed bar radius, the user should set GEOM.max_bar_radius=GEOM.min_bar_radius+delta, where delta is a small positive number.
The code employs three stopping criteria.The first criterion is satisfied if the 2-norm of the change in the vector of design variables falls below a specified value, indicated by the parameter OPT.options.step_tol.The second is satisfied if the norm of the Karush-Kuhn-Tucker optimality conditions falls b e l o w t h e v a l u e s p e c i f i e d i n t h e p a r a m e t e r OPT.options.kkt_tol.The third criterion is exceeding a maximum number of iterations, given by the parameter OPT.options.max_iter.If any of these criteria are satisfied, the optimization is stopped.

Output
The program produces several forms of output.During the course of the optimization, two MATLAB figures are created and updated at every iteration.The first is a plot of the bars drawn by using directly the geometric parameters as rectangles with semicircular ends and colored such that the transparency indicates the value of the size variable.A bar with a size variable value of unity has no transparency, and bars with a size variable value less than 0.05 are removed from the plot.
The second figure is a plot of the combined density of (5).In the case that fmincon is used as optimizer, this figure has two buttons to either pause or stop the optimization.At the end of the optimization, a third figure is produced that plots the iteration history of the objective and constraint functions.The generation of these plots can be turned off by setting the variable plot_cond to false in the master input file; this option is convenient when, for example, the code is executed in a queue in a high-performance computing system.The code also writes information on the optimization iterations to MATLAB's Command Window.
Finally, our code writes Visualization Toolkit (.vtk) files of the combined density of the design to the folder specified in the OPT.options.vtk_output_pathparameter (which by default is set to 'output_files').The user can request VTK output for none, all, or the last iteration by setting the parameter OPT.options.write_to_vtkto 'none', 'all', or 'last', respectively.These files can be visualized with Open Source tools like ParaView (Ahrens et al. 2005;Ayachit 2015), which offer wide post-processing capabilities.This is particularly useful for 3D problems.

Finite difference check
Another utility provided by our code is a forward finite difference check of the analytical sensitivities of the cost and/ or objective functions.The finite difference check is turned on by setting the parameter OPT.make_fd_check to 'true' in the master input file.If this option is chosen, the code performs the finite difference check and subsequently exits.The size of the finite difference step is specified in the parameter OPT.fd_step_size.The user can choose what function should be checked by setting the parameters check_cost_sens and check_cons_sens to 'true' or 'false'.This functionality is useful when developing new optimization functions to ensure the accuracy of the analytical sensitivities.

Distance calculation
The computation of the distance d be of ( 14) from the centroid x e of each element e to the medial axis of each bar b requires, conceptually, a double for loop, as shown in Fig. 3.This computation, which is of order O(n e n b ), can be quite expensive; however, it fortunately is embarrassingly parallel.In distributed memory implementations, the elements in the mesh can be divided among the available compute cores and each core calculates the distance from those elements to all of the bars (see, e.g., Zhang et al. (2016a)).This computation scales linearly with number of cores, and therefore if enough cores are employed for the distance computation (as well as the calculation of the combined density discussed in the following section), eventually the finite element analysis dominates the cost of each function evaluation as the number of elements increases, and the geometry projection only represents a small fraction of the cost.
When running on a single workstation with multiple cores, MATLAB employs shared-memory processing for many of its vectorized operations.Therefore, as discussed in Andreassen et al. (2011), one of the keys to efficient implementations of numerical methods in MATLAB is to vectorize the computations as much as possible.In particular, whenever a loop can be replaced with a vectorized operation, the improvements in performance can be drastic.
In the case of the distance calculation, the strategy we used to vectorize the double for loop is to employ threedimensional arrays, as can be seen in the excerpt of the distance computation script shown in Appendix A. For instance, the array x_e_1b (line 36 in the code shown in the Appendix) contains all the vectors x e/1b , used in (12)(13)(14).The first dimension of the array corresponds to the spatial components in the vector (of size 2 and 3 in 2D and problems, respectively); the second dimension corresponds to the bar; and the last to the element.Note that even the computation of the branched function of 2 is vectorized by using Boolean matrices of dimensions 1 × n b × n e whose components equal true if the condition for the corresponding branch is satisfied and false otherwise (which are equivalent to 1 and 0, respectively).
As the problem size increases, accessing three-dimensional matrices can get quite slow (particularly if they have to be stored out of memory); therefore, eventually this approach becomes inefficient.However, in our numerical experiments (and using our hardware), we have been able to run problems with meshes that have hundreds of thousands of elements and tens of bars in reasonable times, and therefore this code should in many cases be efficient enough for method development (indeed, most publications with methods to perform topology optimization with geometric components use mesh sizes well within this range).
Another important aspect of the code that is not explicitly stated in the pseudo-code of Fig. 3 is that the sensitivities of each quantity are computed in the same script where the quantity is computed, and they are then stored in the global structures or passed to the calling function so that the chain rule is used for subsequent sensitivity calculations.For instance, the derivatives of the distance d be in the script of Fig. 19 with respect to the design parameters are computed in the same script.This makes for a more modular code and makes it easier to verify the correctness of the individual terms of the sensitivities.As with the distance calculation, the computation of the sensitivities makes extensive use of vectorization for efficiency.

Combined density calculation
The computation of the combined density of ( 5) is also vectorized for efficiency.By default, the radius r e of the sample window in (2) is computed as the radius of the circle (in 2D) or sphere (in 3D) that circumscribes the square or cubic element, respectively, with the same volume of the element: where d = 2, 3 for 2D and 3D, respectively.This default can be overridden by uncommenting the line with the parameter OPT.parameters.elem_r in the master input file and assigning to this parameter the actual value of the sample window radius.In this case, the same radius will be used for all elements in the mesh.
The code has two options for the penalization function μ of (4), implemented in the script penalize.m:the SIMP penalization mentioned in Section 2.1 and the rational approximation of material properties (RAMP) (Stolpe and Svanberg 2001).The choice of penalization scheme is indicated by the parameter OPT.parameters.penalization_scheme in the master input file, and the penalization parameter is assigned in OPT.parameters.penalization_param.
Several combination strategies are also available to compute the smooth maximum g max of (5), implemented in the script smooth_max.m.These include the modified p-norm of  m a s t e r i n p u t f i l e ) .S e c o n d , t h e p a r a m e t e r s OPT.make_fd_check and OPT.check_cost_sens must be set to true.The perturbation size for the finite difference check is set as OPT.fd_step_size=1e-6;.
The largest absolute and relative differences between the analytical sensitivities and the finite difference sensitivities are reported to MATLAB's Command Window as − 0.0037 and − 0.0013, respectively, indicating a good agreement.They correspond to the x component of the eighth point, which is the point closest to the right-hand side corner of the top edge.The code also produces a plot comparing the two sensitivities, shown in Fig. 14.

3D cantilever beam
The last example corresponds to a 3D cantilever beam with fixed supports on all four corners at one end, and a tip downward load in the center of the opposite face, as shown in The mesh for this problem is generated automatically, and it consists of 80 × 40 × 40 = 128000 elements.This problem is solved using GPUs on a machine with 24 Intel Xeon CPUs at 2.2 GHz, 32 Gb of RAM, and an NVIDIA Quadro M2000 GPU card; running on Ubuntu 18.04.3LTS; and using MATLAB R2019b.
The problem was solved in 63 min, and it took 106 iterations to convergence.The optimal design for this problem is  A MATLAB code for topology optimization using the geometry projection method

Fig. 2
Fig. 2 Definition of vectors and quantities for distance computation

Fig. 5
Fig. 5 Initial design for MBB beam

Fig. 6
Fig. 6 Optimal design for MBB beam

Fig. 15 .
The magnitude of the force is F = 0.1.We minimize the compliance subject to a volume fraction constraint of v * f ¼ 0:1.The initial design for this example, shown in Fig.16, is made of 16 floating bars with 32 points.Bounds on the bars' radii of 0.5 ≤ r b ≤ 1.0 are imposed and MMA is used as an optimizer.

Fig. 14
Fig. 13 Optimization history of compliance and volume fraction for L-bracket problem

Fig. 19
Fig.19Excerpt of distance computation script (the portion with sensitivities calculation is omitted for brevity) Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.