Topology optimization is a powerful computational tool to design optimal structures under given loads and boundary conditions. Since the seminal work of Bendsøe and Kikuchi (1988), the field of topology optimization has experienced tremendous growth and had a major impact on several areas of engineering and technology. For an overview of this field, we refer the interested readers to textbooks (Christensen and Klarbring 2009;Haftka and Gürdal 2012;Bendsoe and Sigmund 2013) and review paper (Rozvany 2009). Among various topology optimization approaches, the density-based approach is commonly adopted, in which the geometry is parametrized by a material density function and the displacement field is approximated by finite elements. Because of its simplicity and efficiency, the choice of piece-wise constant density parametrization, where each finite element is assigned with a constant density, is typically employed in conjunction with a lower-order Lagrangian-type displacement approximation. However,this choice of density parametrization suffers from numerical instabilities such as checkerboard patterns and one-node connections (Diaz and Sigmund 1995;Sigmund et al. 1998). In a general setting, Jog and Haber (1996) formulated the topology optimization as a mixed variational problem and demonstrated how different choices of displacement and density interpolations affect the appearance of numerical instabilities. Along this direction, different choices of displacement and density approximations were investigated. One popular choice is to approximate material density fields with nodal design variables. For instance, Matsui and Terada (2004) and Rahmatalla and Swan (2004) employed continuous approximations of the density fields interpolated by Lagrangian finite element basis functions. This approach, sometimes known as continuous approximation of material distribution (CAMD), is effective in preventing the checkerboard patterns but may lead to other forms of numerical instabilities such as "islanding" and "layering" with linear finite elements (Rahmatalla and Swan 2004). Guest et al. (2004) introduced nodal design variables and employed a projection map with an embedded length scale to construct an element-wise constant density field. Belytschko et al. (2003) proposed a formulation in which implicit functions are used to describe topologies of the designs. The implicit functions are described by their nodal values and interpolated by C 0 finite element shape functions. Alternatively, Kang and Wang (2011) proposed to decouple the approximations of density and displacement fields into two independent discretizations. In their approach, the density field is interpolated from nodal design variables using non-local Shepard interpolants, and the displacement field is approximated using standard finite elements. This feature helps the method to better handle problems with complex design domains and allows for flexible adaptive representation of the structural topologies (Wang et al. 2013). Other related choices also include the use of higherorder displacement approximations (Diaz and Sigmund 1995;Matsui and Terada 2004;Groen et al. 2017), nonconforming elements (Jang et al. 2003;Jang et al. 2005), and approximating density and displacement fields on different discretizations (Paulino and Le 2009;Guest and Smith 2010). The multi-resolution topology optimization approaches (see, e.g., Nguyen et al. (2010Nguyen et al. ( , 2012)); Filipov et al. 2016;Groen et al. 2017) belong to the last family. In addition to various types of numerical instabilities, topology optimization is typically performed on structured meshes (e.g., uniform grids), which may lead to mesh-dependent solutions (Antonietti et al. 2017) and limited ability to discretize complex design domains.
In recent studies, polygonal finite elements have been shown to be effective in suppressing checkerboard patterns and reducing mesh dependency in the solutions of topology optimization (Talischi et al. 2009(Talischi et al. , 2010(Talischi et al. 2012a;;Gain et al. 2015;Antonietti et al. 2017). Moreover, as compared to standard FEM uniform grids, polygonal elements are more versatile in discretizing complex domains. To this effect, a robust mesh generator named PolyMesher, able to discretize arbitary 2D domains with polygonal elements, has been developed (Talischi et al. 2012b). Other efforts in developing polygonal and polyhedral meshers also include Abdelkader et al. (2018) and Pouderoux et al. (2017). In addition, the geometrical flexibility of polygonal finite elements also makes them attractive for mesh adaptations in topology optimization (see, e.g., Nguyen-Xuan (2017) and Hoshina et al. (2018)). However, most of the abovementioned investigations are in 2D and efficient extensions to 3D problems possess several challenges.
The first challenge comes from the difficulties of polyhedral FEM (Hormann and Sukumar 2018). For polyhedral finite elements, one major difficulty is associated with obtaining the shape functions and their gradients. Although several shape functions exist in the literature with closedform expressions, most of them are limited to certain classes of element geometry. For example, the Wachspress shape functions only work with strictly convex and simple polyhedra (meaning the collection of faces that include each vertex consists of exactly three faces) (Floater et al. 2014), and the mean value coordinates are mainly applicable to polyhedra with simplicial faces (Floater et al. 2005). Those limitations in the element geometry could potentially affect the accuracy of the polyhedral FEM when dealing with degenerated elements, such as the ones with skinny faces or small edges. Other types of shape functions, such as harmonic (Martin et al. 2008;Bishop 2014) and max-entropy (Arroyo and Ortiz 2006;Hormann and Sukumar 2008), allow for more general polyhedra (e.g., concave ones). However, their values and gradients at integration points can only be computed numerically, which could be undesirable from a computational perspective for large-scale problems. Another difficulty of the polyhedral FEM is associated with numerical integration (Talischi and Paulino 2014;Manzini et al. 2014;Bishop 2014). Because there is no iso-parametric mapping for polyhedral finite elements, numerical integration needs to be performed in the physical domain. Due to the nonpolynomial nature of their shape functions, efficient yet consistent numerical integration rules are difficult to construct on general polyhedral finite elements. Ensuring convergence of the numerical solution typically requires a prohibitively large number of integration points in each element (see Talischi et al. (2015) and Chi et al. (2016Chi et al. ( , 2015) ) for some recent works that attempt to overcome this difficulty).
In this work, we also identify another major challenge to use polyhedral meshes in 3D topology optimization, which is related to the computational efficiency. As mentioned in the preceding paragraphs, typically, topology optimization adopts an element-wise constant density approximation and a lower-order displacement approximation with the degrees of freedom (DOFs) located at the vertices of the mesh. As a polyhedral mesh (e.g., centroidal Voronoi tessellation (CVT)) usually contains significantly more vertices than elements, this typical choice of density and displacement approximations becomes considerably less computationally economical on polyhedral meshes. In particular, when compared with uniform grids of a similar number of elements (thus similar numbers of densities and design variables), we need to solve a much larger structural system in each optimization step if a polyhedral mesh is considered.
The virtual element method (VEM) (Beirão da Veiga et al. 2013) is a recently proposed approach that has the potential to overcome the difficulties of the polyhedral FEM. The VEM can be viewed as a generalization of the FEM and is able to effectively handle arbitrary polygonal and polyhedral meshes. One main feature of the VEM is that its shape functions are defined implicitly through a suitable set of partial differential equations (PDEs). Instead of solving the PDEs for the values of shape functions and their gradients at the integration points, the VEM constructs a set of projection operators which project the shape functions and their gradients onto polynomial functions of suitable orders (Beirão da Veiga et al. 2013, 2014). By construction, these projections can be exactly computed using only the DOFs of the unknown fields. Another major feature of the VEM is that it decomposes the weak forms into consistency and stability terms, both of which can be directly formed using the projections of the shape functions and their gradients as well as the DOFs of the unknown fields (Beirão da Veiga et al. 2013). As a result, for any element geometry, the VEM only needs to integrate polynomials (and not non-polynomials as in the polyhedral FEM). For lower-order VEM, which is the focus of this paper, no numerical integration is needed. Because of such attractive features, the VEM has gained significant visibility in the computational mechanics community. For instance, the VEM has been developed for linear elasticity (Beirão da Veiga et al. 2013;Gain et al. 2014;Artioli et al. 2017), small deformation non-linear elastic and inelastic (Beirão da Veiga et al. 2015;Artioli et al. 2017;Taylor and Artioli 2018), finite deformation elasticity and elastoplasticity (Chi et al. 2017;Wriggers et al. 2017Wriggers et al. , 2018;;Wriggers and Hudobivnik 2017), plate bending (Brezzi and Marini 2013;Antonietti et al. 2018;Mora et al. 2018;Zhao et al. 2016), and damage and fracture problems (De Bellis et al. 2018;Benedetto et al. 2018), to name a few. We also mention that the VEM has been adopted to solve the state equations in topology optimization on unstructured polygonal (Antonietti et al. 2017) and polyhedral (Gain et al. 2015) meshes as well as on Escher-based tessellations (Paulino and Gain 2015).
In this work, we propose a new VEM-based topology optimization framework on general polyhedral discretizations.
To address the abovementioned challenges of efficiently formulating topology optimization on polyhedral meshes, we adopt the VEM in both structural and optimization problems. Similarly to Gain et al. (2015), we use VEM to solve the elasticity equation in the structural problem. The capability of VEM in handling arbitrary element geometry allows us to solve the structural problem more efficiently (i.e., no numerical integration is needed) and robustly (i.e., with respect to degenerated elements) on polyhedral meshes. Exploiting the flexibility of VEM in defining local spaces, this work also introduces an enriched continuous approximation of material densities using nodal VEM basis functions. As compared to the standard element-wise constant density approximation, the continuous approximation contains a greater number of DOFs for any given polyhedral mesh and can thus improve the quality of structural topology parameterizations. Moreover, this work explores various approaches of discretizing complex domains in 3D, including regular polyhedra-dominated and unstructured polyhedral meshes, and investigates their influences on the quality and the numerical stability of solutions in the topology optimization. Several design examples are presented on non-Cartesian domains to demonstrate the main features of the VEM-based topology optimization framework. To complement the library of educational codes (e.g., Sigmund (2001), Andreassen et al. (2011), Talischi et al. (2012a), Liu and Tovar (2014), Pereira et al. (2016), Wei et al. (2018), andSutton (2017)) in the topology optimization and VEM literature, the source code for a MATLAB implementation, named PolyTop3D, is provided as (electronic) Supplementary Material.
The remainder of this paper is organized as follows. Section 2 provides an overview of the theory and implementation of VEM for 3D linear elasticity problems. Section 3 introduces the proposed VEM-based topology optimization together with a simple numerical example comparing the performance of the proposed continuous density parametrization with the standard element-wise constant one. In Section 4, we present a set of design examples featuring non-Cartesian domains to highlight the main features of the proposed VEM-based topology optimization framework. Section 5 contains several concluding remarks and future research directions. In Appendix, the implementation of the PolyTop3D is presented and the computational efficiency of the code is demonstrated.
We consider an elastic solid ∈ R 3 with its boundary denoted by ∂ . The solid is subjected to a prescribed displacement u 0 on one portion of the boundary u and a traction t on the other portion t , such that u ∪ t = ∂ and u ∩ t = ∅. In the interior of , the solid is subjected to a body force f. For a given displacement field u, the linearized strain tensor ε(u) is obtained as ε(u) = 1/2[∇u+ (∇u) T ], where ∇ stands for the gradient operator. The stress tensor is given by the constitutive relationship:
where C is the elasticity tensor.
The weak form of the elasticity problem consists of finding the displacement u among the space K of kinematically admissible displacements:
where,
and K 0 stands for the space of kinematically admissible displacements that vanish on u .
For a general polyhedral element E consisting of planar faces, this subsection describes the construction of the local virtual space V(E) following the technique introduced in Beirão da Veiga et al. ( 2018) and Ahmad et al. (2013). In this technique, the construction of the virtual space for E depends on the virtual spaces on the faces of E. Therefore, before introducing the construction of the virtual space on E, we will first describe the construction of virtual spaces on faces of E. We also note that, in the literature, there exist different approaches to construct local virtual spaces on general polyhedral elements. The interested readers are referred to Ahmad et al. (2013), Beirão da Veiga et al. (2017), Gain et al. (2014), andChi et al. (2017) for further information.
Here, we introduce the definition of the virtual space V(f ) on a generic face f of E, which is assumed to be a planar polygon. The basic idea is to first introduce a preliminary space denoted by V(f ) as:
where e denotes a generic edge of face f and P k (•) is the polynomial space of order k. For the preliminary space, V(f ), a set of DOFs consists of Beirão da Veiga et al. (2013) and Ahmad et al. (2013):
• The values of v at vertices of f, (5)
• The mean value of v over f, i.e.,
where |f | stands for the area of face f . By identity:
we can equivalently replace the DOF (6) by the following integral (Beirão da Veiga et al. 2017):
where
being the centroid of face f . In fact, once we know the above integral, we can compute the mean value of v over f using the DOFs (5) and the identity (7).
Having defined the preliminary virtual space V(f ), we can define the formal virtual space V(f ) ⊂ V(f ) on face f such that:
By definition, we can show that P 1 ⊆ V(f ) and ( 5) constitutes a complete set of DOFs of V(f ). Using this set of DOFs, we can exactly compute moment of v on f according to the identity (7) as:
where n e denotes the outward normal vector of edge e.
Noticing that
• n e takes constant value for any points on edge e (which is assumed to be straight), we can simply evaluate it at any point a e on e, i.e., x f • n e = (a ex f c ) • n e , ∀a e ∈ e.
Once we know the virtual space on each face f , we are ready to define the virtual element space V(E) on E. Following the same concept, we define the final form of the virtual space V(E) as:
where x E . = xx E c with x E c being the centroid of E. Similarly to V(f ), we can define the set of DOFs of V(E) as the values of its functions on the vertices of E. Because V(E) includes the polynomial space P 1 (E), i.e., P 1 (E) ⊆ V(E), we can define a projection operator ∇ E : V(E) → P 1 (E) such that for any v ∈ V(E):
where we use x v to denote the position vector of a generic vertex in E. Because both ∇ ∇ E v and ∇p 1 are constant vectors, the first condition in (12) can be simplified as:
where |E| denotes the volume of element E and n f is the (constant) outward unit normal vector of face f . This condition ensures that the gradient of the projection ∇ E v equals the average gradient of v over E. Recalling from the definition of V(f ) that, given the value of v at the vertices of f , we can compute the moment of v over f (see ( 10)), and consequently, we can explicitly compute the boundary integral on the right-hand side of (13) using the DOFs of v and geometric information of E. On the other hand, the second condition in (12) determines the constant component of the projection by ensuring that, when evaluated at the vertices of E, the average value of ∇ E v is equal to the average value of v. Furthermore, making use of the following identity:
and the definition of V(E), we can express the moment of v over E as:
• n f is constant for any points x on planar f , we can simply evaluate this quantity at any point a f on face f , i.e., x E = a fx E c •n f , ∀a f ∈ f . Thus, we can compute the moment of v over E exactly as:
using only the DOFs of v (recalling (10)).
Consider a polyhedron E consisting of m vertices numbered as x 1 , ..., x m and denote F i as the set of faces that are connected to the ith vertex. Suppose that a face f
.., m f , we locally renumber those vertices in a counterclockwise fashion with respect to the outward normal n f which points out of the element. We also utilize a map G f to denote the relation between the global numbering and local numbering of the vertices on face f . If the ith vertex of E (under global numbering) is numbered as the j th vertex of f (under local numbering), we write j = G f (i). Figure 1 shows an illustration of the above notation for a hexahedral element. Following this notation, this subsection provides the implementational details to constructthe virtual space V(E). In particular, we focus on the procedures of calculating the projection ∇ E ϕ i and the moment E ϕ i , where ϕ i is the basis function of V(E) associated with the ith vertex. For vertex i, ϕ i is defined to be a function which belongs to V(E) and takes 1 at the ith vertex and 0 at other vertices of E.
Based on the definition of projection operator ∇ E in (12), we can express ∇ E ϕ i as:
where projection, we first need to know the moment of ϕ i on f, ∀f ∈ F i . Using relation (10) and realizing that ϕ i varies linear on ∂f , we can compute the moment of ϕ i on f as:
where ∧ stands for the cross product and the convention is used:
Substituting the above relation into ( 13), we can show that:
and, by ( 17), we arrive at the following expression for
Moreover, since ∇ E ϕ i ∈ P 1 (E) ⊆ V(E), the set of basis function ϕ i can interpolate any linear function exactly, namely, 1 = m i=1 ϕ i and x = m i=1 ϕ i x i . Therefore, we can rewrite the relation (20) as a linear combination of the set of basis functions ϕ 1 , ..., ϕ m as:
The above expression can be further simplified in matrix form as:
where P E (ij ) is the (i, j )th component of the matrix P E ∈ R m×m of the form:
In terms of implementation, we form P E using matrix multiplication as follows. We first define G E ∈ R m×3 as a matrix collecting the information of ∇ ∇ E ϕ i , i = 1, ..., m:
and compute P E as:
where 1 ∈ R m×1 is a column vector with all components being 1.
For later use, we also provide the expression to compute the moment of ϕ i over E, which, according to expressions ( 16) and ( 18), takes the form
(26)
We consider a discretization, denoted by h , of the solid into non-overlapping polyhedra consisting of planar faces, where h denotes the average element size. We denote t h and u h as the portions of the mesh boundary where the traction and displacement boundary conditions are applied, respectively. We define the discrete global displacement space K h ⊂ K as:
T belongs to the vectorial space [V(E)] 3 , having three displacement DOFs for each vertex of E. In the following discussion, we define the projection operator for vector fields ∇ E : [V(E)] 3 → [P 1 (E)] 3 as the action of ∇ E on every component of the vector field, e.g.,
Having introduced the discretization, we can decompose the continuous bilinear form (i.e., bilinear form with exact integration) a(u h , v h ) into the summation of element-level contributions:
By exploiting the property that projection ∇ E v satisfies the following orthogonality condition:
we can decompose a E (u h , v h ) as:
Because both of its arguments are linear functions, we can evaluate the first term exactly as:
The second term, on the other hand, involves higher-order displacement components and is typically impossible to evaluate exactly. In the VEM framework, this term is approximated by:
where S E (•, •) is a bilinear form that is inexpensive to compute and satisfies the coercivity condition; and α E is a scaling parameter that scales S E (•, •) to the same order of magnitude as a E (•, •). Typical choices of S E (•, •) and α E take the form:
and
where h E . = |E| 1/3 represents the size of element E and x v stands for the vertices that belong to E. This gives the final form of the element-level discrete bilinear form as:
The first and second terms of a E h (u h , v h ) are respectively known as the consistency and stability terms, and they are responsible for the satisfaction of the two key conditions, namely consistency and stability, respectively, to ensure the convergence of the VEM approximation (Beirão da Veiga et al. 2013).
On the other hand, we approximate the continuous loading term (i.e., the loading term with exact integration) (v h ) as Chi et al. (2017):
where We are now ready to state the final form of the VEM approximation for 3D linear elasticity problems, which consists of finding u h ∈ K h such that:
where K 0 h is a subspace of K h with functions that vanish on u h .
For a given element E, we consider a set of basis functions, ϕ 1 , ..., ϕ 3m , for the local displacement space [V(E)] 3 of the form:
where we recall ϕ 1 , ..., ϕ m is the set of basis functions for V(E). With the set of basis functions, any displacement field
3 can be interpolated using its DOFs (the values at the vertices) as:
Moreover, by definition, the projection ∇ E ϕ j , j = 1, ..., 3m is given by:
According to (34), we evaluate the (j, k)th component of the element stiffness matrix k E ∈ R 3m×3m as:
More specifically, we can define matrices B E ∈ R 6×3m and D ∈ R 6×6 of the form:
and rewrite the first term (40) in matrix form as:
For the second term of (40), we can express it in matrix notation as:
where S E ∈ R 3m×3m is given by:
with P E defined in (25).
Having computed the local stiffness matrix for each element, we can obtain the global stiffness matrix K = E k E through the standard assembly procedure and solve the linear system of equation:
for the nodal displacement vector U, where F is the external force vector. For a given discretization, the external force vector F contains the nodal loads computed from applied traction t and body force f based on (35).
In this subsection, the performance of the VEM approximations on convergence and accuracy is assessed through a benchmark problem. In the assessment, we evaluate two error measures of the displacement and stress solutions. The measure of error in the displacement solution is defined as:
where u is the exact displacement solution and u h denotes the displacement field obtained by interpolating the VEM DOFs using the 3D Wachspress shape functions (Floater et al. 2014). On the other hand, we also define the L 2 norm of the stress error as:
where σ is the exact stress solution and σ h is a piecewise constant stress field defined such that:
For both error measures, the integrals are evaluated using a fourth-order integration rule on each tetrahedral subdivision of E.
We consider a boundary value problem in which a cantilever beam is loaded by end shear. As described in Fig. 2 (a), the beam occupies domain = (-1, 1)×(-1, 1)×(0, 10) and is subjected to a constant traction t = [0, -τ, 0] T on its top surface. According to Barber (2010), the stress solution of such a problem is given by the following expressions:
sin(nπx) sinh(nπy)
where E Y and ν are Young's modulus and Poisson's ratio, respectively. For the stress solution, we can see that the beam is traction-free on its four lateral surfaces. Additionally, the displacement solution that corresponds to the above stress distributions, up to the addition of a rigid body motion, is given by:
where w(x) is the anti-derivative of σ yz with respect to y. In our numerical study, the material properties of the solid are taken to be E Y = 25 and ν = 0.3, and the magnitude of the shear load is chosen to be τ = 0.1. We apply the analytical displacement field given in (51) on the bottom surface and analytical traction (calculated from ( 50)) on the top surface. Three families of structured meshes made up of regular space-filling polyhedra (i.e., regular hexahedra, truncated octahedra, and rhombic dodecahedra) and an additional family of CVT meshes are considered, as shown in Fig. 2 b-e. The convergence of both displacement error, u , and stress error, σ , as functions of the average mesh size h is depicted in Fig. 3 a and b, respectively. The results confirm that optimal convergence rates are obtained for both displacement and stress errors-the rate of convergence for u is 2 and that for σ is 1.
This section introduces a new VEM-based topology optimization framework on polyhedral meshes, which features continuous design and material density functions.
We shall focus on the classic compliance minimization problems and remark that the proposed framework is readily applicable to other optimization formations. We also note that, from now on, the body force f is neglected. For a given discretization h consisting of nonoverlapping polyhedra, the topology optimization formulation for the minimum compliance problems is stated as:
where ρ h is a material density function, A h is the space of admissible designs, and V is the allowable volume fraction. To regularize the formulation, the material density function ρ h is defined as the image of a design function η h under a map F h (e.g., the density filter) and the DOFs of the design function η h are the design variables (DVs). Moreover, a material interpolation function m S (ρ h ) is employed to relate the material stiffness to the value of ρ h at any given point. For instance, if the SIMP model is used (Bendsøe (1989), Rozvany et al. (1992), and Bendsøe and Sigmund (1999)), we have m S (ρ
, where ε is the Ersatz parameter and p is the penalization parameter. Incorporating the spatially varying stiffness m S (ρ h ), a ρ h (u h , v h ) is the discrete bilinear form constructed using VEM.
The topology optimization formulation (52) can be viewed as a two-field mixed approximation problem involving a discrete displacement space K h and a discrete design space A h (Jog and Haber 1996). The standard densitybased topology optimization framework in the literature typically employs a continuous displacement field, whose DOFs are the displacements at the vertices of the mesh, and a piecewise-constant design function, whose value in each element is the associated DV. The material density within each element takes a constant value as well. Several examples of mixed elements of this type are shown in Fig. 4a. Roughly speaking, for a given discretization, the number of DVs governs the "resolution" of the topology whereas the number of the displacement DOFs (the size of the state equation) dictates the computational cost. If we consider a unit cube and discretize it with the mixed elements shown in Fig. 4a, we can then plot in Fig. 4b the total numbers of DVs as functions of the total numbers of displacement DOFs when those meshes are refined. We observe that those functions are close to linear and their slopes can be used to quantify the computational efficiency of the topology optimization framework on various discretizations: the larger the slope is, the more computationally efficient the discretization is. As shown from the slopes in Fig. 4b, the mixed approximation in the standard topology optimization framework leads to considerably more displacement DOFs than of DVs on various discretizations, which is undesirable from a computational efficiency perspective. Moreover, Fig. 4b also suggests that polyhedral discretizations yield smaller slopes as compared to the hexahedral ones. This observation indicates that, although polyhedral discretizations exhibit several geometric advantages in topology optimization (Gain et al. 2015), they are less computationally efficient than hexahedral ones in the standard topology optimization framework. Motivated by the above observations and discussions, we propose a new approximation of the design function (as well as the material density function) for topology optimization on general polyhedral meshes. The basic idea is to consider a more enriched local design space with continuous design functions in each element. In terms of DVs, they are placed at the vertices as well as the mid-edge nodes of the meshes. On the other hand, the displacement approximation is kept the same as in the standard case. This leads to a new mixed approximation for topology optimization on general polyhedral meshes. We note that the idea of enriching the design space is conceptually similar to the multi-resolution methodologies (Nguyen et al. 2010) in topology optimization. An illustration of several new mixed elements of this type is shown in Fig. 5a. Again, considering a unit cube discretized with those mixed elements, we plot in Fig. 5b the total numbers of DVs as functions of the total numbers of displacement DOFs with mesh refinement. By comparing the slopes of those functions in Figs. 4b and5b, we conclude that the proposed mixed approximation improves the computational efficiency as compared to the standard topology optimization framework. More specifically, for a given discretization with a fixed number of displacement DOFs (thus with a roughly fixed computational cost), the proposed mixed approximation yields a more enriched design field and, thus, can produce topologies with improved resolutions.
This subsection formalizes the proposed topology optimization framework. In this framework, both the discrete displacement and design fields are constructed using the VEM. As we can see, because of its flexibility in handling any element geometry as well as local spaces, the VEM provides an efficient platform to formulate the proposed topology optimization framework.
To introduce the new topology optimization framework, we define the space A h for admissible density function ρ h as:
In the above definition, we assume that both density function ρ h and design function η h are continuous functions with their DOFs being their values at vertices as well as midedge nodes of the mesh. The density function ρ h is defined as the map of a design function η h by the density filter F h (•). For a discretization h consisting of M elements and N D DVs, we henceforth introduce two vectors, z = [z 1 , z 2 , ..., z N D ] T and y = [y 1 , y 2 , ..., y N D ] T , where z i and y i are values of η h and ρ h at the ith DOF of η h and ρ h , respectively.
Since both functions ρ h and η h can be characterized by their DOFs, we adopt the following approach to construct the density filter map F h (•) between ρ h and η h , which is based on their DOFs and the associated position vectors. If we denote by S(i) the set of the indices of DOFs whose positions fall within a sphere of prescribed radius R centered at x i (the position vector associated with DOF i), y i (the ith DOF of the density field ρ h ) is computed as:
where q is the order of the density filter (Bourdin 2001;Zegard and Paulino 2016). Figure 6 illustrates this density filter on a 2D mesh patch. For easier implementation, we express the density filter mapping between ρ h and η h in a matrix form as:
where P F is a constant and sparse matrix with its (i, j )th component given by:
Moreover, within each element E (assuming E has m vertices), the density function ρ h belongs to the VEM space V(E) defined in (11). Likewise, if we use y E i to denote the ith DOFs of ρ h in E, we can express ρ h | E in terms of the set of basis functions for V(E), ϕ 1 , ..., ϕ m , as:
The volume constraint function in (52) can then be recast as:
where we recall from ( 26) that E ϕ i , i = 1, ..., m can be exactly computed by the definition of V(E). To assist easier implementation, we also express the volume constraint function in matrix notation as:
where
T is a vector collecting element volumes and P V is a constant matrix with its (i, j )th component being:
Fig. 6 An illustration of a quadratic (i.e., q = 2) density filter mapping for a 2D mesh patch. The DVs with solid markers are within the filter radius R and thus have non-zero weights
We also note from the above definition that matrix P V is sparse, and thus it is formed and stored as a sparse matrix in our implementation.
In the state equation, the interpolated stiffness function m S (ρ h ) is utilized, which is assumed to be an element-wise constant function, such that within element E, m S (ρ h )| E . = m S (< ρ h > E ), where < ρ h > E denotes the volume average of ρ h over E:
Having defined the form of the stiffness interpolation function, we propose to incorporate it in the element-level discrete bilinear form (34) as follows:
Notice that the first term (the consistency term) on the righthand side of the above equation is now a constant function of the form:
Furthermore, we can pull m S (< ρ h > E ) out of the above element-level bilinear form and get the global discrete bilinear form a ρ h (u h , v h ) as:
Similarly, we introduce a vector s in the implementation whose ith component s i is the value of m S (< ρ h > E ) for element E i . Utilizing the matrix P V , the vector s can be expressed as:
We can then compute the global stiffness matrix for a ρ h (m S (ρ h ), u h , v h ) as:
Finally, we arrive at the topology optimization formulation considering nodal densities as:
We next describe how to compute the gradients of the objective and volume constraint functions with respect to the design variable z. For the objective function, we first compute its gradient with respect to the vector s as:
and then, using the chain rule, we arrive at:
where J m S (y) .
= diag(m S (y 1 ), ..., m S (y N D )) is the Jacobian matrix of the stiffness interpolation function m S . The
We conclude this subsection by noting that, even with a more enriched space for design and material density functions, the proposed formulation (67), which is formulated in the VEM context using local projections, can be implemented in a similar code structure to the PolyTop software (Talischi et al. 2012b).
In this subsection, we verify the proposed continuous density approximation through a simple numerical example and compare its performance with the standard element-wise density approximation. In particular, we want to compare how these two approximations perform in terms of capturing the varying density in a simple boundary value problem.
The setup of the boundary value problem is as follows. As shown in Fig. 7a, we consider a unit cube and place two (ρ = 0) tubes in it. The rest of the cube is solid (ρ = 1). The tube whose central axis passing through points [0.22, 0.4, 0] T and [1, 0.7, 0.78] T has a radius of R 1 = 0.13, and the other one whose central axis passing through [0.5, 0, 0.3] T and [0.3, 0.6, 1] T has a radius of R 2 = 0.06. For this setup, we can introduce a continuous parametrization ρ(x) of the material density over using smooth Heaviside functions as:
where ε is a small positive number assigned to ensure the positivity of ρ(x), and r 1 (x) and r 2 (x) are the distances of a given point x to the central axes of the first and second tubes, respectively. As illustrated in Fig. 7b, for a tube whose central axis passing through two given points x 1 and x 2 , r(x) is given by: where || • || stands for the Euclidean norm. Making use of ρ(x), the (continuous) distribution of the elasticity modulus in is given by C(x) = ρ(x)C 0 , where C 0 is the elasticity modulus tensor of the solid material. In the present study, the solid material is considered to be isotropic with Young's modulus being E Y = 25 and Poisson's ratio being ν = 0.3. Adopting the method of manufactured solutions, we assume an exact displacement solution u of the form: u x = z 2 e xy u y = 2y 2 z 3 + zx 4 and u z = z sin(2πx)e y , (73) and prescribe it on the entire boundary of the cube. Accordingly, a body force, which is computed using the exact displacement solution u and stiffness distribution C(x), is prescribed in the interior of . Four families of polyhedral meshes are considered: hexahedral, truncated octahedral, rhombic dodecahedral, and CVT; and their examples are shown in Fig. 7c-f, respectively. To assess the accuracy of the solution, we make use of the L 2 error of the displacement defined in ( 47) and a H 1 displacement error, which is defined as:
where the integral is evaluated using the same fourth-order integration rule on each tetrahedral subdivision of E.
In the numerical simulations, the material distribution ρ(x) needs to be approximated in order to compute the stiffness matrix. We consider two approaches here. The first approach assigns a constant density to each element with the value being ρ(x) evaluated at the centroid of that element. This approach resembles the situation in the standard density-based topology optimization where element-wise constant densities are used. On the other hand, the second approach evaluates ρ(x) at the vertices and mid-edge nodes, and interpolates them using VEM basis functions. This approach resembles the situation in the proposed topology optimization framework where the DVs are located at the vertices as well as mid-edge nodes of the mesh. Once the material distribution is approximated, the stiffness matrix of the discretized system for both approaches can be constructed. For the former approach, the stiffness matrix is constructed using the standard procedure, whereas, in the latter approach, the stiffness matrix is formed following (64).
Figure 8 a-d show the convergence of both L 2 and H 1 displacement errors as functions of mesh size h for both approaches on the four families of meshes. The standard elements represent the first approach and the new elements stand for the second approach. As we can see from the comparisons, the second approach (corresponding to the proposed topology optimization framework) always gives more accurate displacement solutions than the first approach (corresponding to the standard density-based framework) on a given mesh. This in turn implies that the material densities interpolated from values on both vertices and mid-edge nodes provide better approximations of the exact density distribution ρ(x) than the one constructed from elementwise constant values. Based on this observation, we argue that, in topology optimization, it is also more favorable to use the proposed material density approximation because it gives a better parametrization of the density distribution of the varying topologies. We also remark that, although Fig. 8 a-d suggest that the regular hexaheral and CVT meshes produce more accurate displacement solutions than other two families of meshes, this conclusion could be problemdependent and might not hold for other density distributions different from the one defined by (71).
In this section, we present several examples involving non-Cartesian design domains to demonstrate the effectiveness and versatility of the proposed VEM-based topology optimization framework. For all the design examples, we select the Ersatz parameter as ε = 10 -9 and take Young's modulus and Poisson's ratio of the solid phase to be E Y = 100 and ν = 0.3, respectively. During optimization, the optimality criterion (OC) (Christensen and Klarbring 2009) is chosen as the design-variable update scheme with the damping parameter and move limit being η = 0.5 and move = 0.3. The maximum tolerance for the change of design variables in the convergence criterion is taken to be 0.1%. Moreover, a continuation scheme of the penalization parameter p is adopted here. We initialize p as p = 1 and increase it every (maximum) 20 optimization iterations by 1 until p = 3. When p reaches 3, we then set the maximum allowable optimization iteration number to be 150. For comparison purposes, the topologies obtained by the standard topology optimization framework (Gain et al. 2015) with element-wise constant densities are also provided. The same VEM formulation is used to solve the state equation in the standard framework as in the proposed framework. Unless otherwise stated, the final topologies are plotted using iso-surface with the cutoff value being 0.5. To distinguish the topologies obtained using the standard framework and the proposed one, we plot the results obtained using the standard framework in blue and those obtained using the proposed framework in red.
We also point out that, for a fixed mesh and fixed set of parameters, the computational times of both proposed and standard frameworks are closely identical because the sizes of state equations are the same.
The first design example is the shear loaded disk problem.
As shown in Fig. 9, the design domain is a disk with an outer radius of 6 and an inner radius of 1. The thickness of the disk is taken to be 2. Fixed in the inner surface, the Fig. 10 The final topologies obtained from the proposed topology optimization without the density filter on a the hexahedral-dominated mesh; b the truncated octahedral-dominated mesh; c the rhombic dodecahedral-dominated mesh; and e the CVT mesh Fig. 11 The final topologies obtained from the standard topology optimization without the density filter on a the hexahedral-dominated mesh; b the truncated octahedral-dominated mesh; c the rhombic dodecahedral-dominated mesh; and d the CVT mesh. In the designs, those elements whose densities are above 0.01 are plotted circumference of this disk is subjected to eight equidistant shear loads of uniform magnitude 1. A volume fraction of 10% is prescribed.
In order to provide a thorough assessment of the performance of the proposed framework on various types of meshes, we consider four meshes in this design example: a hexahedraldominated mesh, a truncated octahedral-dominated mesh, a rhombic dodecaheral-dominated mesh, and a CVT mesh, as shown in Fig. 9b-e, respectively. The first three meshes consist of regular space-filling polyhedra in the interior of the design domain and unstructured polyhedra in the boundary regions, while the CVT mesh is made up of unstructured polyhedra inside the entire design domain. The statistics of this mesh is provided in Table 1. Notice that the four meshes have the similar numbers of elements.
We first design the shear loaded disk problem without applying the density filter. Figures 10 and11 show the final topologies obtained from the proposed and the standard frameworks, respectively, on the four meshes. Notice that, for the results obtained with the standard topology optimization framework, we get almost black and white designs for all the four meshes. Thus, for those results, instead of showing the iso-surface plots of the final topology, we simply plot those elements whose densities are above 0.01. In the topology optimization literature, it is well known that, without density or sensitive filters, quadrilateral or hexahedral meshes will produce checkerboard patterns in the standard framework (e.g., see Fig. 11a). For general discretizations, a previous work has demonstrated that 2D regular hexagonal (Talischi et al. 2009) and 2D and 3D CVT meshes (Talischi et al. 2010(Talischi et al. , 2012a;;Gain et al. 2015) are free of checkerboard patterns. From our numerical results in Fig. 11c for truncated octahedral and Fig. 11d for rhombic Fig. 12 The final topologies obtained from the proposed topology optimization with a quadratic density filter of radius R = 0.45 on a the hexahedral-dominated mesh; b the truncated octahedral-dominated mesh; c the rhombic dodecahedral-dominated mesh; and d the CVT mesh Fig. 13 The final topologies obtained from the standard topology optimization with a quadratic density filter of radius R = 0.45 on a the hexahedral-dominated mesh; b the truncated octahedral-dominated mesh; c the rhombic dodecahedral-dominated mesh; and d the CVT mesh dodecahedra, we observe that, unlike the regular polygonal discretization in 2D, regular polyhedral discretization in 3D tends to generate patterns indicating numerical instability. Those patterns provide artificial stiffness and are impractical from a manufacturing perspective. On the contrary, with the proposed topology optimization Fig. 14 The convergence history of objective functions for both proposed and standard frameworks on a hexahedraldominated mesh; b truncated octahedral-dominated mesh; c rhombic dodecahedraldominated mesh; and d CVT mesh. Both frameworks deliver almost identical convergence rates. We also emphasize that it is unfair to compare the absolute values of objective functions between the two frameworks because different density approximations are adopted Fig. 15 a The geometry, load, and boundary conditions of the hook problem. b A hexahedraldominated mesh consisting of 134,634 nodes and 63,350 elements. c A hexahedraldominated mesh generated by rotating the seeds of the mesh in b counterclockwise 60°in the xy plane. This hexahedral mesh contains 132,457 nodes and 63,320 elements framework, we obtain physical designs on all the four meshes considered, which resemble a flower, although their iso-surface plots exhibit rough boundaries because of the absence of regularization from the density filter.
We then apply a quadratic density filter (i.e., q = 2) with a radius of R = 0.45, and keep the other parameter settings unchanged. The results obtained are shown in Figs. 12 and 13 using the proposed and standard frameworks, respectively. Several observations can be made. First, compared to the ones obtained from the standard framework, the designs obtained from the proposed framework have iso-surface presentations with greatly improved smoothness. This is because, with the proposed framework, we have more DVs as compared to the standard framework on the same discretization, as shown in Table 1. Second, unlike the ones obtained using the standard framework, the optimal topologies obtained using the proposed framework exhibit similar topologies on all meshes, suggesting that the proposed framework is less sensitive to mesh types than the standard framework. Moreover, Fig. 14 a-d depict the convergence history of the objective functions for both the proposed and standard frameworks on hexahedral-dominated, truncated octahedral-dominated, rhombic dodecahedral-dominated, and CVT meshes, respectively. The comparison in those figures suggests that the proposed and standard frameworks deliver almost identical convergence rate on various types of polyhedral meshes. Because, for any given mesh, the size of the stiffness matrix is identical for both frameworks, the major computational costs for the standard and proposed frameworks are identical as well. However, as shown in Table 1, the proposed framework can handle significantly more numbers of design variables than the standard one on all the meshes considered, in this sense, that the proposed framework is more efficient than the standard one. Finally, we have to emphasize that it is unfair to compare the absolute values of the objective function between the standard and proposed frameworks. Even on the same mesh, different approximations of the density field (especially for cases where the SIMP penalization factor p > 1) will surely yield different displacement solutions and, eventually, lead to different absolute values of the objective function.
Having investigated the performance of the proposed topology optimization framework through the last design example, we now apply it to several problems involving complex design domains which are non-Cartesian. Let us first look at the hook problem. As shown in Fig. 15a, the design domain in this problem is fixed in the upper half of the circle and is subjected to a uniformly distributed line load in the negative z direction. The volume fraction of this problem is set as V = 15% and the radius of the quadratic density filter is R = 2.5. Through this design problem, we aim to quantify the mesh bias in both the standard and proposed topology optimization frameworks and demonstrate that the proposed one is less biased to the initial mesh than the standard one.
To that end, we consider two similar polyhedral meshes of the design domain. The first mesh, shown in Fig. 15b, is generated using the Voronoi seeds with Cartesian alignment; and the second mesh, shown in Fig. 15c, is generated using the seeds alignment obtained by rotating the seeds of the first mesh counterclockwise 60 • in the xy plane. For both meshes, the seeds in regions near the boundary are then updated using Lloyd's algorithm (Thedin et al. 2014) in order to capture the boundary geometry of the design domain. As a result, both meshes are made up of regular hexahedra in the interior and unstructured polyhedra in regions near the boundary, and contain similar numbers of nodes and elements, as summarized in Table 2. However, the two meshes possess different preferred orientations, as shown in Fig. 15b andc. In Fig. 16a and b, we show the converged designs from both the proposed and standard topology optimization frameworks for the original and rotated meshes, respectively. Our main conclusion from comparing Fig. 16a andb is that the designs obtained from the proposed framework are less biased to the initial meshes than the standard framework. To visualize the bias, let us compare the fan regions of the designs obtained from the two frameworks. For the designs obtained by the proposed framework on both original and rotated meshes, the fan regions resemble the one obtained from a 2D analogue of this problem (Talischi et al. 2012b). The orientations of the members in the fan regions seem to be not influenced by the preferred orientations of the initial meshes. In contrast, for the designs obtained by the standard framework, it is clear that the orientations of the members in the fan regions are biased toward the preferred orientations of the initial meshes. For example, the original mesh has preferred orientations along the x and y axes. As shown in the blue design Fig. 16 a The final topologies for the hook problem obtained from both the proposed and standard frameworks on the original mesh (c.f. Fig. 15b). b The final topologies for the hook problem obtained from both the proposed and standard frameworks on the rotated mesh (c.f. Fig. 15c) in Fig. 16a, several members in the fan region of its design are clearly biased toward these two orientations. In addition, comparing the final designs obtained by the standard and proposed frameworks, we also notice that the ones obtained from the proposed framework possess smoother iso-surface representations than the ones obtained from the standard framework.
We also investigate the design of a wrench problem. The design domain of this problem is depicted in Fig. 17a. In the domain, the bigger circle is fixed and a half of the other circle is subjected to a distributed line load along the negative y direction. Here, we consider two CVT meshes, a coarse one and a refined one, whose statistics are given in Table 3. Again, a quadratic density filter is applied with a radius of R = 0.05 and the volume fraction is prescribed as V = 15%. In this design example, we prescribe another requirement that the final topology has to be symmetric. To achieve this, we use a matrix P S introduced in Talischi et al. (2012b), which is given by:
1 if 1) the y coordinate of the kth DOF is non-negative; and 2) k = or the kth and th DOFs are y -symmetric 0 otherwise (75)
and the vector y and s are then given by: y = P S P F z and s = m S (P V P S P F z)
respectively. As we can see, the matrix P S enforces symmetry through mapping the admissible topologies in the design space to symmetric configurations. Figure 18 a and b depict the final topologies obtained using the proposed and standard frameworks, respectively, for the coarse CVT mesh; and Fig. 19 a andb show the ones for the refined CVT mesh. As an immediate observation, the proposed and standard frameworks yield symmetric 1 designs similar to each other for both the coarse and refined meshes. With the coarse mesh, the proposed framework yields a clearer and manufacturable design, which resembles the 2D optimization result in Talischi et al. (2012b), than the standard framework, indicating that the 1 Although the CVT meshes are generated from reflected seeds to ensure their symmetry. Certain regions on the mesh boundaries, especially near the two circles, are not fully symmetric due to the limitation of the meshing software used by Thedin et al. (2014) in representing curved boundaries. Hence, results obtained from the proposed framework (i.e., Figs 18 and19a) are slightly asymmetric in those regions. However, we note that this minor issue will not affect the quality of the designs and can be resolved when an improved version of meshing software is used. proposed framework is more effective on relatively coarser meshes.
In the final design example, we perform topology optimization in a serpentine domain. As shown in Fig. 20a, the serpentine domain is fixed on its left face and is subjected to a point load along the negative z direction in the middle of the lower edge of the right face. A volume fraction of V = 10% is prescribed and a quadratic density filter of radius R = 0.4 is used. As in the hook example, we consider a mesh composed of regular hexahedra in the interior and unstructured polyhedra in regions near the boundary. The mesh is plotted in Fig. 20b and its statistics is presented in Table 4. In Figs. 20c and 21, we show both the final design produced by the proposed framework and the manufactured design using FDM 3D printing, respectively. This example demonstrates that the proposed topology optimization framework can lead to designs that are directly manufacturable (Zegard and Paulino 2016).
In this work, we present a 3D VEM-based topology optimization framework on general polyhedral discretizations.
The unique feature of this work is that it takes full advantage of the VEM and applies it to both structural and optimization problems. In terms of structural problems, the VEM is adopted to solve the state equation efficiently and effectively. Because VEM does not require explicit computations of the shape functions and their gradients, it does not need numerical integration and is less sensitive to degenerated polyhedra (e.g., ones with skinny faces or small edges) as compared to the FEM. In terms of optimization problems, exploiting the great flexibility of VEM in element geometries and local space definitions, we introduce an enhanced VEM space for the continuous design and material density functions, which contain DOFs at the vertices as well as the mid-edge nodes of the mesh. The total volume of any topologies in this design space can be computed exactly using a properly defined VEM projection operator. As a result, for a given mesh and under a similar computational cost, the proposed VEM-based topology optimization is shown to produce designs with improved geometrical resolutions as compared to the standard topology optimization framework with element-wise constant DVs and material densities. We also demonstrate that the VEM-based topology optimization framework can be implemented in a similar code structure to the PolyTop software (Talischi et al. 2012a). In terms of discretizing complex domains in 3D, this work explores two approaches: unstructured polyhedral (i.e., CVT) meshes and regular polyhedra-(i.e., hexahedra, truncated octahedra, and rhombic dodecahedra) dominated meshes.
Both of the VEM approximations for the displacement field and material density function are verified through numerical examples. The convergence of the VEM in solving the state equations is verified via a benchmark problem involving beam bending. Through a simple example, we also compare the performance of the proposed continuous density parametrization with the commonly used element-wise constant one in terms of approximating varying material density distributions. Moreover, several design examples involving non-Cartesian domains are presented, showcasing that the proposed VEM-based topology optimization framework produces designs with improved quality and achieves higher computational efficiency. In terms of design quality, we demonstrate through the hook and wrench design examples that the proposed framework is less biased to the initial mesh (see the results in Fig. 16) and can produce designs with smoother iso-surface representations, especially for smaller meshes (e.g., see the comparison results in Fig. 18). In terms of computational efficiency, we demonstrate that the proposed framework delivers the same convergence rate as the standard one on various types of discretizations (see Fig. 14). For a given mesh, because the size of the stiffness matrix is identical for both frameworks, the proposed framework is able to consider a significantly bigger number of design variables than the standard one with almost identical computational effort, and thus achieves a higher efficiency. Additionally, we show that the proposed framework is flexible in imposing various pattern and layout constraints (i.e., symmetry constraints shown in Figs. 18 and 19), and can lead to designs that are directly manufacturable by 3D printing (see Fig. 21).
Finally, we remark that the proposed VEM-based topology optimization offers an effective tool for mesh adaptation in topology optimization. The potential of this research includes developing efficient mesh adaption strategies for topology optimization by exploiting the advantages of polyhedral elements. The times are in seconds with percentage of total runtime provided in parentheses out on a desktop computer with an Intel(R) Xeon(R), 3.00 GHz processor, and 256 GB of RAM running MATLAB R2016a. For all the meshes, the two codes produce almost identical final topologies and are not shown here for the sake of conciseness.
In Table 6, we present a comparison of the total runtimes of PolyTop3D and top3d for the three meshes. In addition, Table 7 shows the breakdown of the total runtime of the PolyTop3D code into major steps. One immediate conclusion from Tables 6 and7 is that the PolyTop3D