An Interface-enriched Generalized Finite Element Method for Levelset-based Topology Optimization
Sanne J. van den Boom, Jian Zhang, Fred van Keulen, Alejandro M. Aragón
AAn Interface-enriched Generalized Finite Element Method forLevelset-based Topology Optimization
S. J. van den Boom, J. Zhang, F. van Keulen, A. M. Arag´onDepartment of Precision & Microsystems Engineering (PME)Faculty of Mechanical, Maritime & Materials Engineering (3ME)Delft University of Technology (TU Delft)Mekelweg 2, 2628 CD Delft, [email protected] 1, 2020
Abstract
During design optimization, a smooth description of the geometry is important, especially for problemsthat are sensitive to the way interfaces are resolved, e.g. , wave propagation or fluid-structure interaction.A levelset description of the boundary, when combined with an enriched finite element formulation, offersa smoother description of the design than traditional density-based methods. However, existing enrichedmethods have drawbacks, including ill-conditioning and difficulties in prescribing essential boundary con-ditions. In this work we introduce a new enriched topology optimization methodology that overcomes theaforementioned drawbacks; boundaries are resolved accurately by means of the Interface-enriched General-ized Finite Element Method (IGFEM), coupled to a levelset function constructed by radial basis functions.The enriched method used in this new approach to topology optimization has the same level of accuracy inthe analysis as standard the finite element method with matching meshes, but without the need for remesh-ing. We derive the analytical sensitivities and we discuss the behavior of the optimization process in detail.We establish that IGFEM-based levelset topology optimization generates correct topologies for well-knowncompliance minimization problems.
The use of enriched finite element methods in topology optimization approaches is not new; the eX-tended/Generalized Finite Element Method (X/GFEM) [1, 2, 3, 4, 5], for example, has been explored in thiscontext. However, the Interface-enriched Generalized Finite Element Method (IGFEM) has been shown to havemany advantages over X/GFEM [6, 7]. In this work we extend IGFEM to be used in a levelset based topologyoptimization framework.Topology optimization, first introduced by Bendsøe and Kikuchi [8], has been widely used to obtain designsthat are optimized for a certain functionality, e.g. , minimum compliance. In the commonly-used density-basedmethods, a continuous design variable that represents a material density is assigned to each element in thediscretization. The design is pushed towards a black and white design by means of an interpolation function, e.g. , the Solid Isotropic Material with Penalisation (SIMP) [9], that disfavors intermediate density values, alsoreferred to as gray values. A filter is then required to prevent checkerboard-like density patterns, and toimpose a minimum feature size. However, due to the filter, gray values are introduced. Density based topologyoptimization is straightforward to implement and widely available in both research codes and commercialsoftware. However, because the topology is described by a density field on a (usually) structured mesh, materialinterfaces not only contain gray values but also suffer from pixelization or staircasing —staggered boundariesthat follow the finite element mesh. Although a postprocessing step can be performed to smoothen the finaldesign, the analysis during optimization is still based on gray density fields and a staircased representation.This may be detrimental to the approximate solution’s accuracy, especially in cases that are sensitive to theboundary description, such as flow problems [10]. Furthermore, because the location of the material boundaryis not well defined, it is difficult to track the evolving boundary during optimization, for example to imposecontact constraints.The aforementioned drawbacks could be alleviated by the use of geometry-fitted discretization methods,which have been widely used in shape optimization [11]. In these methods, the location of the material interfaceis known throughout the optimization, and the analysis mesh is modified to completely eliminate the pixalization1 a r X i v : . [ c s . C E ] M a r nd gray values. Mesh-morphing methods such as the deformable simplex method [12, 13, 14, 15], anisotropicelements [16], and r -refinement [17], have been demonstrated for topology optimization. Nevertheless, adaptingthe mesh in every design iteration remains a challenge. Not only is it an extra computational step, the changingdiscretization also introduces another complication in the optimization procedure because design variables needto be mapped to the new discretization [18].A more elegant option is to define the material interface independently from the FE discretization, e.g.,implicitly by means of the zero-contour of a levelset function. Although the required mapping between thegeometry and the discretization mesh can be done with an Ersatz method [19], this again introduces gray valuesand staircasing into the analysis. Similarly, NURBS-based topology optimization using the Finite Cell Method(FCM) [20] provides a higher resolution boundary description, that is, however, still staircased. Alternatively,there are methods that allow for a one-to-one mapping of the topology to the analysis mesh, resulting in a non-pixalized boundary description. These methods combine the advantages of clearly defined material interfaceswith the benefits of a fixed discretization mesh used in density-based methods. In the literature, levelsetbased topology optimization has been established based on CutFEM [21] and the eXtended/Generalized FiniteElement Method (X/GFEM) [22, 23, 24].In enriched finite element methods such as X/GFEM, the standard finite element space is augmented byenrichment functions that account for a priori knowledge of the discontinuity of the field or its gradient,due to the material interface. Although X/GFEM has been shown to be advantageous in many applications, e.g. , fluid-structure interaction [25] and fracture mechanics [26], the method also has weaknesses: degrees offreedom (DOFs) corresponding to original mesh nodes do not automatically retain their physical meaning, andessential boundary conditions mostly have to be prescribed weakly . Moreover, the X/GFEM may result in ill-conditioned matrices, in which case a Stable Generalized FEM (SGFEM) is needed [27, 28, 29]. Furthermore,the approximation of stresses can be highly overestimated near material boundaries [30, 31, 32]. Finally, as theenriched functions are associated with original mesh nodes, the accuracy of the approximation may degrade inblending elements—elements that do not have all nodes enriched [33].The Interface-enriched Generalized Finite Element Method (IGFEM) [6] was first introduced as a simplifiedgeneralized FEM to solve problems with weak discontinuities, i.e. , where the gradient field is discontinuous.The method overcomes most problems of X/GFEM: In IGFEM, enriched nodes are placed along interfaces, andenrichment functions are non-zero only in cut elements , i.e. , elements that are intersected by a discontinuity.Furthermore, enrichment functions are exactly zero at original mesh nodes. Therefore, original mesh nodesretain their physical meaning and essential boundary conditions can be enforced directly on non-matchingedges [7, 34, 35]. It was shown that IGFEM is optimally convergent under mesh refinement for problems withoutsingularities [6, 36], and stable by means of scaling enrichment functions or a simple diagonal preconditioner [7].The method has been applied to the modeling of fibre-reinforced composites [36], multi-scale damage evolutionin heterogeneous adhesives [37], microvascular materials with active cooling [6, 36, 38, 39], and the transversefailure of composite laminate [40, 41]. Extensions of IGFEM are found in the Hierarchical Interface-enrichedFinite Element Method (HIFEM) [42], that allows for intersecting discontinuities, and in the Discontinuity-Enriched Finite Element Method (DE-FEM) [35], that provides a unified formulation for both weak and strongdiscontinuities. DE-FEM, which inherits the advantages of IGFEM over X/GFEM, has successfully been appliedto problems in fracture mechanics [35, 43] and fictitious domain or immersed boundary problems with stronglyenforced essential boundary condition [7].In the context of optimization, IGFEM has been explored for NURBS-based shape optimization [44], theshape optimization of microvascular channels [45] and their combined shape and network topology optimiza-tion [46], the optimization of microvascular panels for nanosatellites [47], and optimal cooling of batteries [48].Nevertheless, IGFEM has not yet been used for continuum topology optimization. In this paper we show topol-ogy optimization based on a levelset function, parametrized with Radial Basis Functions (RBFs) [49, 50], incombination with IGFEM. We demonstrate the method on benchmark compliance problems. It should be notedthat no significant performance improvement is expected for these cases, as they are not sensitive to the way theboundaries are discretized, but rather by the bulk. This paper should, however, be seen as the necessary proofof concept before considering more complex cases. The sensitivities are derived and the method is compared todensity-based topology optimization and to the levelset-based Ersatz method. In this work we focus on elastostatics and heat conduction problems on solid domains, as represented in Figure 1.A design domain Ω ⊂ R d is referenced by a Cartesian coordinate system spanned by base vectors { e i } di =1 . Thisdomain is decomposed into a solid material domain and a void domain, denoted by Ω m and Ω v , respectively, suchthat the domain closure is Ω = Ω m ∪ Ω v , and Ω m ∩ Ω v = ∅ . The boundary of the design domain, ∂ Ω ≡ Γ = Ω \ Ω,2 m ne e Γ t ¯ t Γ m Ω v Ω m Γ u x j x k φ < φ > u and Γ t , respectively. The materialdomain is referred to as Ω m , while the void region is denoted Ω v . The insert shows the discretization with amaterial interface, defined by the zero-contour of the levelset function φ , that is non-matching to the mesh.Original mesh nodes and enriched nodes are denoted with and symbols, respectively.is subjected to essential (Dirichlet) boundary conditions on Γ u , and to natural (Neumann) boundary conditionson Γ t , such that Γ = Γ u ∪ Γ t and Γ u ∩ Γ t = ∅ . The material boundary, Γ m = (cid:0) Ω m ∩ Ω v (cid:1) \ Γ, is defined implicitlyby a levelset function, φ ( x ) = 0, that is a function of the spatial coordinate x .For any iteration in the elastostatic optimization procedure, the boundary value problem is solved withprescribed displacements ¯ u : Γ u → R d , prescribed tractions ¯ t : Γ t → R d , and body forces b i ≡ b | Ω i : Ω i → R d , where i = m , v. We denote the field u i as the restriction of u to domain Ω i , i.e. u i ≡ u | Ω i . Note thathere the field is solved on both the material domain and the void domain. However, following the techniquesdescribed in [7], it is also possible to completely remove the void regions from the analysis.We define a vector-valued function space, denoted V ≡ (cid:2) H (Ω) (cid:3) d , where components of v ∈ V belong tothe first-order Sobolev space that satisfies homogeneous essential boundary conditions on Γ u . In this work weonly focus on problems with homogeneous Dirichlet boundary conditions. For problems with non-homogeneousessential boundary conditions, the reader is referred to [35]. The weak form of the elastostatics boundary valueproblem can be written as: Find u ∈ V such that B ( u , v ) = L ( v ) ∀ v ∈ V , (1)where the bilinear and linear forms can be written as B ( u , v ) = (cid:88) i =m , v (cid:90) Ω i σ i ( u i ) : ε i ( v i ) d Ω , (2)and L ( v ) = (cid:88) i =m , v (cid:90) Ω i v i · b i d Ω + (cid:90) Γ t v i · ¯ t d Γ , (3)respectively, where the stress tensor σ i ≡ σ | Ω i follows Hooke’s law for linear elastic materials, σ i = C i (cid:15) i , and C i is the elasticity tensor. Small strain theory is used for the strain tensor, i.e. , ε ( u ) = ( ∇ u + ∇ u (cid:124) ), for theelastostatic boundary value problem.For the heat conduction problem both the trial and the weight function are taken from the space V ≡ H (Ω).For a prescribed temperature u : Γ u → R , prescribed heat flux q : Γ q → R , heat source f i : Ω i → R , andconductivity tensor κ i ≡ κ | Ω i → R d × R d , the bilinear and linear forms for each iteration in the heat complianceproblems equation are written as B ( u, v ) = (cid:88) i =m , v (cid:90) Ω i ∇ v i · ( κ i · ∇ u i ) d Ω (4)and L ( v ) = (cid:88) i =m , v (cid:90) Ω i v i f i d Ω + (cid:90) Γ t v i ¯ q d Γ . (5)3 i α i Figure 2: Schematic representation of enrichment function ψ i corresponding to the enriched DOFs α i , wherethe enriched nodes are shown as . This enrichment function is constructed from the standard Lagrange shapefunction of the integration elements.It is worth noting that interface conditions that satisfy continuity of the field and its tractions (or fluxes) donot appear explicitly in Eq. (1) (or Eq. (4)), because they drop out due to the weight function v (or v ) beingcontinuous along the interface.The design domain is discretized without prior knowledge of the topology as Ω h = (cid:83) i e i , where e i is the i thfinite element, resulting in a mesh that is non-matching to material boundaries. The levelset function, whosezero contour defines the interface between void and material, is then evaluated on the same mesh. This is donefor efficiency, as the mapping needs to be computed only once, and results in discrete nodal levelset values.New enriched nodes are placed at the intersection between element edges and the zero contour of the levelset.The locations of these enriched nodes, denoted x n , are found by linear interpolation between nodes x j and x k of the original mesh: x n = x j − φ j φ k − φ j ( x k − x j ) . (6)Here x j and x k have levelset values of opposite sign. Cut elements are then subdivided into integration elements .Following a Bubnov-Galerkin procedure, the resulting finite dimensional problem is then solved by choosingtrial and weight functions from the same enriched finite element space. The IGFEM approximation can thenbe written as u h ( x ) = (cid:88) i ∈ ι h N i ( x ) U i (cid:124) (cid:123)(cid:122) (cid:125) standard FEM + (cid:88) i ∈ ι w ψ i ( x ) α i (cid:124) (cid:123)(cid:122) (cid:125) enrichment , (7)for elastostatics or u h ( x ) = (cid:88) i ∈ ι h N i ( x ) U i (cid:124) (cid:123)(cid:122) (cid:125) standard FEM + (cid:88) i ∈ ι w ψ i ( x ) α i (cid:124) (cid:123)(cid:122) (cid:125) enrichment , (8)for heat conduction problems, where the first term corresponds to the standard finite element approximation,with shape functions N i ( x ) and corresponding standard degrees of freedom U i , and the second term refers tothe enrichment, characterized by enrichment functions ψ i ( x ) and associated enriched DOFs α i . The index setcontaining all enriched nodes is denoted ι w . Enrichment functions ψ i can be conveniently constructed fromLagrange shape functions of integration elements, as illustrated in Figure 2, while the underlying partition ofunity shape functions are kept intact.Subsequently, the local stiffness matrix k e and force vector f e are obtained numerically; elements that arenot intersected follow the standard FEM procedures. For every integration element k e and f e are defined as k e = (cid:90) e (cid:20) ∆ N ∆ ψ (cid:21) D e (cid:2) ∆ N (cid:124) ∆ ψ (cid:124) (cid:3) d Ω , and f e = (cid:90) e (cid:20) Nψ (cid:21) b e d Ω + (cid:90) δe (cid:20) Nψ (cid:21) ¯ t d Γ , (9)where D e is the constitutive matrix, the shape functions vector N and enrichment functions ψ are stackedtogether. The differential operator ∆ is defined as: ∆ ≡ (cid:34) ∂∂x ∂∂y ∂∂y ∂∂x (cid:35) (cid:62) , ∆ ≡ ∂∂x ∂∂y ∂∂z ∂∂y ∂∂x ∂∂z
00 0 ∂∂z ∂∂y ∂∂x (cid:62) , (10)for elastostatics in 2-D and 3-D, respectively, and ∆ ≡ (cid:104) ∂∂x ∂∂y (cid:105) (cid:62) , ∆ ≡ (cid:104) ∂∂x ∂∂y ∂∂z (cid:105) (cid:62) , (11)4 − x y Figure 3: Compactly supported RBF as given by Equation 13 with the coordinates x i = y i = 0 and the radiusof influence r s = 1.for the heat equation in 2-D and 3-D, respectively.In this work, we are concerned with linear triangular elements, for which a single integration point instandard and integration elements is sufficient. The discrete system of linear equations KU = F is finallyobtained through standard procedures, where K = A e k e , F = A e f e , (12)and A denotes the standard finite element assembly operator.For a more detailed description on IGFEM, the reader is referred to [6]. Although it is possible to directly use the levelset values φ j on original nodes of the finite element mesh asdesign variables, we choose to use compactly supported radial basis functions for the levelset parametrizationfor a number of reasons [50]: i) RBFs give control over the complexity of the designs, and as such, they act similarly to a filter in density-based topology optimization; ii)
By decoupling the finite element analysis mesh from the RBF grid, the design space can be restrictedwithout deteriorating the finite element approximation. This can be used to mitigate approximationerrors due to too coarse discretizations; and iii)
By tuning the radius of support of RBFs, we can ensure that the influence of each design variable extendsover multiple elements. This allows the optimizer to move the boundary further and therefore convergefaster, while using fewer design variables.Figure 3 illustrates a compactly supported RBF θ [49] described by θ i ( r i ) = (1 − r i ) (4 r i + 1) , (13)where the radius r i is defined as r i ( x , x i ) = (cid:112) (cid:107) x − x i (cid:107) r s , (14)and r s is the radius of support. In (14) (cid:107)·(cid:107) denotes the Euclidian norm, and x i the coordinate of the center ofthe RBF θ i .The scalar-valued levelset function φ ( x ) is found as a summation of every non-zero RBF θ i , scaled with itscorresponding design variable s i : φ ( x ) = Θ ( x ) (cid:124) s = (cid:88) i ∈ ι s θ i ( x ) s i , (15)where ι s , the index set containing all design variables, and s ∈ D is a vector of design variables with length | ι s | ,the cardinality of ι s . Finally, evaluating this function at the original nodes of the finite element mesh results inthe levelset vector φ = θ (cid:124) s , (16)where θ is a matrix that needs to be computed only once.5 .3 Optimization The optimization problem is chosen as a minimization of the compliance C with respect to the design variables s that scale the radial basis functions. It needs to be emphasized that compliance minimization is merely ademonstrator problem, and the method is not limited to it. The minimization problem is subject to equilibriumand a volume constraint V c . Furthermore, the design variables are bounded between s min and s max . Thisproblem can be written as s (cid:63) = arg min s ∈ D C = U (cid:124) K ( s ) U subject to KU = F V Ω m ≤ V c s min ≤ s ≤ s max . (17)The Method of Moving Asymptotes (MMA) [51] is employed to solve this optimization problem. The compliance minimization problem is self-adjoint [52], resulting in the sensitivity of the compliance C withrespect to the design variables s as ∂C∂ s = − U (cid:124) ∂ K ∂ s U + 2 U (cid:124) ∂ F ∂ s . (18)Applying the chain rule, the sensitivity of the compliance C with respect to design variable s i can be writtenat the level of integration elements in terms of the nodal levelset values φ j : ∂C∂s i = (cid:88) j ∈ ι i (cid:88) e ∈ ι j (cid:88) n ∈ ι e (cid:18) − u (cid:124) e ∂ k e ∂ x n ∂ x n ∂φ j u e +2 u (cid:124) e ∂ f e ∂ x n ∂ x n ∂φ j (cid:19) ∂φ j ∂s i . (19)In (19), a summation is done over all the nodes in the index set ι i which contains all the original mesh nodesthat are in the support of the RBF corresponding to design variable s i . Then, a summation is done over ι j ,which refers to the index set of all integration elements e in the support of original mesh node j , i.e. , the regionwhere the original shape function N j is nonzero. Lastly, a summation is done over the index set ι e , whichcontains all the enriched nodes n in integration element e . The location of these enriched nodes is denoted x n .Note that a number of terms can be identified in the sensitivity formulation: the derivatives of nodal levelsetvalues with respect to the design variables, ∂φ j /∂s i , the design velocities ∂ x n /∂φ j , and the sensitivity of theelement stiffness matrix and force vector with respect to the location of the n th enriched node, ∂ k e /∂ x n and ∂ f e /∂ x n , respectively.First, the sensitivity of the nodal levelset values with respect to the design variables is simply computed bytaking the derivative of (15) with respect to s as ∂ φ ∂ s = θ (cid:124) . (20)The design velocities ∂ x n /∂φ j also remain straightforward as they are computed by taking the derivative ofEq. (6) as ∂ x n ∂φ j = − φ k ( φ j − φ k ) ( x j − x k ) . (21)Note that the enriched nodes remain on the element edges of the finite element mesh, and thus the direction ofthe design velocity is known a priori .More involved is the sensitivity of the e th integration element stiffness matrix k e with respect to the locationof enriched node n , which for a single integration point can be written as: ∂ k e ∂ x n = ∂j e ∂ x n B (cid:124) e D e B e + j e ∂ B (cid:124) e ∂ x n D e B e + j e B (cid:124) e D e ∂ B e ∂ x n , (22)where B e = (cid:2) ∆ N e J − p ∆ ψ e J − e (cid:3) and J p and J e are the Jacobian of the parent and integration element, respec-tively. Recall that the material within each integration element remains constant, and therefore ∂ D e /∂ x n = .The first term in (22) contains the sensitivity of the Jacobian determinant, and represents the effect of thechanging integration element area; the second and third term contain the sensitivity of the element B e matrix,and represents the effect of the changing shape and enrichment functions. The latter is computed as ∂ B e ∂ x n = (cid:20) ψ e ∂ J − e ∂ x n (cid:21) . (23)6 LE E ¯ t Figure 4: Problem description and initial design for the cantilever beam example in § ∂ ∆ N e /∂ n = ). The Jacobian of the parentelement is not influenced by the enriched node location either ( ∂ J p /∂ x n = ). The enrichment functions havea constant value at the integration point of the integration element ( ∂ ∆ ψ e /∂ n = ). Appendix A describeshow to compute the derivative of the Jacobian inverse and determinant, ∂ J − e /∂ x n and ∂j e /∂ x n , respectively,by straightforward differentiation.Finally, the sensitivity of the design-dependent force vector f e is evaluated. Due to the IGFEM discretization,enriched nodes whose support is subjected to a line or body load contribute to the force vector, implying thatthe derivatives of the force vector are nonzero for cases with line loads or body forces. Similarly to the sensitivityof the element stiffness matrix, the sensitivity of the element force vector consists of two terms, one relating tothe Jacobian derivative, and another containing the function derivatives: ∂ f e ∂ x n = ∂j e ∂ x n [ N (cid:124) e ψ (cid:124) ] b e + j e (cid:20) ∂ N (cid:124) e ∂ x n (cid:124) (cid:21) b e . (24)Here, the element right hand side f e is computed as a function of body sources, but line loads and tractions canbe handled completely analogously. In the second term, only the original shape functions have a contribution.This is because enriched functions always have the same value on the Gauss point of an integration element.However, as the location of the Gauss points with respect to the background element changes, ∂ N e /∂ x n isnonzero, and can be evaluated as ∂ N e ∂ x n = ∆ N e A ∂ x e ∂ x n N e , (25)where A is the inverse isoparametric mapping that maps global coordinates to the local master coordinatesystem of the parent element.Although the sensitivity analysis seems involved, every partial derivative is relatively straightforward andcheap to compute. The enriched method outlined above is demonstrated on a number of classical compliance optimization problems.The results generated by this approach are compared to those generated by open source optimization codes,and the influence of the design discretization is investigated. A 3-D compliance optimization case and a heatsink problem are also considered.In this section, no units are specified; therefore, any consistent unit system can be assumed. For the MMAoptimizer [51], the following settings are used unless otherwise specified: • The lower and upper bounds on the design variables s i are given by s min ,i = − s max ,i = 1 for i = 1 , . . . , I ; • The move limit used by MMA is set to 0 . • A value of 10 is used for the coefficient that controls how aggressively the constraints are enforced.
First, our approach to enriched levelset-based topology optimization is compared to the following open sourcecodes: i) the 99-line SIMP-based code by Sigmund [53]; ii) an 88-line code for parameterized levelset optimization7GFEM SIMP [53] Density mapping [54] Discrete levelset [55]21 × × × × × ×
21, 61 ×
31, 81 ×
41, and 101 ×
51 nodes, respectively.using radial-basis functions and density mapping, proposed by Wei et al. [54]; and iii) a code for discrete levelsettopology optimization with topological derivatives by Challis [55].The optimization problem for this comparison is the widely-used cantilever beam problem, as illustrated inFigure 4. It consists of a 2 L × L rectangular domain that is clamped on the left and subjected to a downwardpoint load ¯ t in the middle of the right side. We set L equal to 1, the volume constraint to 55% of the designdomain volume, and use | ¯ t | = 1. The material domain Ω m is assigned a Young’s modulus E = 1, whereas thevoid domain Ω v has Young’s modulus E = 10 − . Both domains have a Poisson ratio ν = ν = 0 .
3. Note thatit is also possible to give the void regions a stiffness of exactly zero by removing DOFs [7]. However, this wouldentail extra overhead, and to ensure a fair comparison with the other models, in this work it is chosen to use asmall value for the void stiffness.Figure 4 shows the initial design that is used for the IGFEM-based optimization, which is the same as thatused in the paper describing the 88-line code [54]. The other two codes do not require an initial design, as theyare able to nucleate holes. The optimization problem is solved on meshes defined on rectangular grids of 21 × ×
21, 61 ×
31, 81 ×
41, and 101 ×
51 nodes. Our proposed method makes use of triangular meshes, whereasthe other methods use quadrilateral meshes. The RBF mesh used in the IGFEM-based solutions is the same asthe analysis mesh, and a radius of influence of r s = √ · a is used, where a is the distance between two RBFs.The results for each code are illustrated in Figure 5. For all methods, the design becomes more detailedwhen the mesh size is increased. Furthermore, the topologies obtained by each method are roughly the same.It is observed that the resulting designs are similar to those given by the code of Wei et al. , especially for thefiner meshes. Indeed, the our proposed method yields results that have clearly defined (black and white) non-staircased boundaries. Figure 6a shows the convergence behavior of the different codes for the finest mesh. It isobserved that our method leads to the lowest objective function value, which again is similar to that obtainedby Wei et alli ’s code, while initially converging faster in the volume fraction.Figure 6b shows the final compliance as a function of the number of DOFs. Initially, the different methodsall find lower compliance values as the mesh is refined, but the method by Wei et al. and our method findslightly higher values for the finest mesh sizes. This may be explained by the optimizer converging to a localoptimum. For each mesh size, the proposed method finds the lowest compliance value at the cost of addingsome enriched DOFs. 8 20 40 60 8050100150200250 Iteration C o m p li a n ce C IGFEMSIMP [53]Wei et al. [54]Challis [55] 00 . . . . V o l u m e f r a c t i o n V Ω m / V Ω C V Ω m /V Ω (a) C o m p li a n ce C IGFEMSIMP [53]Wei et al. [54]Challis [55] (b)
Figure 6: Results of the cantilever beam problem; (a) shows the compliance and volume ratio as a function of thenumber of iterations, showing the convergence of the different methods considered in the study, (b) illustratesthe final compliance of the cantilever beam as a function of the number of DOFs, optimized using the differentmethods considered in § L LE E ¯ t Figure 7: Problem description and initial design for the MBB beam example in § The influence of the number of radial basis functions is investigated on the well-known MBB beam , which isillustrated in Figure 7. The optimization problem consists of a 3 L × L domain with symmetry conditions on theleft. On the bottom right corner, the domain is simply supported, and a downward force ¯ t is applied on the topleft corner. As in the previous example, the volume constraint is set to 55% of the volume of the total designdomain. The initial design is also indicated in the figure, and the same material properties as in the previousexample are used.This optimization problem is solved on a triangular analysis mesh defined on a grid of size 151 ×
51, usinga discretization of the design space consisting of 61 ×
21, 91 ×
31, 121 ×
41 and 151 ×
51 radial basis functions,so that only for the finest design space discretization, both resolutions match, and an RBF is assigned to everynode in the analysis mesh. The support radius r s is changed together with the design grid so that the overlapof RBFs is the same in each case: r s = √ · a , where a is again the distance between two RBFs.Figure 8 shows the optimized designs. As expected, the level of detail in the design can be controlled bythe RBF discretization. However, it is noted that in the finest RBF mesh, artifacts appear on the designboundary. This behavior will be further analysed in § The original Messerschmitt-B¨olkow-Blohm (MBB) beam problem, as introduced by Olhoff [56], also specified that the upperand lower surfaces have to remain planar, in addition to a maximum allowable deflection and maximum stress. Over the years amore free interpretation of the problem formulation has become commonplace. × × ×
61, 31 ×
91, 41 ×
121 and 51 ×
151 RBFs.behavior is similar in all cases. Moreover, as shown in Figure 9b, the compliance no longer significantly improvesfor the finest RBF discretization.
To show that the method is not restricted to 2-D, a 3-D cantilever beam example is also considered. The materialproperties are the same as those of previous examples. The size of this cantilever beam is 2 L × L × . L , anda structured mesh with 40 × × × × × r s = √ · a . Figure 10 shows the initial design, alongwith the boundary conditions; the right surface is clamped, and a distributed line load with | ¯ t | = 0 . .
001 to preventthe optimizer from moving the boundaries too fast, as only a small number of RBFs is used with a large r s compared to the analysis mesh. The objective function is again the structural compliance, and the volumeconstraint is set to 40% of the total design domain.Figure 11a displays the optimized design; the corresponding convergence plot is shown in Figure 11b, whereit can be seen that the volume satisfied the constraint, and the objective function converges smoothly. Lastly, we consider a heat compliance minimization problem, illustrated in Figure 12. In this two-materialproblem, a highly conductive material ( κ = 1) is distributed within an L × L square domain with a lowerconductivity ( κ = 0 . u = 0, whereas thedomain edges are adiabatic boundaries, i.e. , ¯ q = 0. The entire domain is subjected to uniform heat source f = 1. The problem is solved on a 41 ×
41 node analysis mesh, using 31 ×
31 RBFs with r s = √ · a .As this problem considers a case with a body load, the load vector also contains enriched degrees of freedomthat depend on the locations of the enriched nodes. Therefore, the right hand side is design dependent, i.e. , ∂ F /∂ s (cid:54) = , even though the body load is constant throughout the entire domain.The results of this optimization problem are shown in Figure 13. In the optimized design, narrow featurescan be distinguished that follow the edges of original elements in the background mesh. This is an effect causedby how the intersections are detected, and is investigated in more detail in § § ,
000 Iteration C o m p li a n ce C × × × ×
51 00 . . . . V o l u m e f r a c t i o n V Ω m / V Ω C V Ω m /V Ω (a) . . · C o m p li a n ce C (b) Figure 9: Subfigure (a) shows the convergence of the compliance C and volume fraction V Ω m /V Ω of the MBBbeam using different discretizations of the design space; (b) shows the final compliance of the MBB beam as afunction of the number of design variables. The compliance value no longer decreases at the finest design mesh.¯ t Figure 10: Initial design of the 3-D example with a schematic illustration of the boundary conditions: the rightside is fixed and a vertical downward line load is applied on the bottom-left edge. (a) C o m p li a n ce C C . . . . V o l u m e f r a c t i o n V Ω m / V Ω V Ω m /V Ω (b) Figure 11: Optimized design for the 3-D cantilever beam optimization example (a), and the convergence of thecompliance C and volume fraction V Ω m /V Ω (b). 11 u = 0 L Lκ κ ¯ q = 0¯ q = 0 ¯ q = q = Figure 12: Problem description and initial design for the heat sink. A fixed temperature is applied to thebottom right corner, and a uniform heat source is applied throughout the entire square domain. (a) C o m p li a n ce C C . . . . V o l u m e f r a c t i o n V Ω m / V Ω V Ω m /V Ω (b) Figure 13: Results of the heat sink problem: (a) shows the optimized design of the heat sink problem. Observethat narrow features are created by placing the zero-contour of the levelset near the edges of the original meshelement. In the convergence plot shown in (b), initially some small oscillations are observed in both the objectiveand volume convergence, which can be prevented by the use of a smaller move limit.12 k φ j Figure 14: Illustration of the structure disconnecting due to the levelset discretization. On the left the zero-contour of the levelset, shown in red, defines the design shown in gray. The white arrows indicate the updateof the levelset in the next iteration. On the right, the next iteration is shown, where the narrowest part of thezero-contour lies within a single element, and the nodal levelset values φ j and φ k have the same sign. Therefore,the two intersections shown as are not found, and the structure disconnects, as shown by the new gray design. Oscillations in the objective functions are visible in the convergence of the heat sink problem in Figure 13,and in the coarsest RBF mesh of the MBB beam in Figure 9. As these oscillations might point to inaccuratemodeling or sensitivities, the phenomenon is discussed here in more detail.Recall that intersections between the zero contour of the levelset function and element edges are found usinga linear interpolation of nodal levelset values. Because the levelset function is discretized, no intersections canbe found if two adjacent nodes have the same sign. This effect is illustrated in Figure 14. On the left, the zero-contour of a levelset function is shown in red, which defines a design shown in white/gray. The white arrowsindicate the movement of the material boundary during the next design update. On the right, the updatedlevelset contour is shown in red. As the two adjacent original nodes φ j and φ k now have the same sign, the twointersections between them, shown as cannot be found.The sudden disconnection of the structure due to the levelset discretization is a discontinuous event thatcannot be captured by the sensitivity information. Therefore, as soon as such discontinuous event occurs, thesensitivities and the modeling deviate, and oscillations may occur.This problem can be alleviated by using a smaller move limit, as was done in the 3D MBB example. Anotherapproach that could mitigate this issue is to evaluate the parametrized levelset function on a finer grid, so thatmultiple intersections are found on an element edge. However, the procedure that creates integration elementswould also need to allow for these more complex intersections. It should be noted that neither of these methodscompletely eliminates the problem of discontinuous events. Rather, the methods alleviate the problem bylimiting their chance of occurrence.A related observation can be made in the zigzagged features in the heat sink design of Figure 13. Asillustrated in Figure 15, this pattern occurs when the optimizer tries to make a narrow diagonal feature in theopposite direction of the mesh diagonals. The red intersections cannot be detected, and therefore the structureis disconnected. As a result, the optimizer can only create diagonal narrow features by zigzagging them alongelement edges, as illustrated in Figure 15 on the right. In the final designs of some of the numerical examples, zigzagging of the edges occurred where the zero contourof the levelset function is not perfectly smooth, as detailed in Figure 16. To investigate the cause of this artifact,the test problem of a clamped beam loaded axially shown in Figure 17 was investigated. The compliance wascomputed for a varying zigzagging angle α while keeping the material volume constant.The results in Figure 18 show that the minimum compliance is not found at α = 0, as one would expect,but instead it is found at a negative value of α . Furthermore, the compliance is not symmetric with respect to α = 0 due to the asymmetry of the analysis mesh. The cause of this zigzagging is an approximation error, as themesh is too coarse to accurately describe the deformations and stresses in the structure, similarly to the effectdescribed for nodal design variables in [57]. This effect can be resolved by reducing the design space with respectto the analysis mesh, for example with the use of RBFs, or by increasing the element order. Furthermore, as13igure 15: Illustration of the zigzagged pattern that appears in Figure 13. When a narrow diagonal line isdesired in the opposite direction of the diagonal lines of the mesh, the problem illustrated in Figure 14 resultsin a disconnected line, as shown on the left. Instead, the optimizer will create narrow features along elementedges, as illustrated on the right.Figure 16: Detail of zigzagging that might occur when the design space is not reduced with respect to the FEmesh. 14 Figure 17: Schematic for the zigzagging approximation error. A beam with zigzagging angle α is clamped onthe left, while a concentrated axial loading is applied on the right. The angle α is varied without changing thematerial volume, and the compliance is evaluated. − π π α C Figure 18: The compliance of the test case, illustrated in Figure 17, as a function of the zigzagging angle α .The compliance for this coarse test case is non-symmetric with respect to 0.the non-smoothness is confined to a single layer of background elements, mesh-refinement makes the issue lesspronounced. In this work we introduced a new enriched topology optimization approach based on the Interface-enrichedGeneralized Finite Element Method (IGFEM). The technique yields non-pixelized black and white designs,that do not require any postprocessing. We have derived an analytic expression for the sensitivities, and haveshown that they can be computed with relatively low computational effort. Furthermore, the method wascompared to a number of open source topology optimization codes, based on SIMP, the Ersatz approach, anddiscrete levelsets. The influence of decoupling the design discretization from the analysis mesh was investigatedusing the classical MBB beam optimization problem. A 3-D cantilever beam and a heat sink problem were alsodemonstrated. The convergence behavior was provided for each numerical example. Any numerical artifacts,such as approximation errors and discretization errors of the levelset, as discussed in §
4, can be mitigated bymeans of suitable move limits and radial basis functions, where the latter serve as a sort of filter because theycan control the design complexity.A number of conclusions can be drawn from this work: • The combination of IGFEM with the levelset topology optimization based on RBFs results in crisp bound-aries in both the design representation and the modeling. Because the RBF mesh and analysis mesh arecompletely decoupled, the resolution of the design and the modeling can be chosen independently, as isthe case in any parametrized levelset optimization. In addition, the radial basis functions help in reducingnumerical artifacts, as they act like a black-and-white filter. Lastly, as the RBFs may extend over multipleelements, they allow the boundary to move further and the optimizer to converge faster; • As only one intersection can be detected per element edge, due to the mapping of the levelset to theoriginal mesh nodes, features smaller than a single element might not be described correctly. As discussedin § • Due to approximation error, numerical artifacts may occur which may be exploited by the optimizer whenthe RBF mesh is too fine with respect to the analysis mesh. Another known issue in IGFEM and otherenriched methods, that may be exploited by the optimizer, is the fact that the computation of stressesnear material interfaces may yield inaccurate results [59, 60]; • In this work, we chose to model the void along with the material domain for a number of reasons, includingthe ease of implementation, and the ease of comparing to other methods. However, we could have chosento completely remove the void from the analysis [7], which would reduce computation times and eliminatethe artificial stiffness in the void.The benefits of using an enriched formulation are expected to be more pronounced for problems that relyheavily on an accurate boundary description, such as fluid-structure interaction and wave scattering. In fact,the optimization of the latter is the subject of an incoming publication. Furthermore, by having an enrichedformulation, the technique has the potential to easily handle problems where boundary conditions have to beprescribed on evolving boundaries. In particular, and unlike other immersed methods, with IGFEM it is evenpossible to prescribe strongly non-homogeneous essential (Dirichlet) boundary conditions [35, 7].
Acknowledgements
The authors would like to thank Krister Svanberg for providing us with the MMA implementation.
References [1] J.T. Oden, C.A.M. Duarte, and O.C. Zienkiewicz. A new cloud-based hp finite element method.
ComputerMethods in Applied Mechanics and Engineering , 153(1):117 – 126, 1998.[2] Nicolas Mo¨es, John Dolbow, and Ted Belytschko. A finite element method for crack growth withoutremeshing.
International Journal for Numerical Methods in Engineering , 46(1):131–150, 1999.[3] N. Mo¨es, M. Cloirec, P. Cartraud, and J.-F. Remacle. A computational approach to handle complexmicrostructure geometries.
Computer Methods in Applied Mechanics and Engineering , 192(28):3163 –3177, 2003.[4] Ted Belytschko, Robert Gracie, and Giulio Ventura. A review of extended/generalized finite element meth-ods for material modeling.
Modelling and Simulation in Materials Science and Engineering , 17(4):043001,2009.[5] Alejandro M. Arag´on, C. Armando Duarte, and Philippe H. Geubelle. Generalized finite element enrichmentfunctions for discontinuous gradient fields.
International Journal for Numerical Methods in Engineering ,82(2):242–268, 2010.[6] S. Soghrati, A.M. Arag´on, C. Armando Duarte, and P.H. Geubelle. An interface-enriched generalizedfem for problems with discontinuous gradient fields.
International Journal for Numerical Methods inEngineering , 89(8):991–1008, 2012.[7] Sanne J. van den Boom, Jian Zhang, Fred van Keulen, and Alejandro M. Arag´on. A stable interface-enrichedformulation for immersed domains with strong enforcement of essential boundary conditions.
InternationalJournal for Numerical Methods in Engineering , 0(0).[8] Martin Philip Bendsøe and Noboru Kikuchi. Generating optimal topologies in structural design using ahomogenization method.
Computer Methods in Applied Mechanics and Engineering , 71(2):197 – 224, 1988.[9] M. P. Bendsøe. Optimal shape design as a material distribution problem.
Structural optimization , 1(4):193–202, Dec 1989.[10] C. H. Villanueva and K. Maute. CutFEM topology optimization of 3D laminar incompressible flow prob-lems.
Comput. Methods Appl. Mech. Eng. , 320:444–473, 2017.[11] Matthew L. Staten, Steven J. Owen, Suzanne M. Shontz, Andrew G. Salinger, and Todd S. Coffey. Acomparison of mesh morphing methods for 3d shape optimization. In William Roshan Quadros, edi-tor,
Proceedings of the 20th International Meshing Roundtable , pages 293–311, Berlin, Heidelberg, 2012.Springer Berlin Heidelberg. 1612] Marek Krzysztof Misztal and Jakob Andreas Bundefinedrentzen. Topology-adaptive interface trackingusing the deformable simplicial complex.
ACM Trans. Graph. , 31(3), June 2012.[13] Asger Nyman Christiansen, Morten Nobel-Jørgensen, Niels Aage, Ole Sigmund, and Jakob AndreasBærentzen. Topology optimization using an explicit interface representation.
Structural and Multidis-ciplinary Optimization , 49(3):387–399, Mar 2014.[14] Asger Nyman Christiansen, J. Andreas Bærentzen, Morten Nobel-Jørgensen, Niels Aage, and Ole Sigmund.Combined shape and topology optimization of 3d structures.
Computers & Graphics , 46:25 – 35, 2015.[15] Mingdong Zhou, Haojie Lian, Ole Sigmund, and Niels Aage. Shape morphing and topology optimizationof fluid channels by explicit boundary tracking.
International Journal for Numerical Methods in Fluids ,88(6):296–313, 2018.[16] K.E. Jensen. Anisotropic mesh adaptation and topology optimization in three dimensions.
Journal ofMechanical Design, Transactions of the ASME , 138(6), 2016.[17] S. Yamasaki, S. Yamanaka, and K. Fujita. Three-dimensional grayscale-free topology optimization us-ing a level-set based r-refinement method.
International Journal for Numerical Methods in Engineering ,112(10):1402–1438, 2017.[18] N. P. van Dijk, K. Maute, M. Langelaar, and F. van Keulen. Level-set methods for structural topologyoptimization: a review.
Structural and Multidisciplinary Optimization , 48(3):437–472, Sep 2013.[19] G. Allaire, C. Dapogny, and P. Frey. Shape optimization with a level set based mesh evolution method.
Computer Methods in Applied Mechanics and Engineering , 282:22 – 53, 2014.[20] Yuliang Gao, Yujie Guo, and Shijie Zheng. A nurbs-based finite cell method for structural topologyoptimization under geometric constraints.
Computer Aided Geometric Design , 72:1 – 18, 2019.[21] Carlos H. Villanueva and Kurt Maute. Cutfem topology optimization of 3d laminar incompressible flowproblems.
Computer Methods in Applied Mechanics and Engineering , 320:444 – 473, 2017.[22] T. Belytschko, S. P. Xiao, and C. Parimi. Topology optimization with implicit functions and regularization.
International Journal for Numerical Methods in Engineering , 57(8):1177–1196, 2003.[23] Carlos H. Villanueva and Kurt Maute. Density and level set-xfem schemes for topology optimization of3-d structures.
Computational Mechanics , 54(1):133–150, Jul 2014.[24] Pai Liu, Yangjun Luo, and Zhan Kang. Multi-material topology optimization considering interface behaviorvia xfem and level set method.
Computer Methods in Applied Mechanics and Engineering , 308:113 – 133,2016.[25] U. M. Mayer, A. Popp, A. Gerstenberger, and Wolfgang A. Wall. 3d fluid–structure-contact interactionbased on a combined xfem fsi and dual mortar contact approach.
Computational Mechanics , 46(1):53–67,Jun 2010.[26] Thomas-Peter Fries and Ted Belytschko. The extended/generalized finite element method: An overview ofthe method and its applications.
International Journal for Numerical Methods in Engineering , 84(3):253–304, 2010.[27] I. Babuˇska and U. Banerjee. Stable generalized finite element method (sgfem).
Computer Methods inApplied Mechanics and Engineering , 201-204:91 – 111, 2012.[28] V. Gupta, C.A. Duarte, I. Babuˇska, and U. Banerjee. A stable and optimally convergent generalized fem(sgfem) for linear elastic fracture mechanics.
Computer Methods in Applied Mechanics and Engineering ,266:23 – 39, 2013.[29] K. Kergrene, I. Babuˇska, and U. Banerjee. Stable generalized finite element method and associated iterativeschemes; application to interface problems.
Computer Methods in Applied Mechanics and Engineering , 305:1– 36, 2016.[30] Laurent Van Miegroet and Pierre Duysinx. Stress concentration minimization of 2d filets using x-fem andlevel set description.
Structural and Multidisciplinary Optimization , 33(4):425–438, Apr 2007.[31] Lise No¨el and Pierre Duysinx. Shape optimization of microstructural designs subject to local stress con-straints within an xfem-level set framework.
Structural and Multidisciplinary Optimization , 55(6):2323–2338, Jun 2017. 1732] Ashesh Sharma and Kurt Maute. Stress-based topology optimization using spatial gradient stabilized xfem.
Structural and Multidisciplinary Optimization , 57(1):17–38, Jan 2018.[33] Thomas-Peter Fries. A corrected xfem approximation without problems in blending elements.
InternationalJournal for Numerical Methods in Engineering , 75(5):503–532, 2008.[34] A.C. Ramos, A.M. Arag´on, S. Soghrati, P.H. Geubelle, and J.-F. Molinari. A new formulation for imposingdirichlet boundary conditions on non-matching meshes.
International Journal for Numerical Methods inEngineering , 103(6):430–444, 2015.[35] A.M. Arag´on and A. Simone. The discontinuity-enriched finite element method.
International Journal forNumerical Methods in Engineering , 112(11):1589–1613, 2017.[36] Soheil Soghrati and Philippe H. Geubelle. A 3d interface-enriched generalized finite element method forweakly discontinuous problems with complex internal geometries.
Computer Methods in Applied Mechanicsand Engineering , 217-220:46 – 57, 2012.[37] A.M. Arag´on, S. Soghrati, and P.H. Geubelle. Effect of in-plane deformation on the cohesive failure ofheterogeneous adhesives.
Journal of the Mechanics and Physics of Solids , 61(7):1600–1611, 2013.[38] S. Soghrati, P.R. Thakre, S.R. White, N.R. Sottos, and P.H. Geubelle. Computational modeling anddesign of actively-cooled microvascular materials.
International Journal of Heat and Mass Transfer , 55(19-20):5309–5321, 2012.[39] S. Soghrati, A.R. Najafi, J.H. Lin, K.M. Hughes, S.R. White, N.R. Sottos, and P.H. Geubelle. Computa-tional analysis of actively-cooled 3d woven microvascular composites using a stabilized interface-enrichedgeneralized finite element method.
International Journal of Heat and Mass Transfer , 65:153–164, 2013.[40] Xiang Zhang, David R. Brandyberry, and Philippe H. Geubelle. Igfem-based shape sensitivity analysis ofthe transverse failure of a composite laminate.
Computational Mechanics , May 2019.[41] Maryam Shakiba, David R. Brandyberry, Scott Zacek, and Philippe H. Geubelle. Transverse failure ofcarbon fiber composites: Analytical sensitivity to the distribution of fiber/matrix interface properties.
International Journal for Numerical Methods in Engineering , 0(0).[42] S. Soghrati. Hierarchical interface-enriched finite element method: An automated technique for mesh-independent simulations.
Journal of Computational Physics , 275:41 – 52, 2014.[43] Jian Zhang, Sanne J. van den Boom, Fred van Keulen, and Alejandro M. Arag´on. A stable discontinuity-enriched finite element method for 3-d problems containing weak and strong discontinuities.
ComputerMethods in Applied Mechanics and Engineering , 355:1097 – 1123, 2019.[44] Ahmad R. Najafi, Masoud Safdari, Daniel A. Tortorelli, and Philippe H. Geubelle. Shape optimizationusing a nurbs-based interface-enriched generalized fem.
International Journal for Numerical Methods inEngineering , 111(10):927–954, 2017.[45] Marcus Hwai Yik Tan and Philippe H. Geubelle. 3d dimensionally reduced modeling and gradient-basedoptimization of microchannel cooling networks.
Computer Methods in Applied Mechanics and Engineering ,323:230 – 249, 2017.[46] Reza Pejman, Sherif H. Aboubakr, William H. Martin, Urmi Devi, Marcus Hwai Yik Tan, Jason F. Patrick,and Ahmad R. Najafi. Gradient-based hybrid topology/shape optimization of bioinspired microvascularcomposites.
International Journal of Heat and Mass Transfer , 144:118606, 2019.[47] Marcus Hwai Yik Tan, Devin Bunce, Alexander R. M. Ghosh, and Philippe H. Geubelle. Computationaldesign of microvascular radiative cooling pasonels for nanosatellites.
Journal of Thermophysics and HeatTransfer , 32(3):605–616, 2018.[48] Marcus Hwai Yik Tan, Ahmad R. Najafi, Stephen J. Pety, Scott R. White, and Philippe H. Geubelle. Multi-objective design of microvascular panels for battery cooling applications.
Applied Thermal Engineering ,135:145 – 157, 2018.[49] Holger Wendland. Piecewise polynomial, positive definite and compactly supported radial functions ofminimal degree.
Advances in Computational Mathematics , 4(1):389–396, Dec 1995.[50] Shengyin Wang and Michael Yu Wang. Radial basis functions and level set method for structural topologyoptimization.
International Journal for Numerical Methods in Engineering , 65(12):2060–2090, 2006.1851] Krister Svanberg. The method of moving asymptotes—a new method for structural optimization.
Inter-national Journal for Numerical Methods in Engineering , 24(2):359–373, 1987.[52] Martin Bendsøe and Ole Sigmund.
Topology optimization. Theory, methods, and applications. 2nd ed.,corrected printing . 01 2004.[53] O. Sigmund. A 99 line topology optimization code written in matlab.
Structural and MultidisciplinaryOptimization , 21(2):120–127, Apr 2001.[54] Peng Wei, Zuyu Li, Xueping Li, and Michael Yu Wang. An 88-line matlab code for the parameterizedlevel set method based topology optimization using radial basis functions.
Structural and MultidisciplinaryOptimization , 58(2):831–849, Aug 2018.[55] Vivien J. Challis. A discrete level-set topology optimization code written in matlab.
Structural andMultidisciplinary Optimization , 41(3):453–464, Apr 2010.[56] Niels Olhoff, Martin P. Bendsøe, and John Rasmussen. On cad-integrated structural topology and designoptimization.
Computer Methods in Applied Mechanics and Engineering , 89(1):259 – 279, 1991.[57] V. Braibant and C. Fleury. Shape optimal design using b-splines.
Computer Methods in Applied Mechanicsand Engineering , 44(3):247 – 267, 1984.[58] Alejandro M. Arag´on, Bowen Liang, Hossein Ahmadian, and Soheil Soghrati. On the stability and in-terpolating properties of the hierarchical interface-enriched finite element method.
Computer Methods inApplied Mechanics and Engineering , page 112671, 2020.[59] Soheil Soghrati, Anand Nagarajan, and Bowen Liang. Conforming to interface structured adaptive meshrefinement: New technique for the automated modeling of materials with complex microstructures.
FiniteElements in Analysis and Design , 125:24 – 40, 2017.[60] Anand Nagarajan and Soheil Soghrati. Conforming to interface structured adaptive mesh refinement: 3dalgorithm and implementation.
Computational Mechanics , 62(5):1213–1238, Nov 2018.[61] Jan R. Magnus and Heinz Neudecker.
Matrix Differential Calculus with Applications in Statistics andEconometrics . 01 2007.
A Derivatives of the Jacobian inverse and determinant
In the sensitivity computation discussed in § j e , the derivative can thus be computed as: ∂j e ∂ x n = tr (cid:18) adj ( J e ) ∂ J e ∂ x n (cid:19) , (26)The sensitivity of the Jacobian inverse can be computed by realizing that J e J − e = I : ∂ J e J − e ∂ x n = ∂ J e ∂ x n J − e + J e ∂ J − e ∂ x n = ∂ I ∂ x n = 0 , (27)and solving for ∂ J − e /∂ x n : ∂ J − e ∂ x n = − J − e ∂ J e ∂ x n J − e . (28)For both (26) and (28), the sensitivity of the Jacobian is required; as the Jacobian of the integration element iscomputed as J e = x e ∆ ψ e it can be computed as ∂ J e ∂ x n = ∂ x e ∂ x n ∆ ψ e , (29)where ∂ x e /∂ x n is simply a selection matrix consisting of zeros except for a one on the coordinates of interestfor enriched node nn