Density and Level SetXFEM Schemes for Topology Optimization of 3D Structures
Abstract
As the capabilities of additive manufacturing techniques increase, topology optimization provides a promising approach to design geometrically sophisticated structures which can be directly manufactured. Traditional topology optimization methods aim at finding the conceptual design but often lack a sufficient resolution of the geometry and structural response, needed to directly use the optimized design for manufacturing. To overcome these limitations, this paper studies the viability and characteristics of the eXtended Finite Element Method (XFEM) in combination with the LevelSet Method (LSM) for topology optimization of three dimensional structural design problems. The LSM describes the geometry by defining the nodal level set values via explicit functions of the optimization variables. The structural response is predicted by a generalized version of the XFEM. The LSMXFEM approach is compared against results from a traditional Solid Isotropic Material with Penalization (SIMP) method for twophase “solidvoid” and “solidsolid” problems. The numerical results demonstrate that the LSMXFEM approach can describe crisply the geometry and predict the structural response of complex threedimensional structures with acceptable accuracy even on coarse meshes. However, the LSMXFEM studied here lacks a robust and intuitive formulation to control the minimum feature size, and the optimization results may depend on the initial design.
Keywords:
eXtended Finite Element Method Topology Optimization Solid Isotropic Material with Penalization Level Set Methods Additive Manufacturing 3D Printing∎
1 Introduction
Recent advances in additive manufacturing allow the precise placement of one or multiple materials at micrometer resolution with essentially no restrictions on the geometric complexity of the spatial arrangement. Complex three dimensional solids can be created with highly nonregular material distributions in a near optimal fashion, enabling the fabrication of structures with enhanced performance. Topology optimization has emerged as a promising approach to utilize the benefits of additive manufacturing (Ning and Pellegrino, 2012; Meisel et al, 2013). Structural topology optimization seeks to find the optimal geometry and/or the material layout of a body within a given design domain. The geometry is represented by the spatial distribution of two or more material phases; in structural problems, one of these material phases may represent void.
Originally topology optimization methods were developed primarily to create conceptual designs in the early stage of the design process (Bendsøe and Sigmund, 2003; Rozvany, 2009). Later, topology optimization was applied to directly design microelectromechanical systems (MEMS), utilizing the ability of thinfilm fabrication techniques, such as photolithography in combination with chemical etching, to create geometrically complex devices at low cost (Sigmund, 2001a, b). As the MEMS design problem is essentially twodimensional, traditional approaches were sufficient to achieve the necessary geometric resolution at acceptable computational cost. Motivated by the availability of affordable additive fabrication methods for three dimensional structures, this paper focuses on topology optimization of three dimensional structures and introduces a new approach for finding optimized designs with high geometric resolution on rather coarse computational meshes.
Most approaches for structural topology optimization are density methods. For a twophase problem, the density is considered a design variable and can assume intermediate values between the density of the material phase “A” and the density of the material phase “B”. The most popular density method is the Solid Isotropic Material with Penalization method introduced by Bendsøe (1989) and Zhou and Rozvany (1991). It features great versatility, robustness, efficiency, and ease of use and implementation for a broad range of applications (Sigmund and Maute, 2013; Deaton and Grandhi, 2013).
Density methods typically describe the boundaries between the material phases either via intermediate density values or by discrete material distributions leading to jagged boundary geometries. In both cases, the enforcement of boundary and interface conditions at the material interface is hampered and may result in nonphysical responses, such as premature yielding (Maute et al, 1998). Often this issue can be mitigated by mesh refinement or adaptive remeshing (Maute and Ramm, 1995, 1997). However, for problems that require an accurate description of the boundaries, such as boundary layer problems in fluids and skindepth issues in electromagnetics, it was reported that density methods like SIMP fail (Sigmund and Maute, 2013).
The shortcomings of density methods have promoted the development of the Level Set Method (LSM) for topology optimization. The LSM allows a crisp representation of the phase boundaries and the accurate enforcement of boundary conditions on fixed meshes. The material interface in the LSM is described implicitly by the isocontours of a Level Set Function (LSF), usually the zero levelset contour (Allaire et al, 2004; Sethian and Wiegmann, 2000; Wang et al, 2003).
The LSF is typically discretized by the same mesh used for the physical field and is updated in the optimization process via the solution of the HamiltonJacobi equations. Alternatively, the parameters of the discretized LSF are defined as explicit functions of the optimization variables, and the resulting parameter optimization problem is solved by standard nonlinear programming methods (van Dijk et al, 2013). The key challenges for the LSM include (a) controlling the spatial gradients of the LSF in the vicinity of the zerolevel set contour to avoid illconditioning of the optimization problem, (b) controlling local feature sizes, and (c) accelerating the convergence of the geometry in the optimization process. For a detailed discussion of the LSM, the reader is referred to the comprehensive review by van Dijk et al (2013) and Gain and Paulino (2013).
In LSMs, the structural geometry can be represented in the mechanical model via an Ersatz material approach, immersed boundary techniques, or by adaptive geometry conforming meshes. The first two approaches allow the use of fixed, design independent meshes while the last approach requires local or global remeshing as the structural geometry evolves in the design process.
Using an Ersatz material approach, the void phase is modeled by a soft material and the material properties in elements intersected by the zerolevel set contour are interpolated between the “void” and solid phase, proportional to the volume ratio of the individual phases. However, this approach faces the same issues as density methods in regards to enforcing boundary conditions across the material interface. Other approaches to model the mechanical response include generalized and adaptive finite element schemes such as the SuperImposed Finite Element Method (SFEM) (Wang and Wang, 2006), the eXtended finite element method (XFEM) (van Miegroet and Duysinx, 2007; Wei et al, 2010; Kreissl and Maute, 2012), and local remeshing schemes (Yamasaki et al, 2011).
In this paper, we focus on the LSM in combination with the XFEM. The XFEM does not require a mesh that conforms to the material interfaces and reduces the complexity of mesh construction. Spatial discontinuities in the structural response are captured by augmenting the standard finite element interpolations with additional shape functions. This approach is similar to the SFEM (Wang and Wang, 2006) as the solution is obtained by superimposing the standard and enriched shape functions. However, unlike the SFEM, the XFEM can combine multiple types of shape functions and thus allows for greater flexibility.
The XFEM builds upon the partition of unity concept developed by Babuška and Melenk (1997). The idea of XFEM was originally proposed by Belytschko and Black (1999) to model crack propagation. The reader is referred to Yazid et al (2009) for an overview of the application of XFEM to problems in fracture mechanics. The XFEM has been used for a variety of problems in computational mechanics, such as fluidstructure interaction (Gerstenberger and Wall, 2008b, a), multiphase flows (Fries, 2009), and nanoscale heat transfer (Lee et al, 2011). A general overview of the method is presented by Fries and Belytschko (2010).
Duysinx et al (2006) originally introduced the XFEM into shape optimization using a simplified XFEM formulation. This formulation does not use additional enrichment functions and is limited to problems where one of the material phases represents void and the geometric configuration is “simple”, i.e. does not contain geometric features that are smaller than the size of two elements (Makhija and Maute, 2013). In this instance, the weak form of the governing equations is only integrated over the solid material in each element. In addition, if the interface between the material phases is traction free, this simplified version of the XFEM only differs from the traditional finite element method with respect to the domain of integration. The simplified XFEM version was applied to structural shape optimization, for example, by Duysinx et al (2006), van Miegroet et al (2005), and van Miegroet and Duysinx (2007), and to topology optimization of three dimensional structures by Li et al (2012).
An XFEM approach based on a standard enrichment strategy allows to model twophase problems with a simple intersection pattern. This formulation is considered, for example, by Wei et al (2010) to solve “solidvoid” structural topology optimization problems, modeling the “void” as a soft material. Maute et al (2011) used a standard enrichment strategy to discretize the phonon Boltzmann transport equation and optimize the thermal conductivity of nanocomposites. However, this enrichment strategy is not guaranteed to consistently approximate the state variable fields for configurations with complex intersection patterns.
An enhanced version of the XFEM was proposed by Makhija and Maute (2013), who presented a generalized enrichment based on the step enrichment of Hansbo and Hansbo (2004) and applied it to two dimensional structural topology optimization. This formulation captures consistently the mechanical response for complex geometries and intersection patterns for general multiphase problems.
This paper will expand the combination of the LSM and the XFEM onto general twophase, three dimensional problems. We will compare results of the proposed LSMXFEM framework with SIMP results for structural topology optimization examples. The numerical results will show that the LSMXFEM combination is a promising approach for three dimensional problems and allows for the use of coarse meshes to represent the structural geometry and to describe the structural response with acceptable accuracy.
The main challenges of expanding the previous LSM with XFEM approaches to three dimensions stem from the increased complexity of possible intersections patterns. Such patterns include elements that are intersected multiple times and elements containing only a small volume of a particular phase. In contrast to Li et al (2012), we will adopt the generalized enrichment strategy of Makhija and Maute (2013) to (a) accurately model the structural response on complex three dimensional patterns and (b) solve solidsolid material distribution problems. Further, we will expand the preconditioning scheme of Lang et al (2013) onto three dimensions to mitigate illconditioning issues in the XFEM analysis problems due to elements with small volume fractions.
The main contributions of the paper are: (a) we present a numerically robust and computationally viable approach for solving general twophase, three dimensional topology optimization problems, and (b) we provide a direct comparison of LSMXFEM and SIMP results for three dimensional problems, highlighting key features of the two methods.
The remainder of this paper is structured as follows: the geometry models of the LSMXFEM and SIMP methods are described in Section 2. The mechanical model and the XFEM formulation are summarized in Section 3. Details of the LSMXFEM and SIMP optimization approaches are presented in Section 4. Section 5 highlights specific computational challenges of the XFEM approach for three dimensional problems. Section 6 presents structural topology optimization examples, comparing the LSMXFEM and SIMP approaches. Section 7 summarizes the main conclusions drawn from this study.
2 Geometry Modeling
In topology optimization, the geometry of a body is defined via its material distribution. In density methods, such as SIMP, the material distribution is discretized by finite elements, with either elemental or nodal parameters defining the distribution within the element. Most often the same mesh is used to approximate the density distribution and the structural response. Alternatively, the state and density fields can be discretized by different meshes with different refinement levels; see for example the Multiresolution Topology Method (MTOP) by Nguyen et al (2010). The optimization variables define analytically or by means of auxiliary partial differential equations the nodal or element density parameters (Sigmund and Maute, 2013). For twophase problems, the density is continuously varied between “0” (phase “A”) and “1” (phase “B”). Implicit or explicit penalization schemes, optionally combined with projection methods, are used to encourage “01” solutions (Guest et al, 2004; Sigmund, 2007).
The crispness of the interface geometry, as described via the optimized material distribution, depends on (a) the formulation of the optimization problem, i.e. the objective and constraints, (b) regularization techniques, such as density or sensitivity filters, and (c) the optimization algorithm. In general, the resolution of the phase boundaries increases as the mesh is refined. For three dimensional problems, often coarse meshes are used to limit the computational costs. In this case, the boundary geometry either lacks crispness due to the presence of elements with intermediate densities or is approximated by a spatially discontinuous material distribution, leading to jagged interfaces.
Alternatively, the material distribution can be described via a level set function . Typically the zero level set contour defines the phase boundaries. Considering a twophase problem, the interface is defined as follows:
(1)  
where the vector collects the spatial coordinates, is the domain of material phase “A”, is the domain of material phase “B”, and defines the material interface between phase “A” and phase “B”. For example, to model a sphere in a three dimensional mesh centered at (,,), the value of the level set function at a grid point (,,) is:
(2) 
where represents the radius of the sphere, and the sign value of at each node determines if the node is inside or outside the circle.
The level set field is typically discretized by shape functions with either local or global support; the reader is referred to the review paper by van Dijk et al (2013). In the simplest and most common approach, the levelset field is approximated on the same mesh used for discretizing the governing equations. In this study, we follow this approach and define the nodal levelset values explicitly as functions of the optimization variables (see 4.1). The resulting parameter optimization problem is solved by a standard nonlinear programming method.
While LSMs provide a crisp definition of the phase boundaries, they require seeding the initial design with inclusions and/or introducing inclusions in the course of the optimization process, for example via topological derivatives (Eschenauer et al, 1994; Sokolowski and Zochowski, 1999; Burger et al, 2004; Norato et al, 2007). An example of an initial design with a regular pattern of spherical inclusions is shown in Fig. 1. The images in the upper row and the image in the lower left corner show only phase “A”. The material layout of both phases is depicted in the lower right image. This layout can be generated by superposing Eq. 2 for spheres at uniformly spaced center locations.
The optimization results of the LSM are typically dependent on the initial layout. Furthermore, it is nontrivial to generate an initial design that satisfies geometric design constraints, such as volume or perimeter constraints. In our experience, severe constraint violations may cause the optimization process to converge to a feasible design with otherwise poor performance.
3 Structural Analysis
In this section, we briefly discuss the structural model and the XFEM analysis used in this paper. The governing equations are presented first, followed by a summary of the XFEM discretization and analysis.
3.1 Structural Model
This paper considers the topology optimization of structures using the LSM approach described above and the XFEM to predict the structural response, assuming infinitesimal strains, a linear elastic material behavior, and static conditions. We consider the twophase problem depicted in Fig. 2, where denotes the surface where traction forces are applied, denotes the surface with prescribed displacements, and is the normal at the material interface pointing from phase “A” to phase “B”. The weak form of the governing equations can be decomposed into the following terms:
(3) 
where collects the contributions from the static equilibrium, including body forces and surface tractions, models the interface conditions along the phase boundaries for “solidsolid” problems, and is due to a fictitious spring model to pin free floating material in “solidvoid problems”.
The weak form of the structural equilibrium equations is:
(4) 
where is the kinematically admissible test function, is the strain tensor associated with the test function , is the displacement vector, is the stress tensor, is the applied body force, and is the external traction applied along .
To enforce continuity of the displacements along the phase boundaries, the static equilibrium equations are typically augmented by either an enhanced Lagrange multiplier or penalty formulations, such as the stabilized Lagrange multiplier and the Nitsche method. Note that the standard Lagrange multiplier approach is not suitable for the XFEM as it suffers from numerical instabilities. The reader is referred to Stenberg (1995), Juntunen and Stenberg (2009), and Dolbow and Harari (2009) for more details.
In this paper, we enforce displacement continuity along phase boundaries for “solidsolid” problems via the following stabilized Lagrange multiplier method:
(5) 
(6) 
(7) 
where is the Lagrange multiplier, and is the associated test function. The higher the weight is, the better the interface condition is satisfied, at the cost of numerical stability.
The formation of free floating solid particles surrounded by void material is possible in “solidvoid” topology optimization, leading to a singular analysis problem. This issue does not exist in an Ersatz material approach as the void phase is modeled via a soft material. A similar approach can be applied to the XFEM to suppress singularities by modeling the “void” phase via a soft material (Wei et al, 2010). However, this approach requires accounting for the interface contributions (5) and integrating the governing equations over the void phase. To avoid the associated complexity and computational costs, we extend the approach of Makhija and Maute (2013) onto three dimensions and assume that the solid phase is supported by fictitious springs. This model leads to the following contribution to the governing equations, assuming that phase “A” is the solid phase:
(8) 
where denotes the stiffness of the distributed system of springs.
For a more detailed explanation of this XFEM formulation, the reader is referred to the paper by Makhija and Maute (2013).
3.2 Discretization
To capture the discontinuities in the strain and stress fields along the phase boundaries, we enrich the standard finite element approximation with additional shape functions. We adopt the generalized enrichment strategy of Makhija and Maute (2013) which resolves consistently the displacement fields in the presence of small features and does not suffer from artificially coupling disconnected phases. Considering a twophase problem, the displacement field is approximated as follows:
(9) 
where is the enrichment level, is the maximum number of enrichment levels used for each phase, are the shape functions, is the vector of nodal displacement components at node for phase , is the level set value evaluated at the integration point, and denotes the Heaviside function. The enrichment level is chosen such that the displacements in disconnected volumes of the same phase are interpolated by separate sets of degrees of freedom. When interpolating the level set field by elementwise linear functions, a maximum of enrichment levels is needed in three dimensions. This enrichment strategy will be revisited in Section 5.
The Heaviside function depends on the level set function and is defined as follows:
(10) 
The Heaviside functions “turns on/off” the standard finite element interpolations in the particular phases. The approximation (9) allows for discontinuities of the displacements along the phase boundaries. Therefore the continuity is enforced weakly via the stabilized Lagrange multiplier method (5).
Following a BubnovGalerkin scheme, we test the governing equations with the same subspace as we use for the trail functions; see Eq. 9. The weak form of the governing equations is integrated numerically over the individual phases, using the Delaunay triangulation of the element along the phase boundaries.
3.3 Preconditioner
As described above, the degrees of freedom interpolate the structural displacements in topologically connected subdomains of phase in the elements connected to node . As the total volume of these subdomains vanishes, the discretized structural model becomes increasingly illconditioned; i.e. the condition number of the stiffness matrix rapidly increases. This phenomenon is more pronounced for three dimensional problems than those in two dimensional ones.
To mitigate this illconditioning issue, we extend the geometric preconditioning scheme of Lang et al (2013), which was introduced and studied for two dimensional heat conduction and flow problems, onto three dimensional problems in structural mechanics. The goal of this preconditioning scheme is to balance the influence of all degrees of freedom in the system, as the volumes in which the subset of these degrees of freedoms interpolates the solution approach zero. To this end, we introduce the following projection:
(11) 
where is the vector of displacement degrees of freedom according to Eq. 9, is a transformation matrix, and is the solution vector in the transformed space. The residual, , and stiffness matrix, , in the transformed space are defined as:
(12) 
where the residual, , and the stiffness matrix, , result from integrating the weak form of the governing equations using the XFEM approximation (9).
The preconditioner is a diagonal matrix built by integrating the spatial derivatives of the shape functions over the nodal support of nodes connected to an intersected element. The diagonal components of the matrix are defined as
(13) 
where corresponds to the degree of freedom , is the node index, is the material phase “A” or “B”, is the enrichment level, is the set of elements connected to node , and is the element domain of phase . The components of the matrix increase as the region of influence of a degree of freedom decreases. The entries of nodes that are not connected to at least one intersected element are set to one.
To avoid numerical issues due to large values for the components of , the degrees of freedom associated with the diagonal entry are constrained to zero if the following condition is satisfied:
(14) 
where is a specified tolerance. As studies by Lang et al (2013) have shown, the above preconditioning scheme is rather insensitive to the value of and is typically set to a value larger than . For more details about the specifics of the formulation, the reader is referred to the paper by Lang et al (2013).
4 Optimization Model
The design optimization problems considered in this paper can be written as follows:
(15)  
where denotes the vector of design variables, the objective function, and the th design constraint. In general, the objective and constraints depend on the optimization and state variables. The optimization problem (15) is solved by nonlinear programming methods, and the gradients of the objective and constraint functions are computed via the adjoint method.
In this paper, we compare the proposed LSMXFEM against the wellknown SIMP method, augmented by a projection scheme. In the following subsection we briefly outline the models that define the discretized level set field (LSM) and the material properties (SIMP) as function of the optimization variables.
4.1 Xfem
The nodal values of the discretized level set field are defined as analytical functions of the optimization variables via the following linear filter:
(16) 
with
(17) 
where is the number of nodes in the discrete model, is the location of the node at which the design variable is defined, is the filter radius, is the weight of node with respect to design variable , and and are the position vector and the computed level set value of node , respectively.
The above linear filter was used previously in the studies of Kreissl and Maute (2012) and Makhija and Maute (2013), and was shown to improve the convergence rate in the optimization process. However, in contrast to density or sensitivity filters used in SIMP methods, the filter above is not guaranteed to control the minimum feature size. This issue will be revisited in Section 6.
4.2 Simp
Here, the material distribution is parameterized by nodal density values, , which are treated as optimization variables, i.e. . Following the work of Guest et al (2004), we compute the elemental density by combining a linear density filter and a projection scheme as follows:
(18) 
where
(19) 
where is the number of elements in the discrete model, is the location of the node at which the design variable is defined, is the filter radius, is the factor of element with respect to design variable , and and are the position vector of the centroid and the computed elemental density of element , respectively.
(Guest et al, 2004) proposed a density projection method to reduce the volume occupied by material with intermediate densities. The projection is based on a smoothed Heaviside function and applied to the elemental densities as follows:
(20) 
where is the projected elemental density, and the parameter controls the crispness of the projection. For the projection turns into an identity operator, i.e. .
The Young’s modulus, , is defined as a function of the density, , using the standard SIMP interpolation:
(21) 
where and are the Young’s moduli for material phase “A” and “B”, and is the SIMP penalization factor. To model a “solidvoid” optimization problem, is set to value much smaller than .
The filter (18) prevents the formation of features smaller than , typically at the cost of generating intermediate density values along the phase boundaries. This effect is mitigated by the projection (20) which, in the limit for , maps nonzero values into “1”. The reader is referred to the papers by Guest et al (2004) and Guest et al (2011) for further details of the scheme presented above, and to Sigmund and Maute (2013) for a comprehensive discussion of projections schemes.
5 Computational Considerations
Expanding the LSMXFEM combination onto three dimensional problems faces both algorithmic and computational challenges which are briefly discussed below.
The XFEM requires integrating the weak form of the governing equations separately in each phase. To this end, an element intersected by the zero level set contour is subdivided. For two dimensional problems and using a linear interpolation of the level set field within an element, there are only 8 intersection configurations which can be tabulated; see Fig. 3. In three dimensions, there are 127 intersection configurations. To handle this complexity, we compute the intersection point of the zero level set contour with the element edges and use a Delaunay triangulation to subdivide the element. In numerical experiments, this approach has proven robust and computationally inexpensive.
Previous studies on topology optimization for three dimensional structures with the XFEM (Li et al, 2012; Herrero et al, 2013) have employed a simplified enrichment scheme which is limited to “solidvoid” problems and may suffer from artificial coupling of disconnected material. Our work overcomes these issues by adopting the generalized enrichment scheme summarized in Section 3.2. The key challenge of this scheme is to identify the enrichment levels needed to consistently interpolate the displacements in elemental subdomains with the same phase.
To this end, the subdomains in all elements connected to a node need to be considered. This naturally leads to an algorithm which loops over all nodes and, in an inner loop, over all elements connected to the current node. As this approach processes an element repeatedly, the following simple and efficient twostep scheme is introduced:

A temporary, elemental enrichment level is assigned to the subdomains in each element. Recall that the enrichment level defines the set of degrees of freedom used to interpolate the displacements in an elemental subdomain. Note that because this assignment is done individually for each element, the continuity of the interpolation across elements is not guaranteed. Fig. 4 shows the triangulation and enrichment level for the red phase in two and three dimensions.

The nodal enrichment levels are constructed to ensure that the displacement field is interpolated continuously across elements, and by a different set of shape functions for each disconnected elemental subdomain of the same phase. To this end, the cluster of elements connected to a node is considered, and the elemental enrichment levels assigned in step 1 are adjusted to satisfy the continuity and consistency conditions.
This process is illustrated in Fig. 5. The node of interest is the one located in the center of the element cluster. In step 1 each subdomain of the red phase is assigned an enrichment level of . Applying this enrichment level to the degrees of freedom for the red phase at the center node would incorrectly couple the displacement fields in the red phase subdomains. Analyzing the element cluster around the center nodes shows that these subdomains are disconnected and individual enrichment level are assigned.
Topology optimization in three dimensions leads to FEM or XFEM models with a large number of degrees of freedom, and typically requires using iterative solvers and parallel computing. The stabilized Lagrange multiplier formulation of the interface conditions (5) leads to a nonsymmetric stiffness matrix. Numerical experiments have shown that the XFEM problems considered in this study can be robustly and efficiently solved by a generalized minimal residual (GMRES) method preconditioned by incomplete LU (ILU) factorization. Note that the ILU preconditioner operates on the projected XFEM system (12).
6 Numerical Examples
We study the features of the proposed LSMXFEM topology optimization approach with numerical examples. The LSMXFEM results of “solidvoid” and “solidsolid” problems are compared against the ones of the SIMP approach outlined in Section 4.2. In all examples we seek to minimize the strain energy subject to a constraint on the volume of the stiff phase. This problem formulation is chosen because it is well studied in the literature and the numerical experiments can be easily repeated. The following numerical studies will provide insight into (a) the convergence of the geometry and the structural response as the meshes are refined and (b) the influence of regularization techniques on the optimized results, such as the filter radii in (16) and (18), and perimeter constraints.
In all examples, the optimization problems are solved by the Globally Convergent Method of Moving Asymptotes (GCMMA) of Svanberg (2002). The sensitivities are computed by the adjoint method. The design domains are discretized by 8node linear elements. The linear systems of the forward and adjoint problems are solved by a parallel implementation of the GMRES method (Heroux et al, 2003). The problems are preconditioned by an ILU factorization with a fill of and an overlap of . The convergence tolerances for both, the GCMMA and the GMRES solver, are chosen sufficiently low such that the optimization results do not depend on the tolerance values. In the LSMXFEM examples, the spring stiffness value, , is .
While the LSFXFEM results can be directly used to fabricate the structure, for example by 3D printing, the SIMP results need to be postprocessed. From a practitioner perspective, only the postprocessed SIMP results should be compared against the LSFXFEM results. To this end we postprocess the SIMP results with a lumping method that uses the isocontour of the density distribution. To obtain a “01” density distribution with smooth phase boundaries, we construct isocontours for from the nodal density values, ; see Section 4.2. The volume enclosed by the isocontour with is considered solid; the remaining volume is considered “void”. The threshold, , is determined such that the volume of the solid domain satisfies the volume constraint. The structural response of the results postprocessed with the isocontour density lumping (IDL) approach is analyzed conveniently with the XFEM.
To gain further insight into the crispness of the SIMP results and the influence of the postprocessing methods above on their performance, we measure the volume fraction, , occupied by elemental densities with as follows:
(22) 
where denotes the design domain.
6.1 Cube with center load
We consider the “solidvoid” optimization problem depicted in Fig. 6. With this example we will illustrate the basic features of the LSMXFEM approach for three dimensional problems and show that the proposed LSMXFEM approach and the SIMP formulation may exhibit comparable convergence behaviors as the mesh is refined.
The cubical design domain is pinned at its four bottom corners in the vertical direction and a unit force is applied at the center of the bottom face. The Young’s modulus of the stiff phase is set to and the Poisson ratio to . The maximum volume of the stiff phase is . We compare LSMXFEM and SIMP results for two mesh sizes: and . The problem is solved on the full mesh.
First, we apply the SIMP approach with a penalization factor of . The size of the smoothing radius is mesh dependent, and is set to for the coarse mesh and for the fine mesh; the projection parameter is set to . Note that the smoothing radius is intentionally set relative to the element size ( the element edge). The Young’s modulus of the void phase is set to . While this approach does not ensure meshindependent optimization results, it still prevents the formation of checkerboard patterns and provides insight into the dependency of the geometry resolution of SIMP as the mesh is refined.
The design domain is initialized with a uniform material distribution of . The optimized material distributions are shown in Figure 7 where material with a density lower than is considered void. The strain energies are reported in Tab. 1. For both meshes the volume constraint is active in the converged designs. As expected, the optimized geometry is smoother and the strain energy is lower for the refined mesh.
The SIMP results for the coarse and fine mesh are postprocessed with the IDL approach described above. The strain energies for varying threshold values, , are plotted in Fig. 8. The volume constraint is met for for the coarse mesh and for the fine mesh. The difference between the SIMPIDL and the SIMP results is of for the coarse mesh and for the fine mesh. The value of is higher for the coarse mesh because it cannot converge to a design with void inclusions. In both cases the strain energies of the postprocessed results match well the SIMP predictions as the density distributions converged well to “01” solutions.
The volume fractions of intermediate densities are and for the coarse and fine mesh, respectively. The postprocessed designs have lower strain energies because the postprocessing counteracts the effect of the density filter (18).
Mesh size  Strain energy  

SIMP  24x24x24  9.1456e01 
65x65x65  3.5244e02  
XFEM  24x24x24  1.0082e+00 
65x65x65  3.5519e02 
The same optimization problem is solved with the proposed LSMXFEM approach. The smoothing radius is set to for the coarse mesh and for the fine mesh. No perimeter constraint is imposed. We seed the initial design with two different configurations of void inclusions to study the influence of the initial layout on the optimization results. For both configurations we impose an equally spaced array of squareshaped holes with rounded corners, by modifying Eq. 2 into the following:
(23) 
One configuration has equally spaced holes with radius , the other holes with radius , as shown in Fig. 9. In both cases, the volume constraint is not satisfied with the initial design. Note that no inclusions are placed at the four bottom corners where the boundary conditions are applied.
Both level set configurations converge to nearly indistinguishable designs and strain energy values, for both the coarse and fine meshes. The optimized designs are shown in Fig. 10. The strain energies of the optimized designs are given in Tab. 1. The convergence history for the coarse meshes in SIMP and LSMXFEM is shown in Fig. 11.
For the example consider here, the SIMP and LSMXFEM results match well, both in regards to the geometry and the strain energy values. The LSMXFEM approach shows a faster convergence as the mesh is refined, for both the coarse and fine meshes. Comparing the optimized geometries, the SIMP results contain more structural features for both mesh resolutions. For example, considering the fine mesh, SIMP generates two small holes in the webs connecting the supports to the load point, while XFEM leads to only one larger hole, independent of the initial design configuration. However, these small differences have only a minor impact on the structural performance, i.e. the strain energy, of the optimized designs.
Considering the conceptual structural layout, both, the SIMP and the LSMXFEM approach, display only minor mesh dependencies for the problem studied here. Although the strain values show significant differences, the optimized geometries obtained with the coarse and fine meshes differ insignificantly for the SIMP and LSMXFEM approach. The following example will demonstrate a less benign convergence and identify more pronounced differences between SIMP and LSMXFEM.
6.2 Cuboid under torsion
The second “solidvoid” example is taken from Nguyen et al (2012) and reveals differences in the SIMP and LSMXFEM approaches. We will show that, without imposing a meshindependent minimum feature size constraint, the proposed LSMXFEM approach may feature a significantly greater convergence than the SIMP method employed in this paper. However, we will also illustrate that our LSMXFEM approach suffers from a lack of a robust and intuitive shape control technique.
The design domain is a cuboid of size , as shown in Fig. 12. A torque moment is generated via unit loads acting at the centers of the edges of the top face. The design domain is clamped at the bottom face. The Young’s modulus is set to and the Poisson ratio to . The volume of the stiff phase is constrained to % of the total volume. The problem is solved on the full mesh.
6.2.1 Mesh convergence study
The optimization problem is solved with the SIMP approach for four different mesh sizes: , , , and . The Young’s modulus of the void phase is set to . The design domain is initialized with a uniform material distribution of . The penalization factor is . First we consider a projection parameter of and scale the smoothing radius with the element size: the element edge length.
The optimized material distributions are shown in Fig. 13 where material with a density lower than is considered void. The strain energies are reported in Tab. 2 and display the expected decrease in strain energy as the mesh is refined. For all meshes the volume constraint is active in the converged designs.
Mesh size  Strain energy  

SIMP  40x10x10  7.5195e+03 
60x15x15  4.2076e+03  
80x20x20  4.0298e+03  
120x30x30  2.6555e+03 
As the mesh is refined, the evolution of the SIMP results shows an interesting discontinuity which is typically not observed for two dimensional problems. The optimized material layout switches abruptly from a gridtype structure, which conceptually agrees with the results of Nguyen et al (2012), to a hollow square prism design. In contrast to two dimensional structures, where refining the mesh with a meshdependent filter radius leads to an ever increasing number of holes, in this example the opposite is the case. As the element size drops below a threshold, it is more advantageous to form a continuous thin outer wall rather than a gridtype structure. This behavior is a direct consequence of the combination of SIMP penalization and density smoothing. We will revisit this issue again later.
The LSMXFEM results for a smoothing radius of the element edge length are shown in Fig. 14. No perimeter constraint is applied to this problem. Here only the results for the coarsest and the finest meshes of the SIMP study above are shown. Note that in contrast to the SIMP results, the LSMXFEM approach leads to conceptually equivalent design on both meshes. Refining the mesh only improves some local details. This feature is due to the ability of the LSM to represent thin structural features on coarse meshes. The thickness of the walls is measured by manually comparing the coordinates of two points at half the height of the design in the visualization tool. The thickness value is for the coarse mesh and for the fine mesh.
The strain energies of the LSMXFEM results are given in Tab. 3. The strain energy for the fine mesh is slightly larger than the one of the coarse mesh. This effect is due to the tendency of coarse finite element discretization over predicting the stiffness.
Mesh size  Strain energy  

XFEM  40x10x10  8.7551e+02 
120x30x30  9.8262e+02 
The differences between the SIMP and LSMXFEM results are significant. Although the discrepancy in strain energy decreases as the mesh is refined, the difference is large even for the two finer meshes where the SIMP and LSMXFEM results are similar. As the following investigation will show, the poorer performance of the SIMP results is primarily caused by the density filter, which prevents the material distribution to converge to a “01” result.
First we study the influence of the projection scheme (20) on the SIMP results for the most refined mesh. The optimized material distributions for and are shown in Fig. 15, where material with a density lower than is considered void. For convenience the result for is shown again. Table 4 reports on the strain energies and the volume fractions of intermediate densities, , as the projection parameter, , is increased. The higher , the lower and the lower the strain energy, approaching the one of the LSMXFEM result. Note that as increases, the more holes emerge. The thickness of the walls for is , which is smaller than the value for , , and closer to the XFEM value.
projection  strain energy  utilization  

SIMP  0  2.6555e+03  4.0688e02 
4  2.0264e+03  2.5131e02  
8  1.9039e+03  1.9509e02 
Instead of enforcing a better convergence toward a “01” solution by increasing the projection parameter , we postprocess the SIMP results for by the IDL postprocessing method. Figure 16 shows the strain energy of the postprocessed design for the coarsest and the finest mesh. The volume constraint is satisfied for a threshold value of for the coarse mesh, and for the fine mesh. The difference between the SIMPIDL and the SIMP results is of for the coarse mesh and for the fine mesh at the threshold . The associated strain energies are given in Tab. 5.
mesh  strain energy  

SIMP  40x10x10  0.4634  1.5601e+03 
120x30x30  0.5174  1.1530e+03 
The strain energy of the postprocessed results of the fine mesh is rather similar to the result obtained for SIMP with in Tab. 4 and the LSMXFEM results in Tab. 3. For the coarse mesh, the strain energy is well below the raw SIMP results from Tab. 2 but still above the results for the LSMXFEM approach in Tab. 3. As we will see below, this is because of the larger smoothing radius which prevents the formation of smaller features and thinner walls.
6.2.2 Feature size control
The mesh refinement study above suggests that the results of the LSMXFEM approach are less sensitive to mesh refinement than the SIMP method without meshindependent filtering. Geometric features, such as the thin walls, can be represented on coarse and fine meshes, independent of their size. This observation is in agreement with studies for twodimensional problems, see for example Kreissl and Maute (2012), but the phenomena is more pronounced and of greater importance for three dimensional problems. While this feature of the LSMXFEM approach may be considered an advantage, the ability to control the minimum feature size is of importance for many applications and to account for manufacturing constraints and costs. The following study will show that the proposed LSMXFEM approach currently lacks the ability to efficiently and intuitively control the local feature size.
We first show that applying the same absolute filter radius in the SIMP formulation efficiently controls the feature size. Figure 17(a) shows the SIMP results on the mesh for a projection parameter , a penalization factor of , and a smoothing radius of which is the same radius applied earlier for the coarsest mesh in Fig. 13(a). Comparing the SIMP results in Fig. 17(a) and Fig. 13(a) confirms the finding of numerous studies (Bendsøe and Sigmund, 2003) that the SIMP approach leads to the same conceptual layout independent of the mesh refinement level if a meshindependent filter is used. The strain energy of the design in Fig. 17(a) is given in Tab. 6.
Mesh size  Strain energy  

SIMP  120x30x30  6.5772e+03 
XFEM  120x30x30  8.2077e+02 
A similar effect is not observed in the LSMXFEM approach when we apply the same filter radius, , used for the coarse mesh to the fine mesh. Figure 17(b) shows the outcome of this procedure. The overall design is unchanged, and increasing the smoothing radius results in a less smooth design. The strain energy of this design is reported in Tab. 6.
To control the overall structural complexity in the LSM, the formulation of the optimization problem (15) is often augmented by a perimeter constraint (van Dijk et al, 2013). While this approach does not directly control the minimum feature size, reducing the maximum feasible perimeter often removes small features which do not alter much the structural performance. To study the influence of a perimeter constraint on the torsion problem, we perform the following two numerical experiments on the mesh using the LSMXFEM approach. We measure the perimeter of the SIMP result shown in Fig. 17(a) and impose this value as an upper bound on the perimeter. One problem uses the SIMP result in Fig. 17(a) as the initial design, and the other one uses the LSMXFEM result in Fig. 14(b). The results are shown in Fig. 18 and the strain energies are given in Tab. 7.
Depending on the initial designs, the LSMXFEM problems converge to different designs. While the design in Fig. 18(b) displays a trusslike design in the bottom half of the design domain, the perimeter constraint does not prevent the formation of thin walls in the upper half. The thickness of the walls in the upper half of the design is . Thus, the perimeter constraint does not control the local feature size. The design in Fig. 18(a) resembles closely the SIMP result from which it was restarted. However, considering the strain energy in Tab. 7, this design has a larger strain energy than the one in Fig. 18(b).
The study above has shown that neither smoothing the level set field nor imposing a perimeter constraint allows controlling the minimum feature size. Further, the effect of a perimeter constraint is nonintuitive as the result in Fig. 18(b) shows. The design has more structural features than the design without perimeter constraint in Fig. 14(b).
Initial design  Strain energy  

XFEM  SIMP  1.3321e+03 
XFEM  1.3185e+03 
6.3 Twophase Cantilevered Beam Design
The examples in the two previous subsections were concerned with “solidvoid” problems. Here we study a “solidsolid” problem to demonstrate the applicability of the proposed LSMXFEM approach to this class of problems. Note that the simplified XFEM formulation discussed in Section 1 is not applicable to such problems. The generalized enrichment strategy of Section 3.2 is required and the interface conditions of Eq. 5 need to be satisfied.
We study the optimal twophase layout of a cantilevered beam subject to a tip load; see Fig. 19. The stiff phase “A” has Young’s modulus of ; three values of Young’s moduli for the soft phase are considered: . Both phases have a Poisson ratio of . The maximum volume of the stiff phase is limited to of the total volume. The design domain is discretized by elements. Because of the symmetry condition, only one half of the cuboid is numerically analyzed. We compare the SIMP and LSMXFEM results.
The optimization problem is solved by a SIMP approach with a penalization factor of , a smoothing radius of ( the element edge) and the projection parameter of . The design domain is initialized with a uniform material distribution of . The optimized material distributions are shown in Fig. 20 where material with a density lower than is transparent. The strain energies are reported in Tab. 8.
The LSMXFEM results for a smoothing radius of ( the element edge) are shown in Fig. 20 and the strain energies are given in Tab. 8. Considering the full design domain, the level set field is initialized with a array of equally spaced holes with radius of . The initial design is shown in Fig. 1 and satisfies the volume constraint for the stiff phase. Note that the interface condition is enforced via the stabilized Lagrange multiplier method (5) with an element wise constant Lagrange multiplier, , and a consistency factor of .
Comparing the SIMP and LSMXFEM results, the same trends can be observed for this “solidsolid” problem as for the “solidvoid” ones studied earlier. The LSMXFEM approach leads to three dimensional structures with thinner walls and higher stiffness. In contrast, the SIMP method generates trusstype structures, in particular if the discretization is too coarse and the optimum wall thickness is less than the size of an element.
For illustration purposes only, we show a realization of the LSMXFEM optimized design for in Fig. 21. The structure was fabricated with a polyjet 3D printing process on a Connex Objet 260 printer. White material represents phase “A”, black represents phase “B”. The left and center pieces show the individual phases printed separately, the printed twophase design is shown on the right.
Density ratio  Strain energy  

SIMP  4.4081e05  
6.2862e05  
7.8627e05  
is void  7.6721e05  
XFEM  4.3221e05  
5.9192e05  
6.4448e05  
is void  6.6283e05 
7 Conclusions
In this paper we presented an optimization approach combining a level set method (LSM) for describing the geometry and an extended finite element method (XFEM) for predicting the structural response. Building upon generalized enrichment and preconditioning schemes, previously developed for twodimensional problems, the proposed optimization scheme was applied to twophase “solidvoid” and “solidsolid” problems in three dimensions. In all examples, the strain energy was minimized subject to a volume constraint on the stiff phase. The results of the LSMXFEM approach with and without perimeter constraints were compared with the ones of a SIMP method which employs density filtering and projection.
The numerical studies suggest that the LSMXFEM method features an improved convergence as the mesh is refined and is able to represent thinwalled structures on coarse meshes. The SIMP approach may require a strong projection to achieve clear “01” results with comparable strain energies. While density filtering is an efficient and intuitive method to control the local feature size, neither level set smoothing nor imposing a perimeter constraint achieve a similar effect on LSMXFEM results.
The current lack of a feature size control and the significant improved complexity of the LSMXFEM formulation limit the attractiveness of this scheme. However, for problems where a high mesh resolution is not tolerable and/or interface conditions need to be enforced with high accuracy, the LSMXFEM approach might be an interesting alternative to SIMPtype methods. The advantages of the LSMXFEM problem have been shown by Kreissl and Maute (2012) for fluid problems at high Reynolds numbers in two dimensions. The authors plan to study three dimensional flow problems with the LSMXFEM approach in the future.
References
 Allaire et al (2004) Allaire G, Jouve F, Toader AM (2004) Structural optimization using sensitivity analysis and a levelset method. Journal of Computational Physics 194(1):363–393
 Babuška and Melenk (1997) Babuška I, Melenk JM (1997) The partition of unity method. International Journal for Numerical Methods in Engineering 40(4):727–758
 Belytschko and Black (1999) Belytschko T, Black T (1999) Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering 45
 Bendsøe (1989) Bendsøe M (1989) Optimal shape design as a material distribution problem. Structural and Multidisciplinary Optimization 1(4):193–202
 Bendsøe and Sigmund (2003) Bendsøe MP, Sigmund O (2003) Topology Optimization: Theory, Methods and Applications. Springer
 Burger et al (2004) Burger M, Hackl B, Ring W (2004) Incorporating topological derivatives into level set methods. Journal of Computational Physics 194(1):344–362
 Deaton and Grandhi (2013) Deaton J, Grandhi R (2013) A survey of structural and multidisciplinary continuum topology optimization: post 2000. Structural and Multidisciplinary Optimization pp 1–38, URL http://dx.doi.org/10.1007/s001580130956z
 van Dijk et al (2013) van Dijk N, Maute K, Langelaar M, Keulen F (2013) Levelset methods for structural topology optimization: a review. Structural and Multidisciplinary Optimization pp 1–36
 Dolbow and Harari (2009) Dolbow J, Harari I (2009) An efficient finite element method for embedded interface problems. Int J Numer Meth Engng 78:229–252
 Duysinx et al (2006) Duysinx P, Miegroet L, Jacobs T, Fleury C (2006) Generalized shape optimization using xfem and level set methods. In: IUTAM Symposium on Topological Design Optimization of Structures, Machines and Materials, Springer, pp 23–32
 Eschenauer et al (1994) Eschenauer H, Kobelev V, Schumacher A (1994) Bubble method for topology and shape optimization of structures. Structural and Multidisciplinary Optimization 8(1):42–51
 Fries and Belytschko (2010) Fries T, Belytschko T (2010) The extended/generalized finite element method: an overview of the method and its applications. International Journal for Numerical Methods in Engineering 84(3):253–304
 Fries (2009) Fries TP (2009) The intrinsic xfem for twofluid flows. Int J Numer Meth Fluids 60(4):437–471
 Gain and Paulino (2013) Gain AL, Paulino GH (2013) A critical comparative assessment of differential equationdriven methods for structural topology optimization. Structural and Multidisciplinary Optimization 48(4):685–710
 Gerstenberger and Wall (2008a) Gerstenberger A, Wall WA (2008a) Enhancement of fixedgrid methods towards complex fluidstructure interaction applications. Int J Numer Meth Fluids 57(9):1227–1248
 Gerstenberger and Wall (2008b) Gerstenberger A, Wall WA (2008b) An extended finite element method/Lagrange multiplier based approach for fluidstructure interaction. Computer Methods in Applied Mechanics and Engineering 197:1699–1714
 Guest et al (2004) Guest J, Prévost J, Belytschko T (2004) Achieving minimum length scale in topology optimization using nodal design variables and projection functions. International Journal for Numerical Methods in Engineering 61(2):238–254
 Guest et al (2011) Guest J, Asadpoure A, Ha SH (2011) Eliminating betacontinuation from heaviside projection and density filter algorithms. Structural and Multidisciplinary Optimization 44(4):443–453
 Hansbo and Hansbo (2004) Hansbo A, Hansbo P (2004) A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Computer Methods in Applied Mechanics and Engineering 193(3335):3523 – 3540
 Heroux et al (2003) Heroux M, Bartlett R, Hoekstra VHR, Hu J, Kolda T, Lehoucq R, Long K, Pawlowski R, Phipps E, Salinger A, Thornquist H, Tuminaro R, Willenbring J, Williams A (2003) An Overview of Trilinos. Tech. Rep. SAND20032927, Sandia National Laboratories
 Herrero et al (2013) Herrero D, Martínez J, Martí P (2013) An implementation of level set based topology optimization using gpu. In: 10th World Congress on Structural and Multidisciplinary Optimization
 Juntunen and Stenberg (2009) Juntunen M, Stenberg R (2009) Nitsche’s method for general boundary conditions. Mathematics of Computation 78:1353–1374
 Kreissl and Maute (2012) Kreissl S, Maute K (2012) Levelset based fluid topology optimization using the extended finite element method. Structural and Multidisciplinary Optimization pp 1–16
 Lang et al (2013) Lang C, Makhija D, Doostan A, Maute K (2013) A simple and efficient preconditioning scheme for xfem with heaviside enrichments. Computational Mechancis
 Lee et al (2011) Lee P, Yang R, Maute K (2011) An Extended Finite Element Method for the Analysis of Submicron Heat Transfer Phenomena, Multiscale Methods in Computational Mechanics, Lecture Notes in Applied and Computational Mechanics, vol 55, Springer, pp 195–212
 Li et al (2012) Li L, Wang M, Wei P (2012) Xfem schemes for level set based structural optimization. Frontiers of Mechanical Engineering 7(4):335–356
 Makhija and Maute (2013) Makhija D, Maute K (2013) Numerical instabilities in level set topology optimization with the extended finite element method. Structural and Multidisciplinary Optimization
 Maute and Ramm (1995) Maute K, Ramm E (1995) Adaptive topology optimization. Structural and Multidisciplinary Optimization 10(2):100–112
 Maute and Ramm (1997) Maute K, Ramm E (1997) Adaptive topology optimization of shell structures. AIAA Journal 35(11):1767–1773
 Maute et al (1998) Maute K, Schwarz S, Ramm E (1998) Adaptive topology optimization of elastoplastic structures. Structural and Multidisciplinary Optimization 15(2):81–91
 Maute et al (2011) Maute K, Kreissl S, Makhija D, Yang R (2011) Topology optimization of heat conduction in nanocomposites. In: 9 World Congress on Structural and Multidisciplinary Optimization, Shizuoka, Japan
 Meisel et al (2013) Meisel N, Gaynor A, Williams C, Guest J (2013) Multiplematerial topology optimization of compliant mechanisms created via polyjet 3d printing. In: Twenty Forth Annual International Solid Freeform Fabrication Symposium – An Additive Manufacturing Conference, August 1214, 2013, Austin, TX
 van Miegroet and Duysinx (2007) van Miegroet L, Duysinx P (2007) Stress concentration minimization of 2d filets using xfem and level set description. Structural and Multidisciplinary Optimization 33(45):425–438
 van Miegroet et al (2005) van Miegroet L, Moës N, Fleury C, Duysinx P (2005) Generalized shape optimization based on the level set method. In: 6 World Congress of Structural and Multidisciplinary Optimization
 Nguyen et al (2010) Nguyen T, Song J, Paulino G (2010) Challenges and advances in system reliabilitybased optimization of structural topology. In: Structures Congress 2010, pp 480–491
 Nguyen et al (2012) Nguyen T, Paulino G, Song J, Le C (2012) Improving multiresolution topology optimization via multiple discretizations. International Journal for Numerical Methods in Engineering 92(6):507–530
 Ning and Pellegrino (2012) Ning X, Pellegrino S (2012) Design of lightweight structural components for direct digital manufacturing. In: 53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference
 Norato et al (2007) Norato J, Bendsøe M, Haber R, Tortorelli D (2007) A topological derivative method for topology optimization. Structural and Multidisciplinary Optimization 33(4):375–386
 Rozvany (2009) Rozvany G (2009) A critical review of established methods of structural topology optimization. Structural and Multidisciplinary Optimization 37(3):217–237
 Sethian and Wiegmann (2000) Sethian J, Wiegmann A (2000) Structural boundary design via level set and immersed interface methods. Journal of Computational Physics 163(2):489–528
 Sigmund (2001a) Sigmund O (2001a) Design of multiphysics actuators using topology optimization – Part I: One–material structures. Computer Methods in Applied Mechanics and Engineering 190(49–50):6577–6604
 Sigmund (2001b) Sigmund O (2001b) Design of multiphysics actuators using topology optimization  part II: Twomaterial structures. Computer Methods in Applied Mechanics and Engineering 190(4950):6605–6627
 Sigmund (2007) Sigmund O (2007) Morphologybased black and white filters for topology optimization. Structural and Multidisciplinary Optimization 33(45):401–424
 Sigmund and Maute (2013) Sigmund O, Maute K (2013) Topology optimization approaches: A comparative review. Structural and Multidisciplinary Optimization
 Sokolowski and Zochowski (1999) Sokolowski J, Zochowski A (1999) Topological derivatives for elliptic problems. Inverse problems 15:123
 Stenberg (1995) Stenberg R (1995) On some techniques for approximating boundary conditions in the finite element method. Journal of Computational and Applied Mathematics 63(13):139 – 148
 Svanberg (2002) Svanberg K (2002) A class of globally convergent optimization methods based on conservative convex separable approximations. SIAM J on Optimization 12(2):555–573
 Wang et al (2003) Wang MY, Wang X, Guo D (2003) A level set method for structural topology optimization. Computer Methods in Applied Mechanics and Engineering 192(12):227–246
 Wang and Wang (2006) Wang S, Wang MY (2006) A moving superimposed finite element method for structural topology optimization. Int J Numer Meth Engng 65:1892–1922
 Wei et al (2010) Wei P, Wang M, Xing X (2010) A study on XFEM in continuum structural optimization using a level set model. ComputerAided Design 42(8):708–719
 Yamasaki et al (2011) Yamasaki S, Nomura T, Kawamoto A, Sato K, Nishiwaki S (2011) A level setbased topology optimization method targeting metallic waveguide design problems. International Journal for Numerical Methods in Engineering 87(9):844–868
 Yazid et al (2009) Yazid A, Abdelkader N, Abdelmadjid H (2009) A stateoftheart review of the xfem for computational fracture mechanics. Applied Mathematical Modelling 33(12):4269 – 4282
 Zhou and Rozvany (1991) Zhou M, Rozvany GIN (1991) The COC algorithm, part II: Topological, geometrical and generalized shape optimization. Computer Methods in Applied Mechanics and Engineering 89(13):309–336