Ambiguous phase assignment of discretized 3D geometries in topology optimization

02/20/2020 ∙ by Jorge L. Barrera, et al. ∙ University of Colorado Boulder 0

Level set-based immersed boundary techniques operate on nonconforming meshes while providing a crisp definition of interface and external boundaries. In such techniques, an isocontour of a level set field interpolated from nodal level set values defines a problem's geometry. If the interface is explicitly tracked, the intersected elements are typically divided into sub-elements to which a phase needs to be assigned. Due to loss of information in the discretization of the level set field, certain geometrical configurations allow for ambiguous phase assignment of sub-elements, and thus ambiguous definition of the interface. The study presented here focuses on analyzing these topological ambiguities in embedded geometries constructed from discretized level set fields on hexahedral meshes. The analysis is performed on three-dimensional problems where several intersection configurations can significantly affect the problem's topology. This is in contrast to two-dimensional problems where ambiguous topological features exist only in one intersection configuration and identifying and resolving them is straightforward. A set of rules that resolve these ambiguities for two-phase problems is proposed, and algorithms for their implementations are provided. The influence of these rules on the evolution of the geometry in the optimization process is investigated with linear elastic topology optimization problems. These problems are solved by an explicit level set topology optimization framework that uses the extended finite element method to predict physical responses. This study shows that the choice of a rule to resolve topological features can result in drastically different final geometries. However, for the problems studied in this paper, the performances of the optimized design do not differ.



There are no comments yet.


page 18

page 20

page 22

page 24

page 25

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In topology optimization, a parametrized field that describes a problem’s geometry can freely evolve in search of an optimal design. The level set method ([1, 2, 3]) provides a convenient framework for this process. Traditionally, immersed analysis techniques smear the interface over elements (e.g., the Ersatz material method [4]). However, a rapidly growing class of these techniques are able to capture and track the interface explicitly (e.g., cut finite element method [5], extended finite element method [6]). In such techniques, embedding the geometry in a computational domain via a discretized level set field (LSF) may allow for ambiguous definitions of the discretized geometry at any stage of the design optimization process. Therefore, methods for detecting and resolving these ambiguities are needed. Yet, with only a few exceptions ([7, 8, 9]), work on level set topology optimization omits mentioning and addressing this issue.

Ambiguities of discretized geometries can be resolved via user-defined rules. The effect of a given rule on the optimization process and the optimized designs, i.e. whether some user-defined rules can lead to suboptimal designs, is still an open question. This scenario is illustrated in the 2D schematics of Fig. 1, where two different rules, and , are used to resolve geometric ambiguities. The sketch shows that, starting from the same initial design, not only different intermediate designs are possible, but also that the optimized designs may differ in geometry and performance. Hence, in this work we investigate the influence of these rules used to determine ambiguous geometrical configurations in the context of topology optimization.

Figure 1: Probable effect of phase assignment using either rule or to resolve ambiguous geometrical features in the optimization process.

In this paper, we focus on level set topology optimization approaches that use fixed background meshes to discretize the geometry and/or the state variable fields. These approaches typically split elements in which the interface is embedded ([3, 10]). The intersected elements are decomposed into sub-elements such that a face in 3D, or an edge in 2D, aligns with the interface defined by the isocontour of the LSF. A solid phase needs to be assigned to each of the generated sub-elements. However, certain geometrical configurations suffer from ambiguous definitions of the interface, and thus the phase cannot be uniquely determined for all sub-elements in these cases. This scenario is displayed in the zoomed-in region of the 2D embedded topology in Fig. 2. The level set values at the corner nodes of the center background element allow for two different interpretations of the geometry. Although both configurations are consistent with the discretized LSF, one of them creates connections between the void inclusions.

Figure 2: Topology consistency issue in 2D.

In a 2D domain discretized by quadrilateral elements, ambiguous topological features exist only if an element is intersected by an isocontour more than once. In 3D domains composed of hexahedral background elements, this issue is also found in simpler intersection configurations that describe a single interface, and is aggravated if multiple interfaces are discretized within an element ([11, 12, 13]). Furthermore, special care may be needed for neighboring background elements that contain sub-elements with ambiguous phases assignment. Inadequate handling of sub-elements that share common faces across adjacent background elements may introduce inconsistencies in the geometry description ([8]).

Although ambiguous geometrical configurations are avoided if a LSF is discretized using tetrahedral elements, hexahedral elements remain the preferred choice ([14]). This is in part because hexahedral elements provide a good compromise between low computational costs and accurate approximations. Furthermore, the need of mesh refinement is reduced in hexahedral elements with complex intersection configurations (i.e., multiple interfaces) if higher order shape functions and/or B-splines are used to adequately predict physical responses ([15]).

Geometrical aspects of the aforementioned ambiguities have been studied in applications ranging from medical data visualization to computer graphics problems. Techniques developed to resolve ambiguous topological features have been formulated for a wide range of requirements such as preserving topology, filtering noise, pleasing visual appearance, among others (

[16, 17, 9]). Some of these techniques seek to improve continuity under continuous deformations, which can restrict the evolution of the geometry in an optimization process. For example, a class of methods that aim at preserving the topology perform nodal analyses to identify critical nodes that generate ambiguous configurations. In these techniques, only shape changes are allowed by preventing sign changes in the level set value of these critical nodes; see [18, 19]. Applying them in a level set topology optimization problem would prevent the merge/split of geometric features. Hence, not all techniques for resolving ambiguous topologies are compatible with the requirements of topology optimization.

Freedom in the evolution of a smooth topology can be achieved using methods with various degrees of complexity. If the generation of the interface is viewed as a boundary-value problem, the solution of a system of partial differential equations can be used to define a smooth interface (

[20, 21, 22, 23]). Here, a set of parameters control the geometry through boundary conditions along the interface. Implementing these methods introduces significant additional complexities associated with building a new system of equations and mapping the resulting surface onto the discretized domain. A simpler approach consists of locally refining intersected subdomains such that a background element is only intersected once ([24]). These methods reduce the complexity of the intersection configurations and improve geometric resolution. Yet, ambiguities may still exist, and the computational framework employed requires handling hanging nodes in the case of hierarchical mesh refinement. As an alternative, intersected elements can be discretized using boundary conformal polyhedral instead of tetrahedra; see [8]. Inter-element ambiguities in this case are resolved by the so-called asymptotic decider. This technique uses the nodal level set values of the hexahedron shared face to linearly interpolate the level set value at the saddle point; see [25] for an in-depth explanation. Though, implementing this approach requires the use of non-traditional finite elements, i.e. polyhedral elements. Another way to circumvent ambiguous topological features consists of subdividing intersected hexahedra into structured tetrahedra, and discretizing the LSF on these generated sub-elements (e.g., marching tetrahedra [26, 27]). However, not only a considerably larger number of tetrahedra are required in these methods, but the approximated LSF may compromise smoothness in the resulting topologies ([28]).

Abstaining form mesh refinement methods, techniques used to resolve topological ambiguities either resort to geometric indicators or user-defined criteria for assigning phases to sub-elements. In this case, a phase is assigned to tetrahedralized hexahedra based on measurements of area, curvature, etc., as shown in, for example, [29, 30, 12]; or by consistently selecting a phase for all ambiguous sub-elements. This class of methods is explored in this paper.

The study presented here focuses on analyzing topological ambiguities in the description of embedded geometries from discretized LSFs on hexahedral meshes. Ambiguous topological features are categorized, and a set of rules to determine a phase in ambiguous tetrahedral sub-elements based on either user-defined or geometric indicators is introduced. Local (element-wise) and global (inter-element) rules are tested using different optimization formulations to assess their influence on the evolution of the geometry of the optimization process and performance of the optimized design. Two-phase, solid-void design optimization problems are considered and are solved with an explicit level set topology optimization framework. The extended finite element method (XFEM) is used to discretize the weak form of the governing equations; see [6]. A linear elastic behavior is considered in all structural design problems.

The remainder of the paper is organized as follows: Section 2 details the characteristics of the geometrical configurations affected by the topological ambiguities. In Section 3, a set of rules for resolving such ambiguities, as well as simple algorithms as guidelines for implementation, are provided. The level set XFEM optimization framework used in this study is summarized in Section 4. The influence of the rules introduced in Section 3 is examined with structural design problems in Section 5. Finally, Section 6 summarizes this paper and provides directions for future work.

2 Topology Ambiguities in Level Set Methods

2.1 Geometry description

The layout of a two-phase solid-void problem in a design domain, , is described by a LSF, , as follows:


The design domain is subdivided into the solid and void phases, and , respectively, such that . The interface, denoted by , is identified by the zero level set isocontour.

The LSF is discretized by tri-linear shape function on a structured mesh,

. This mesh is composed of regular eight-node hexahedra that are classified into two groups,

and , with . The subdomain contains all the un-intersected hexahedra. A hexahedron, , is part of this set if all their nodal level set values are either positive or negative. Furthermore, the subset consists of intersected hexahedra defined as all hexahedra for which at least one of their eight nodal level set values differs in sign from the remaining nodal level set values.

All intersected hexahedral elements, i.e. , are split into tetrahedral sub-elements via a tetrahedralization algorithm. The vertices of the tetrahedra are constructed from vertices of the background element and intersection points along the edges. Hence, some of the faces of these generated tetrahedra contain intersection points only. Tetrahedralization algorithms employed for this subdivision receive the nodal level set values and the intersection points of a hexahedron as input and return the geometrical information of the set of generated tetrahedra, . In this work, all the generated tetrahedra are contained in the set . The Delaunay method ([31, 32]) is used for tetrahedralizing the intersected hexahedral elements in this study.

The intersection points along the edges are defined by the points with a zero level set value. The level set values along an edge are linearly interpolated from the level set values of the edge’s corner nodes. Thus, the linear interpolation of the LSF limits the possible configurations to one intersection per edge. Since each node can have a positive or negative level set value, a total of nodal phase assignment configurations (NPACs) are possible in an eight-node hexahedral element. After accounting for reflective and rotational symmetries, the number of unique NPACs for intersected elements is reduced to the 14 cases shown in Fig. 3. The reader is referred to [11, 12, 13] for a case-by-case analysis, as well as detailed visual schematics.

Figure 3: All unique NPACs using a discretized LSF onto an intersected eight-node hexahedron.

In this work, the discretized geometry is constructed by connecting intersection points along intersected edges. Thus, a small discrepancy between the discretized geometry and the approximated level set field using linear interpolation functions is expected. The left side of Fig. 4 illustrates this in a simple 2D schematic. It can be seen that the zero isocontour of the LSF, described by a blue dashed curve, is approximated by a straight blue line between two intersected edges.

Figure 4: Two contiguous hexahedra which discretized geometry present ambiguous phase assignment (center). The highlighted faces show (left) the difference between approximated level set field and discretized geometry, and (right) level set field contours on such face and its saddle point.

Once intersected hexahedra have been split into tetrahedra, one or multiple interfaces must be defined to approximate the problem’s geometry. This implies that either a solid or a void phase must be assigned to every tetrahedron, . In some cases, however, one or more tetrahedra can be defined as solid or void while keeping consistency with the associated NPAC. The lack of uniqueness inherent in such cases leads to multiple possible interface definitions. In other words, different discretized domains can be generated from the same discretized LSF depending on the phase assigned to these ambiguous tetrahedra, as shown in Fig. 2 in 2D.

2.2 Ambiguous tetrahedra (ATs)

An ambiguous tetrahedron, AT, is defined as any tetrahedron whose vertices are all intersection points. Hence, a minimum of four intersection points is needed for an AT to exist. For this reason, all but the NPAC-1 configuration can contain ambiguities; see Fig. 3. Despite the number and shape of the generated tetrahedra may change depending on the tetrahedralization algorithm employed ([32, 33]), ATs cannot be avoided.

Figure 5: Illustrations of the four cases of intersection configurations for a tetrahedralized hexahedron: (a) no ambiguous sub-elements; (b) IATs with one interface; (c) IATs with multiple interfaces; (d) IATs and BATs with multiple interfaces.

To simplify the analysis of topology ambiguities, the NPACs described in Fig. 3 are categorized as shown in Fig. 5. Examples of continuous LSFs to be embedded in a discretized domain are shown on the left column. The discretized geometries are shown in dark gray, together with the ATs highlighted in light gray in the center column. The right column shows only the ATs for better visualization. In the first case, illustrated in Fig. 5(a), no ATs are present since only three intersection nodes exist. Fig. 5(b) shows the most innocuous NPAC with ATs, which occurs in hexahedra with a single interface. In this configuration, assigning different phases to the ATs results in different shapes but not topologies. If the background element is intersected more than once, as shown in Fig. 5(c), then ATs can connect/disconnect interfaces, and thus alter the topology. Note that in both Figs. 5(b) and (c), none of ATs faces is part of the background hexahedron faces. We classify them as Internal Ambiguous Tetrahedra, herein labelled as IATs. However, NPACs can also include ATs with one of their faces overlapped with one of the hexahedron’s faces, as is depicted in Fig. 5(d). We term these Boundary Ambiguous Tetrahedra as BATs.

2.2.1 Internal ambiguous tetrahedra (IATs)

A key characteristic of this type of ambiguities is that at least one of the faces of an IAT lays on an interface. If the NPAC leads to only one interface, all ambiguous sub-elements are IATs, as shown in Fig. 5(b). Furthermore, the difference between the interface geometry approximated by the IATs and the one described by the linearly interpolated LSF is typically negligible.

2.2.2 Boundary ambiguous tetrahedra (BATs)

This type of AT can only exist in intersected hexahedra that are able to describe multiple interfaces. Furthermore, BATs are always accompanied by IATs. A well-known issue of NPACs with BATs is guaranteeing consistency in phase assignment across two adjacent hexahedral background elements that share a face on which a BAT is defined; see [25, 8]. The sketch at the center of Fig. 4 shows this scenario. To achieve a consistent phase assignment, BATs sharing a hexahedron face must have the same phase.

2.2.3 Effect of ambiguous tetrahedra

Internal ambiguities are ubiquitous and, in most cases, harmless since they do not alter significantly the geometry. However, hexahedral background elements with multiple interfaces can have IATs which phase assignment may have a strong influence on the topology of a given design. In most of those scenarios, the analysis has to include BATs.

Inappropriate treatment of ambiguous configurations that combine IATs and BATs can result in undesired geometries, especially when thin features are present in the design. Shell-like structures, which are not uncommon in engineering design problems, can be greatly affected since the appearance of dents may be favored. This is illustrated in the example in Fig. 6. Here, the influence of ATs on the geometry representation of a spherical shell is shown for three different wall thicknesses. The shell has a constant inner radius and its thickness is reduced by varying the external radius. An eighth of the domain is discretized by a hexahedral background mesh and is highlighted in dark gray for visualization purposes. In this example, void phase is assigned to all ATs. The discretized geometry, together with the ratio between volumes of AT and solid subdomain, are provided for the three configurations on the right side of Fig. 6. It can be observed that once the difference between the inner and outer radii is smaller than the edge length of a background element, the volume of the ambiguous subdomain increases up to 26.71% of the solid subdomain. Note that the extreme misrepresentation of the spherical shell geometry can be avoided by assigning a solid phase to all ATs. However, assigning the same phase to ATs is not necessarily the most beneficial option for all intersection configurations, as is shown in Section 5.

To assess whether it is more beneficial in the optimization process to use a simple criterion, such as the one used for the spherical shell, or geometrical information, a group of rules that systematically assigns a phase to AT is presented in next section.

Figure 6: Spherical shell with varying thickness in an immersed boundary geometrical description. An eighth of the domain is highlighted in dark gray to enhance contrast.

3 Rules for Resolving Topology Ambiguities

The ambiguities described in Section 2.2 can be resolved by applying phase assignment rules either locally (i.e., separately for each hexahedral background element) or globally (i.e., on groups of adjacent hexahedra). The local and global rules explored in this paper are denoted by and , respectively. These rules collect a subset of techniques that allow for the free evolution of the geometry in a shape or topology optimization process, but this subset is by no means complete. The reader is referred to [14, 25, 27, 30] for more information on rules used to this end.

3.1 Local rules ()

The local rules described in this section resolve IATs and BATs separately, using one criterion for each ambiguity type. Consistent phase assignment across hexahedra (i.e., BAT phase assignment) is achieved by enforcing the asymptotic decider [25]. This technique uses the nodal level set values of the hexahedron face containing a face of a BAT to compute the level set value at the saddle point, , assuming a bi-liner interpolation of the level set field over the hexahedron face. The right side of Fig. 4 shows a schematic of the location of on a face that contains BATs. The level set value at the saddle point is computed as follows:


where , , , and correspond to the nodal level set values of the sharing hexahedron face. The sign of determines the phase to be assigned. If , then all BATs in the adjacent hexahedra sharing the face are assigned the solid phase. Conversely, if , BATs are assigned the void phase. The reader is referred to [25] for details of this technique.

After BATs have been resolved via the asymptotic decider, the phase of the remaining IATs is determined using one of the rules described below (in order of increasing complexity in terms of implementation).

3.1.1 Fixed phase assignment (/)

In this rule, all IATs are assigned the same phase throughout the optimization process. If ambiguous sub-elements are assigned the solid phase, the rule is denoted as ; conversely, if void phase is assigned, it is named . This approach requires no additional computations for phase determination, as can be seen in Algorithm 1.

// Loop over all intersected hexahedra of the background mesh
for  do
       // Check if any of the generated tetrahedral sub-elements are ambiguous
       if  then
             // Assign predefined phase to all IAT
Algorithm 1 Resolving IATs using the scheme

3.1.2 Previous main phase in optimization process ()

An easy-to-implement alternative that promotes a smoother evolution of the design during the optimization process is choosing the phase assigned to the IATs in the current design iteration, , based on an elemental phase, , function of the previous design iteration, . If, for example, a hexahedron with the configuration NPAC-4 shown in Fig. 5(b) is part of the void phase at , this rule prevents connecting the two interfaces by setting the IATs to void phase. In this scheme, all internal ambiguities are assumed to be solid phase in the first design iteration (i.e., at ) to promote a topology preserving evolution of the design at the initial stage of the optimization process. An implementation of this logic is detailed in Algorithm 2.

// Loop over all hexahedra of the background mesh
for  do
       // Check if element is not intersected to initialize elemental phase
       if  then
             // Store phase of un-intersected hexahedron
             if  then
                   // If this is the first design iteration, set elemental phase to solid phase
             // Assign phase to the ambiguous sub-elements
             if  then
                   IATs phase =
Algorithm 2 Resolving IATs using the scheme

3.1.3 Based on level set value at averaged centroid ()

Expanding the idea of the asymptotic decider to 3D, IATs can be resolved based on the sign of a single interpolated level set value in . The decider point (equivalent to the saddle point of the asymptotic decider) is defined as the average of the centroids of the IAT within a background element. This is illustrated in 2D in Fig. 7, where the centroid of the ambiguous subdomain, , is used to determine a phase. Note that information of the previous optimization step is not needed in this rule. Also, in hexahedra with a single interface, the phase is assigned in accordance to the unambiguous LSF used to construct the interface. A scheme that implements this rule is shown in Algorithm 3.

Figure 7: 2D schematics of the rule.
// Loop over all intersected hexahedra of the background mesh
for  do
       // Check if hexahedron contains ambiguous sub-elements
       if  then
             // Compute centroids of every IAT using the coordinates of their vertices,
             // Compute the average of the centroids of the ambiguous sub-elements
             // Interpolate the level set value at centroid using linear shape functions, , and the nodal level set values, , of the containing hexahedron
             // Assign elemental phase based on the sign of the level set value
             if (
             // Assign a phase to the ambiguous sub-elements
             IATs phase =
Algorithm 3 Resolving IATs using the scheme

3.1.4 Shared area between Iat and neighboring tetrahedra per phase (/)

To avoid interpolating the LSF to resolve the IATs, a rule based on a geometric indicator is proposed. In this rule, the shared surface area between each IAT and its unambiguous neighboring tetrahedra are computed per phase. The shared surface areas of two neighboring IATs are omitted from this calculation. The phase of the unambiguous tetrahedra with whom the IATs share the maximum () or minimum () surface is chosen for all IAT. Figure 8 shows a 2D schematic of this rule.

Note that implementing this scheme requires knowledge of the phase of adjacent tetrahedra. Hence, the connectivity of the generated tetrahedra must be constructed per intersected hexahedron. An approach that employs this rule is summarized in Algorithm 4.

Figure 8: 2D schematics of the rule.
// Loop over all intersected hexahedra of the background mesh
for  do
       // Check if hexahedron contains ambiguous sub-elements
       if  then
             // Initialize areas for both phases
             // Loop over all ambiguous tetrahedra in intersected hexahedron
             for   do
                   // Loop over the faces, , of current ambiguous tetrahedron
                   for  do
                         // Get tetrahedra sharing current face
                         and sharing
                         // Skip face computations if both tetrahedra are ambiguous.
                         if  ()  then
                               // Compute face area,
                               // Update area of shared tetrahedra based on the phase of unambiguous tetrahedron
            // Assign elemental phase based on shared area
             if   then
                   for , or for
                   for , or for
            // Assign a phase to the ambiguous sub-elements
             IATs phase =
Algorithm 4 Resolving IATs using the scheme

3.2 Global rules for resolving topology ambiguities ()

In this set of rules, the asymptotic decider is no longer used to enforce consistent phase assignment to BATs across intersected hexahedra. Here both IAT and BAT are resolved using the same criterion. Consequently, chains of connected ATs across multiple elements can be formed. The criteria described below are applied either on the entire domain () or per cluster of connected ambiguous tetrahedra connected across intersected tetrahedra ().

3.2.1 Fixed phase assignment (/)

This rule enforces either solid (), or void () phase to all ambiguous tetrahedra. Similar to , no additional implementation is needed once the ATs have been identified. However, in contrast to , here BATs are also included. Note that for this rule the choice of using or may have a major impact on the resulting geometry compared to the other rules presented in this paper.

3.2.2 Shared area between ambiguous and uniquely defined tetrahedra per phase (/)

This global rule is an extension of the local rule and is enforced on clusters of connected tetrahedra. These clusters can exist across intersected hexahedra instead of being contained within a single hexahedron; see Section 3.1.4. Hence, knowledge of the connectivity of all tetrahedra, which changes for every design iteration, is needed to identify these clusters. A graph search scheme to generate this connectivity that avoids issues associated with memory and/or computational time consumption is provided in Algorithm 5. Note that implementing this rule is laborious, especially in parallelized software. Existing third party libraries, such as MOAB ([34]) or the STK package that is part of Trilinos ([35]), provide robust frameworks to manage mesh databases and can be used to ease the search of connected ATs.

// Initialize container of clusters, = []
// Initialize flag for tetrahedra of intersected background elements, : = false,
// Loop over all intersected hexahedra of the background mesh
for All tetrahedra faces,  do
       // Get tetrahedra sharing current face, ,
       // Check if any tetrahedra sharing is ambiguous and has not been processed
       if () & ()  then
             // Initialize a face list with current face,
             // Initialize an empty cluster of AT, = []
             // Process tetrahedra in face list
             while  do
                   // Loop over tetrahedra sharing current face,
                   for  do
                         // If tetrahedron is ambiguous, process it
                         if  then
                               // Add tetrahedron to current cluster
                               = [ , ]
                               // Loop over faces in current tetrahedron,
                               for   do
                                     // Add faces to list
                              // Flag tetrahedron as processed to avoid repeating searches in first loop
                  // Remove current face in list
       // If a new cluster was generated from the processed tetrahedra
       if  then
             // Add new cluster to container of clusters
             = [, ]
// After all clusters of connected ATs have been generated, loop over them
for  do
       // Initialize areas for both phases:
       // Loop over faces of current cluster
       for  do
             // Get tetrahedra sharing current face,
             // Check if only one of the shared tetrahedra is ambiguous
             if Only one  then
                   // Compute face area,
                   // Update area of phase of shared tetrahedra (either or )
      // Assign phase based on shared area
       if   then
             for , or for
             for , or for
      // Loop over all ambiguous sub-elements in current cluster
       for  in  do
             Assign phase
Algorithm 5 Resolving ATs using the scheme

3.3 Local versus global rules

In terms of implementation, local rules , and , and global rule , are more attractive as they require little effort. Additionally, these rules have a low computational overhead. However, rules and provide a more consistent/smooth geometry description at the cost of higher computational effort. The advantage of these last two rules comes at the cost of involved implementations of problem specific geometric indicators to resolve ambiguities.

Local rules , and allow for imperfections in the geometry description due to disparity in phase assignment between BAT and IAT. Either solid floating pieces embedded in a mostly void subdomain, or small holes in a mostly solid subdomain, are likely to occur in hexahedra with multiple interfaces, see [28]. In the former case, the numerical technique used for the physical analysis must include appropriate treatment of these disconnected pieces, such as consistent generalized enrichment when using XFEM and adding stiffness to avoid zero energy modes; see Section 4.3. The local rule does not suffer from the issues described above since the phase of the interior clusters is determined from the phase of the surrounding unambiguous tetrahedra.

Global rules have a more pronounced effect on the resulting discretized geometry compared to local rules. Although only rarely observed, assigning different phases to large clusters of connected ATs can significantly alter the topology. Differences in topologies are less evident if local rules are used because of the asymptotic decider. The influence (if any) of using the rules described in this section on the smoothness of the evolution of the geometry in a topology optimization process, as well as the performance of the optimized designs, has been unexplored to date. In this paper, a study of these rules is presented in Section 5.

4 Explicit Level Set Topology Optimization Framework

The design, i.e. the geometry, is controlled by a vector of design variables,

bounded between user-defined lower, , and upper, , bounds. A design variable is assigned to each node of the structured hexahedral background mesh; thus, equals the number of nodes of the background mesh; see [3, 10] for details.

To increase numerical stability and enhance convergence of the optimization problem, the explicit relationship between the design variables and the nodal coefficients of the discretized LSF is defined by the linear filter described in [36]. This distance-based filter is defined as:


where is the number of nodes within the filter radius ; and is the Euclidean distance between nodes i and j. Within each hexahedral background element, the LSF is interpolated tri-linearly. Furthermore, the regularization scheme described below in Section 4.1 is adopted to promote a uniform spatial gradient of the LSF near the solid-void interface while converging to either a positive or negative target value away from the interface.

4.1 Optimization problem formulation

The optimization problems considered in this work are formulated as:


The first component of the objective, , represents the quantity of interest to be minimized, , such as strain energy, target displacement, and mass. The second term is the normalized perimeter control penalty, added to avoid the emergence of irregular geometric features:


which is normalized by the perimeter of the initial design, . The last component is a penalty that aims at regularizing the LSF. This penalty ensures (i) an approximately uniform spatial gradient near the solid-void interface; and (ii) that the LSF away from the interface convergence to target upper and lower bounds. The normalized formulation of reads:


with , , and as weighting factors. The first term penalizes the difference between the LSF, , and a target value, , away from the interface. The second term penalizes the difference between the norm of the spatial gradient of the LSF, , and a target gradient norm, , in the vicinity of the interface. The third term enforces a gradient of zero away from the interface. The parameter is a measure of the distance of a node to the interface defined as:


where the parameter determines the proximity of the interface. In this regularization scheme, the distance of a node to the interface is approximated by a neighborhood-level field , as described in [15], and illustrated in Fig. 9. The parameter represents a measure of the of the region along the interface over which the approximation is evaluated. The higher , more background elements surrounding the interface are considered. In this work, we set , , , and the weighting factors in the range of . Furthermore, the weights in the objective function are selected such that the contributions of the regularization and perimeter control penalty components are significantly lower than .

Figure 9: Neighborhood-level field, , using .

The design needs to satisfy a set of problem dependent inequality constrains, [, …,

] that consider limits on, e.g., mass, maximum allowable stress, and minimum eigenvalue. The number of design variables are denoted by

and the number of state variables is ; see Eq. 4. The state variables are governed by a set of discretized partial differential equations, which are satisfied at each design in the optimization process.

The parameter optimization problems are solved using a gradient-based nonlinear programming scheme, namely the Globally Convergent Method of Moving Asymptotes (GCMMA), with no inner iterations; see [37].

4.2 Sensitivity analysis

The design sensitivities of the objective function with respect to the design variables are defined as


The first term of the right-hand side represents the explicit dependencies. The second term accounts for the implicit dependencies through the state variables. The implicit term is evaluated using the adjoint method:


where is the adjoint response and is the residual defined later in Eq. 13. The adjoint linear problem is formulated as:


For more information on the sensitivity analysis within a level set XFEM topology optimization framework, the reader is referred to [38].

4.3 Physical analysis

A generalized Heaviside enriched XFEM

strategy where the degrees of freedom within each unique sub-phase are approximated using standard finite element shape functions is adopted. The flexibility of this

XFEM approach has been demonstrated in single and multiphysics topology optimization problems; see for example, [39, 40]. Here, the response is consistently interpolated in elemental subdomains with the same phase using the formulation shown below.

4.3.1 XFEM formulation

The state variable vector, , at node i of a two-material problem (phase I, phase II) is approximated as:


where is the Heaviside function depending on the LSF:


In Eq. 11, denotes the maximum number of enrichment levels, represents the nodal shape function, and is the Kronecker delta which selects the active enrichment level for node . The Kronecker delta ensures that displacements at node are only interpolated by a single set of degrees of freedom defined at node position such that the partition of unity principle is satisfied ([41]). The number of nodes per element is given by . The reader is referred to [42, 43, 44] for an in-depth explanation of the generalized Heaviside enrichment strategy employed.

Note that other immersed boundary techniques with a crisp interface definition, such as CutFem [5] and HIFEM [45], are equally applicable without loss of generality.

4.3.2 Governing equations

Static equilibrium is described by the following residual equation augmented with stabilization terms:


in which the weak form of the linear elastic governing equation, , is defined as:


The displacement field and the test function are denoted by and , respectively. Traction forces, , are applied along the boundary

. The material tensor for isotropic linear elasticity,

, together with the infinitesimal strain tensor, , define the Cauchy stress .

The unsymmetric version of Nitsche’s method is employed to weakly enforce Dirichlet boundary conditions ([46, 47]). These conditions are applied through the following residual component:


The normal vector on a domain boundary is denoted by . The jump operator is defined as


The integrals of Eq. 15 correspond to the standard consistency, adjoint consistency, and the Nitsche penalty terms, respectively. The third term is scaled by the Young’s modulus, , the element size, , and the penalty factor . This last term provides additional control over the accuracy at which a boundary condition (BC) is enforced. Our numerical experiments show that choosing a penalty factor within a range of provides an adequate enforcement of BCs while avoiding an ill-conditioned linear system.

To prevent numerical instabilities due to vanishing zones of influence of certain degrees of freedom when the interface approaches a node, face-oriented ghost penalization is used in the vicinity of the interface; see [48] and [5]. Approximating and by tri-linear shape functions, ill-conditioning is mitigated by applying the following penalty:


Element faces in the vicinity of the interface for which at least one of the two adjacent elements is intersected are included in , as explained in [49]. Normals of these element faces are denoted by . The influence of the ghost penalty term presented above is controlled by the penalty factor , which is usually chosen to be within based on numerical experiments; see [50, 51].

To suppress rigid body motion of disconnected solid subdomains that may emerge and develop, selective structural springs are added to the governing equations; see [50, 52]. An additional stiffness term is applied only to solid subdomains that are disconnected from the support via the following residual component:


The parameter denotes the spring stiffness scaling and is non-zero only for the free-floating pieces of solid material. These pieces are identified by an auxiliary indicator field constructed from the solution of a diffusion problem, as described in [52].

In the final example, a gradient stabilized scalar stress field, , is post-processed via the XFEM informed smoothing procedure described in [53]; see Section 5.3. The additional set of state variables is computed by solving the residual equation:


The field and the parameter represent the von Mises stress and the ghost penalty weight, respectively. Spatial oscillations of are mitigated by penalizing the jump in spatial stress gradients across elemental faces (second term in Eq. 19).

5 Numerical Examples

The local and global rules presented in Section 3 are studied in this section to assess their influence on the evolution of the topology, and in the optimized design. Furthermore, the influence of the initial seeding pattern of holes on the optimization process is investigated for different rules. These rules are tested on structural linear elastic single material, solid-void problems solved by the explicit level set XFEM topology optimization framework detailed in Section 4. A strain energy optimization problem is formulated in Example 1, minimizing the displacement in a portion of the structure is investigated in Example 2, and in the last example, the objective consists of minimizing the total mass while satisfying a stress constraint. The optimization problem is considered converged if the constraints are satisfied and the relative change in objective between two consecutive design iterations is less than . The initial seeding of the design domains with holes is constructed using patterns of cuboid or spherical primitives.

The governing equations are discretized using the XFEM approach outlined in Section 4.3. Relevant parameters used for the numerical examples 1 and 2 are listed in Table 1 in self-consistent units. Table 5 summarizes the parameters used in the last example in SI units. The example problems consist of a one-way coupled set of governing equations; i.e., the diffusion problem describing the auxiliary indicator field, the stabilized linear elasticity equations (Eq. 13) and, for the third problem, the stabilized stress projection equations (Eq. 19). The systems of equations of the first and second examples are solved using an ILU preconditioner ([54]). The third example was solved iteratively via a Trilinos algebraic multigrid solver ([55]).

Parameter Value
Target LSF
LSF regularization control
Filter radius
Nitsche penalty factor
Ghost penalty factor
Spring stiffness factor
Stress ghost penalty factor (Ex. 3)
Table 1: Parameters common to all numerical examples function of the element size .

5.1 Beam

A beam problem is used to investigate the effect of the topology consistency rules detailed in Section 3 as the mesh is refined. In this example, we focus on analyzing differences in the optimized designs that typically favor thin-walled shear-webs when using a level set topology optimization approach. Three levels of refinement are considered.

The problem setup is shown in Fig. 10. The design domain of size is simply supported at all four corners of the bottom face, as shown on the left side of Fig. 10. A traction load is applied on the center region of the top face on a area. Only one quarter of the design domain is simulated by applying symmetry boundary conditions along the and planes at the center of the design domain. The loading and support regions highlighted in dark gray are excluded from the design domain. Structured hexahedral meshes with element sizes of are considered. The initial design for this problem is constructed using cuboids as void inclusions. The right side of Fig. 10 shows the hole seeding arrangement for the finest mesh used in this example.

Figure 10: Beam problem: (left) design domain with boundary conditions, symmetry planes and dimensions; and (right) initial design with void inclusions.

The optimization problem seeks to minimize the strain energy, , subject to a mass, , constraint; and is formulated as follows:


The first component of the objective is normalized by the strain energy of the initial design, denoted by . The two remainder components are described in Section 4.1.

In absence of proper hole seeding mechanisms (e.g., using topological derivatives ([56, 57, 4]) or a density field to nucleate holes in regions of low densities ([58])), the optimization process can converge to sub-optimal designs. For this reason, continuation schemes are typically used to gradually enforce constraints. In this example, a mass constraint with is enforced using the continuation strategy described in Table 2 to avoid premature removal of mass resulting in suboptimal designs. This constraint is reduced every 12 design iterations in a total of 8 steps (i.e., for ). The mass of the design is normalized by the mass of the entire design domain, . The optimization problem starts feasible for all mesh sizes, i.e., the initial hole seeding satisfies the initial mass constraint. The weights of the objective are =[ 0.85, 0.05, 0.10].

0.66 0.58 0.50 0.42 0.35 0.30 0.25 0.20
0.40 0.36 0.33 0.30 0.27 0.24 0.22 0.20
0.27 0.26 0.25 0.24 0.23 0.22 0.21 0.20
Table 2: Continuation parameters of beam problem for and three mesh sizes.

5.1.1 Effect of ambiguity decider on optimization problem

The evolution of the strain energy, , and mass constraint, , for the results using the finest mesh (i.e., ) are shown in Fig. 11. The sudden violations of the constraint, and to a lesser extent in the objective, are consequence of the mass constraint being updated at every continuation step. Overall, a smooth behavior is observed in all cases. The strain energies of the optimized designs are summarized in Table 3 and show that the relatively small differences in performance are reduced as mesh is refined. In this example, the rules used to resolve ambiguities have a negligible effect on both the smoothness of the evolution of the design, as well as the performance of the optimized designs.

Figure 11: Evolution of (left) objective and (right) mass constraint of beam problem using all ambiguity decider rules.
Mesh size
Table 3: Final strain energy of beam problem using three different mesh sizes.
Figure 12: Optimized designs of selected rules using meshes with element size .

The optimized designs shown in Fig. 12 for all three refinement levels and , , , and rules show that differences in final topologies are imperceptible in this example. Although not shown, almost identical designs where obtained using the remaining rules. Continuous walls connecting the two flanges at the top and bottom of the beam with a thickness of about the element edge length are observed. Ambiguous tetrahedra, highlighted in red, are observed throughout the optimization process and stem exclusively from hexahedra with one intersection. Regardless of the topology consistency rule employed, the thickness of the shear webs in the final structures is slightly larger than the element lengths. Hence, holes in such walls that may be consequence of resolving ambiguities using some of the rules described in this paper are not observed in the optimized designs. The lack of ambiguities that can create connections throughout the evolution of the design explains the little variations in the optimized designs.

Figure 13: Evolution of ambiguous elements throughout the beam optimization problem: (left) ratio between volume of all ATs and total volume of solid phase; and (right) ratio between volume of ATs defined as solid and total volume of solid phase.

The evolution of the ratio of ATs with respect to the total solid volume, , is plotted on the left side of Fig. 13. A slight increase in the disparity between rules of these ratios is observed as the optimization problem converges. Nevertheless, the AT volume fraction remains under in all cases. The plot on the right of Fig. 13 shows the ratio of the volume of ATs for which a solid phase was assigned with respect to the total solid volume, . As expected, global rules and provide upper and lower limits for the cases considered; see Section 3.2.1. Furthermore, it can be seen that rules , , , and tend to distribute the phase assignment more evenly between the void and solid subdomains. In all cases, the volume of the ATs is, at all times, considerably smaller than the uniquely defined domain.

5.1.2 Shear web structures and mass constraint

To demonstrate that the avoidance of ambiguities along the shear webs are not specific for the particular mass constraint studied above, lower target mass constraints, i.e., , are examined next. The only difference with respect to the previous problem setup are the continuation parameters; see Table 4. The analysis for this modified setup focuses on rules and and a mesh size of since they provide the widest variations in the problem’s geometry at any particular design iteration.

0.27 0.25 0.23 0.21 0.19 0.17 0.16 0.15
0.27 0.24 0.21 0.18 0.16 0.14 0.12 0.10
Table 4: Updates of mass constraint of beam problem for .

Figure 14 shows the optimized designs using the and phase assignment rules, for decreasing mass constraint from left to right and highlighting ATs in red. It can be seen that the reduction in mass requirements is satisfied by generating new holes rather than decreasing the thickness of the vertical thin walls. A small influence on the topology consistency rule employed is observed on the shape and size of these holes; see onsets of right column of Fig. 14. In this example, all ATs are internal, i.e., the configuration shown in Fig. 5(b). Hence, determine ambiguities as solid or void in the optimized designs can only generate shape changes while keeping the same topology. However, optimized designs can also include BATs resulting from hexahedra with multiple intersections. This is shown in Example 3 using a different objective formulation.

Figure 14: Optimized designs of beam problem using the and rules, and for .

5.2 Structural suspender

In the second problem, we investigate the evolution of a design when the initial topology is significantly affected by the rule chosen to resolve ambiguities. The design domain of size is subject to the boundary conditions specified on the left side of Fig. 15. All four top corners of the domain are clamped in all three directions within the volumes highlighted in Fig. 15. A prescribed traction of is applied in the vertical direction in a rectangular region at the center of the bottom face. Symmetry conditions along the and planes at the center of the domain are considered. The domain is discretized using a mesh with an element edge length of .

Figure 15: Structural suspender problem: (left) design domain with boundary conditions, symmetry planes and dimensions; and (right) hole seeding arrangement in initial design.

The design problem in this example consists of minimizing the displacement at a center region, which is identical to the loaded region, at the bottom of the domain. The optimization problem, which resembles the compliant minimization problem studied in the previous example, reads:


where is defined as:


and represents the average displacement in the direction on a surface , which coincides with the surface on which the traction load is applied. The remainder second and third components of the objective function are explained in Section 4.1. The objective weights are =[0.90, 0.05, 0.05], and a mass constraint with is enforced through a continuation scheme. The continuation step size is 12, and the maximum mass allowed decreases by setting . Both the supports and loading subdomains are excluded from the design domain.

Figure 16: Difference in the initial topology of the structural suspender problem (void phase) as a result of the phase assignment rule used.

The initial seeding for this problem consists of spherical holes placed such that elements of the structured hexahedral mesh have multiple intersections, as shown on the right side of Fig. 15. In this scenario, both internal and external ambiguities are present and connections between the seeded holes may occur. The analysis of this problem focuses on rules and since these rules provide the greatest difference in the resulting geometries. To better visualize this, the additional volume defined as void when using the rule is highlighted in red on the right side of Fig. 16. The void subdomain is shown in this case. It can be seen that rule generates connections between the void spheres. These connections reduce the solid volume by 3%. Also, note that there is approximately a15% difference between the initial average displacements.

Figure 17 shows a 3D view and two plane views of the optimized designs obtained using the and rules on the top and bottom halves, respectively. Both designs converge to four legs spreading from the top corner clamped subdomains to the region at the bottom where the distributed load is applied. Different arrangements of truss-like connections between the legs and the center of the domain on the top are observed. Overall, the optimized structures differ only locally, and similar performances are achieved despite starting from significantly different initial objectives. Hence, in this example, the optimization framework compensates for the phase assignment rule and develop similar discretized geometries.

Figure 17: Optimized designs of suspender problem using the global rules.

The evolution of the ratio is plotted in Fig. 18. Note the drastic reduction in the total volume of ATs at the initial stages of the optimization process. The ATs in the optimized designs represent approximately in the extreme cases examined. Similar to the previous example, the majority of these ambiguities are originated from intersected hexahedra with a single interface. However, unlike in the previous example, here the ATs on the surface produce more pronounced differences in smoothness.

Figure 18: Evolution of ratio between volume of all ATs and total volume of solid phase throughout the structural suspender optimization problem.

Note that in problems with a strict mass requirement such as the one shown here, a large number of holes may need to be seeded relatively close to each other since they have to cover a large portion of the design domain. In such scenarios, rules that favor void phase assignment, i.e. and , can make the initial design unfeasible if the relative difference in volume of ATs is significant. For this reason, using such rules can result in premature merging of holes and convergence to sub-optimal designs.

5.3 Bracket

In the final example, all local and global rules described in this paper are studied with a realistic engineering design problem. This problem is characterized by a non-standard 3D geometry of the design domain, and complexities associated with multiple objective components and a stress constraint. The objective of this optimization problem consists of finding a structure that supports a payload given a set of supports and bolts for attaching the structure to the payload.

The problem design and non-design subdomains are illustrated in Fig. 19. The design subdomain () is colored in light gray, while the non-design subdomains for the payload box (), the supports (), and bolts () are colored in dark gray. A 3D view together with side and bottom views are provided for a better visualization of the constrained subdomains. A uniform pressure load acts on the top surface of the payload box. The entire structure is subject to a body force in direction, representing an equivalent shock loading. The material properties and load conditions for this static problem can be found in Table 5.

Figure 19: Bracket problem domain highlighting: (left) payload design subdomain; and (center, right) supports and bolts non-design subdomains.
Property Value
Young’s Modulus (, solid) x [N/]
Young’s Modulus (, void) x [N/]
Material Density (, solid) x [kg/]
Material Density (, void) [kg/]
Poisson Ratio ( and ) = 0.342 [-]
Pressure load () 1.2x [N/]
Maximum stress () 398.7 [N/]
Table 5: Material properties and load conditions for bracket problem.

The goal of this optimization problem is to minimize the mass, , of the structure that supports the payload box. The optimization problem is formulated as follows:


The first term in the objective is the mass of the structure in the design subdomain. The second term is a strain energy component added to promote a sufficiently stiff structure and prevent an overly aggressive removal of mass early in the design process. The penalty term, , in and enforces constraints on the geometry and is formulated as follows:


with the field defined as:


based on the previously defined subdomains highlighted in Fig. 19. This penalty (i) forces the support structure to stay within the design domain, (ii) prevents contact between the payload box and the support structure except for the bolts, and (iii) avoids altering the LSF in the payload box, attachment bolts and the supports. The remaining penalty terms in the objective ( and ) are defined in Section 4.1. The weights employed in this problem in the objective components and constraints are , , , , , and . The objective scaling parameter is set to