A shape function based interpolation approach for nodal design variable based structural optimization

In this article, a novel density interpolation scheme for topological optimization based on nodal density variables is presented. The method relies on the definition of an independent interpolation mesh for the description of the density field. This mesh is used for the interpolation of the density field at the Gauss points of the design domain mesh. Provided that the interpolation mesh is realized by linear elements, the proposed scheme is range restricted and assures physically meaningful values of the interpolated density. Explicit analytical sensitivity expressions can be derived. The utilization of an independent interpolation mesh allows for the utilization of any element type and geometry on the body for the computation of the displacements. It is also shown that the chosen interpolation scheme has interesting properties of volume and material volume conservation with respect to the mesh utilized on the body and the domain of the interpolation mesh can differ from the domain of the mesh for the displacements computation. This allows for the realization of very regular interpolation meshes even in case of very complex shapes of the domain of the body to be optimized. Examples are provided showing the performance of the proposed interpolation scheme.

element densities are no longer design variables but are computed by the values of densities at the nodes. Element densities can be considered uniform or nonuniform inside each element. 10 In the continuous approximation of material distribution (CAMD) approach, 6 the density design variables at the nodes are interpolated by means of shape functions. The resulting density field is C 0 and has implicit gradient limitations that avoid the checkerboard instability. Different schemes can be employed for the interpolation of the density field with respect to the mesh used for the interpolation of the displacements. In the Q4/Q4 scheme, the two meshes share the same nodes, but schemes with different nodes arrangement can be considered. 6,10,11 However, in all the published papers, the mesh used for the interpolation of the density field is derived from the base mesh used for the interpolation of the displacements. CAMD is free from checkerboard instabilities, but it suffers from instabilities due to "layering" and "islanding." 6 These instabilities seem to be related to the utilization of meshes too coarse to catch the material distribution and can be avoided with mesh refinement. However, in some cases, no convergence with mesh refinement has been found. 6 Convergence problems can be avoided by considering constraints on the density distribution, such as perimeter constraints. However, the possibility to avoid the utilization of filters, allows for the appearance of fine structures that can provide important details otherwise lost. An important feature of the CAMD approach is the possibility to obtain a higher resolution in the topology design without refining the mesh used for the computation of the displacements. 10,12 In the projection of the nodal design variable, the density value at each node is projected to the elements contained into a influence volume. The density value at each element is then computed by considering the density values of all the node influencing that element. Different interpolation schemes can be used, usually linear or nonlinear function of the distance between the element and each influencing nodes. 8,9,13,14 The influence distance performs similarly to a density filter and avoids checkerboard instability. These methods are usually quite mesh insensitive, but the distance parameter influences the solution. 9 This approach can be applied also to meshless optimization methods, which do not require mesh connectivity nevertheless maintaining a sufficient accuracy and numerical stability. 15,16 Interpolation or projection methods can be used to realize mesh adaptive [17][18][19] or multiresolution 12,20-22 topology optimization approaches to increase the results definition while reducing computational time. These methods decouple the discretization on the design variables from the discretization for the displacement field to obtain high resolution results with relatively coarse meshes. Isogeometric topology optimization approaches are considered in References 21 and 23 to improve the geometrical consistency between the FE model and the CAD model.
The mentioned approaches are able to effectively solve the optimization problem by considering nodal design variables and effectively reduce or eliminate the numerical instabilities. However, if the interpolation or projection scheme is based on the same mesh that is used for the computation of the displacements, very regular elements have to be employed. In this case, even simple geometries require partitioning or complex mesh transformations. 23 Projection schemes which use distance-based weights, such as Shepard 9 or Heaviside 13 interpolation, present some interpolation issues such as the non-consideration of the direction of interpolation and numerical problems for particular neighborhoods. 24 Also, these methods rely on the choice of the distance term which strongly influences the results.
In this article, a novel interpolation approach is presented. In this proposed approach, the nodal design variables are interpolated by shape functions on a mesh completely independent from the mesh used for the computation of the displacements. This kind of interpolation allows for a complete freedom in the definition of the displacement mesh. Also, the formulation will be given as to allow the realization of regular interpolation meshes in a very simple way and with no geometrical restrictions on the optimization domain. The proposed interpolation can be applied to a series of design problems characterized by complex design spaces where a non-structured high order mesh is required for stress computation, in case with body remeshing. Other envisaged design applications are optimization of bodies sharing the design space, 25,26 combined shape and structural optimization, 27 embedded design areas, 28 and embedded bodies. 29 The article is organized as follows. First, the formulation of the topology optimization problem is recalled. Then, the proposed interpolation scheme is presented and the sensitivities are computed. In the last part of the article, some numerical examples are given.

TOPOLOGY OPTIMIZATION PROBLEM FORMULATION
The most common topology optimization problem consists in the definition of the best material distribution to achieve the minimum compliance of a component for a given mass constraint. 1,2 Among the various approaches for topological optimization, in this article, the well-known solid isotropic microstructure (or material) with penalization (SIMP) problem formulation 4 is considered. The main focus of this article is the interpolation of the density field, therefore the SIMP solution algorithm proposed by Sigmund 2 is implemented with only the adjustments required by the nodal density variable approach. The numerical formulation of the problem in the framework of the finite element analysis reads 2 Find min ∶ C ( ) = U T KU, where C is the compliance, U the nodal displacements, K the stiffness matrix of the body, V 0 the design volume, f the volume fraction, the density, and min the (small) value of the minimum considered density, greater than zero to avoid numerical singularities. If a constant density is considered for each element, the compliance can be computed as with n e number of elements, u e nodal displacements of each element, K 0 e reference element stiffness matrix, V ( ) material volume, e is the element density. K 0 e is computed for the reference values of the material properties. If the density is not constant over each element, the compliance formulation can be modified as where GP is the density value at each Gauss point used for the integration of the element stiffness matrix K e . In this article, the problem of Equation (1) is solved by the optimality criteria (OC) method and the new density field is computed as 2,3 where m is a positive move limit and B represents the optimality condition The Lagrange multiplier is computed by a trust-region Dogleg algorithm. is a numerical damping coefficient assumed equal to 0.5. 2 For the implementation of the OC method, the two sensitivities V (sensitivity of the material volume with respect to the design variables) and C (sensitivity of the compliance with respect to the design variables) have to be computed. The computation of these sensitivities has to be completed considering the interpolation scheme adopted for the nodal densities.
It can be observed that the interpolation of the nodal densities does not affect the optimization algorithm. Therefore, mesh dependency 3 and dependency of the solution from the initial conditions 3 typical of the SIMP method are expected. Also numerical artifacts are expected. However, instead of the checkerboard instabilities of the usual formulation of the SIMP method, "islanding" and "layering" are expected due to use of a nodal density field with interpolation. 6

INTERPOLATION OF THE DENSITY FIELD
The interpolation scheme for the nodal values of the density field proposed in this article is based on the shape functions approach used for the finite elements. However, the interpolation is not based on the elements of the displacement mesh as usually employed with shape functions, 6,8 but it is independent from this mesh. The interpolation is obtained by considering two completely different and independent meshes, one for the finite element analysis of the body and one for the description of the density field. The two meshes do not share any node and can be made of elements of different type and dimension. For clarity, the following nomenclature is considered. The mesh on the body and used for the finite element analysis is called body mesh. The quantities referred to the elements of the body mesh are indicated with the subscript e. The density of each element of the body mesh is indicated as e if it is referred to its center or GP if it is referred to one of its Gauss points. The mesh for the interpolation of the density field is called interpolation mesh. The quantities referring to each node of the interpolation mesh are indicated with the subscript i. The density at each node of the interpolation grid is indicated as i . The densities i are the design variables.
Any element type can be used for the body mesh. For the interpolation mesh, to have physically meaningful values of interpolated density values, only linear, bilinear (2D), or trilinear (3D) elements have to be employed. 9 In fact, higher order elements do not have range-restriction properties and the interpolated value of density can exceed the range 0-1 even if the nodal value is in this range.
The two considered meshes are depicted in Figure 1. In this figure, the body mesh is realized by quadratic triangles, while the interpolation mesh by bilinear quadrilaterals. The characteristic dimension of the body mesh is more than three times the characteristic dimension of the interpolation mesh. The nodes of the interpolation mesh are marked with red points and represent the location of the design variables. The blue squares represent the point where the density field is evaluated on the body. These values are then averaged inside each element of the body mesh and this value of density is assigned to the center of the element. Any arrangement of points could be used for the evaluation of density, however, to have a simpler averaging, the points are located accordingly to the Gauss point quadrature. To have a reasonable number of evaluation points in correspondence of any element of the interpolation mesh, in the example, a scheme with 13 Gauss points has been considered for the body mesh. In general, the number of Gauss points used for the evaluation of the density field can be different from the number of Gauss points used for the construction of the stiffness matrix of the element.
For any point x of the body, the density can be computed as where N i (x) are the shape functions of the element of the interpolation mesh containing x, n is the number of nodes of the element, and i are the pseudo-densities at the nodes of the element. Considering the coordinates x GP of the Gauss points of the body mesh, the density at each Gauss point reads F I G U R E 1 Interpolation scheme. The design variables (density) are evaluated at the nodes of the interpolation mesh (red points) In order to apply Equation (7), the element of the interpolation mesh containing the considered Gauss point has to be identified. In a general 3D case with unstructured interpolation meshes of any order this search can be quite complex, time consuming and special algorithms are required. However, the proposed interpolation scheme has been developed to have a regular interpolation mesh and to have physically meaningful values only first order elements have to be used. In this situation, for simple problems even a series of "if" conditions can be effectively used for the search. For large 3D meshes, the use of regular meshes speeds up the search process. Also, this search has to be realized once when the body mesh is generated at the beginning of the analysis and, in case, for any remeshing of the body and not at any iteration.
The element material volume is given by integration over the n GP Gauss points as where J GP is the Jacobian of the shape functions computed at each Gauss point and w GP is the corresponding weight of the Gauss quadrature rule. The density of each element can then be computed as being the element volume computed as The proposed interpolation scheme presents the following properties.
for ∀x, if linear or bilinear interpolation elements are used, that is, the interpolation values at any point are bounded between the minimum and maximum value of the design variables. This property ensures a physical meaning of the density at any point.
These properties can be also found in distance based interpolation schemes, such as the Shepard interpolation. 9 However, while distance based interpolation schemes consider only the distance of the interpolated points with respect to the data points, 24 the proposed approach considers also the direction. Also, distance based interpolations can present some oscillations in the interpolated data, while the proposed approach is constrained to linear interpolation between data which prevents any oscillation.

SENSITIVITY ANALYSIS
For the computation of the sensitivities, the set Ω i of the elements of the interpolation mesh sharing the node i has to be considered. Ω i is depicted in Figure 2. If a node of the interpolation mesh is considered (black point in Figure 2), its value of density i is used for the computation of the densities at all the Gauss points contained in the elements sharing this node. The densities interpolated at these Gauss points are then used for the computation of the density at the center of the elements of the body mesh containing the considered Gauss points (Equation 9). Therefore, i influences the density and the stiffness matrix of all the elements that have at least one Gauss point in Ω i . These elements have been highlighted in red in Figure 2. These dependencies have to be considered when computing the sensitivities.

Material volume sensitivity
The material volume of each element of the body mesh is computed by Equation (8). Let us consider each term of the summation on the Gauss points of Equation (8). Each of these terms represents the fraction of material volume pertaining to each Gauss Point of the body mesh and has expression The material volume of the body influenced by the density value at each node of the interpolation grid can be computed as The sensitivity of V Ω i with respect to i is given by where By replacing Equation (14) into Equation (13), the sensitivity of the material distribution with respect to the design variables reads It can be easily verified that and where Ω is the domain of the interpolation mesh and Ω b is the domain of the body. Equations (16) and (17) state that the body volume (v) and material volume (V) computed on the body mesh or on the interpolation mesh are equal. These two conditions can be exploited for the computation of the actual material volume of the body directly on the interpolation mesh, without any approximation. That is, the volume fraction constraint can be exactly imposed on the interpolation mesh. Moreover, it can be easily checked that by considering these equations, the domain Ω of the interpolation mesh does not have to be coincident with the domain Ω b of the body, but it suffices that Ω b ⊆ Ω. This condition can have a dramatic effect on the complexity of the optimization algorithm. In fact, regardless of the complexity of the geometry of the body, it is always possible to define an interpolation domain of very simple shape containing the domain of the body. Therefore, it is always possible to have very regular elements on the interpolation mesh even in the case of irregular body meshes. As a consequence, the construction of the interpolation mesh, the imposition of manufacturing constraints and the application of symmetries can be done in an easy way.
Finally, Equations (16) and (17) imply the possibility of remeshing the body without losing any information on the material distribution of the body. These properties can be very useful if a remeshing phase is included in the optimization algorithm.

Compliance computation and derivation
The compliance of each element of the body mesh can be computed as (Equation 3) where K 0 GP is the stiffness matrix contribution computed at each Gauss point for the reference material stiffness. The total compliance of the body reads The sensitivity of the compliance of each element with respect to each design variable i can be derived by the chain rule as Finally, the sensitivity of the total compliance of the body with respect to each design variable can be obtained by summing up the sensitivity of each element.

Alternative compliance computation and derivation
In the previous section, the compliance has been derived considering that the Gauss points used for density interpolation and element integration coincide. If a different number of Gauss points are used, the compliance can be derived as follows.
The compliance of each element of the body mesh can be computed as (Equation 2) The sensitivity of the compliance of each element with respect to each design variable i can be derived by the chain rule as where GP i is given in Equation (22) and the sensitivity of the total compliance of the body with respect to each design variable can be obtained by summing up the sensitivity of each element by Equation (23).
This alternative formulation for the sensitivity, even if more prone to lose some detail of the density field, can be used with any finite element without requiring special element formulation to perform the integration with a higher number of Gauss points. In fact, in standard elements, the number of Gauss points is usually small to reduce the computational time required for element integration and to avoid over-stiff elements. 30

NUMERICAL EXAMPLES
In this section, two numerical examples are considered to show the performances of the proposed interpolation scheme.
The first examples is the well-known MBB-beam problem. 31 In the second, the MBB-beam problem is modified to include a hole in the beam geometry.
In the examples, for brevity, the considered mesh elements will be referred according to the nomenclature reported in Table 1. In all the examples, the initial condition is a uniform density field with value equal to the desired volume fraction. This initial condition is the usual initial condition for this kind of problems. 2

MBB-beam
The MBB-beam problem 31 is depicted in Figure 3, where only half of the structure is considered for symmetry. The considered volume fraction for the problem is 50% and the material is steel (elastic modulus 210,000 MPa, Poisson coefficient 0.33). In Figure 4, two meshes are used for the solution of the problem. The body mesh (Figure 4 top) is realized by Q8 elements with characteristic dimension of 1 mm. The mesh strictly describes the geometry of the body. The interpolation mesh is realized by Q4 elements of characteristic dimension of 0.6 mm. In the example, as the interpolation mesh is not required to describe the geometry of the body (see Equations 16 and 17), the interpolation mesh is larger than the considered body. The solution of the MBB-problem 31 with the two meshes depicted in Figure 4 is reported in Figure 5. In Figure 5 top, the solution has been obtained by considering the compliance sensitivity formulation of Section 4.2, while in Figure 5  bottom the compliance sensitivity formulation of Section 4.3. In both cases, 16 Gauss points have been used for the integration of the density field. For the compliance sensitivity formulation of Section 4.2, 9 Gauss points have been used for the integration of the element stiffness matrix. The left column of Figure 5 shows the mean element densities of the body mesh, the central column the densities of the nodes of the interpolation mesh, and the right column the densities computed at the Gauss points of the body mesh. In both cases, no filter has been applied to densities or sensitivities. The initial and final compliance values computed for the two formulations are reported in Tables 2 and 3.
In both solutions of Figure 5, the final computed densities of the nodes of the interpolation mesh outside the domain of the body and not influencing the density of any Gauss point of the body mesh is equal to the minimum density value. These nodes do not influence the solution of the problem, however, the possibility to locate nodes outside the body allows for simpler and more regular interpolation meshes. The two solutions show a similar topology. The solution obtained with the formulation of the compliance sensitivity of Section 4.2 (i.e., the upper line of Figure 5) shows a more  Note: Compliance sensitivity formulation in Section 4.2. In all considered cases, the interpolation mesh has Q4 elements with characteristic dimension of 0.6 mm. a Mesh type definition in Table 1.

TA B L E 3
Results of the MBB-beam problem 31 of Figure 3 for different body meshes Note: Compliance sensitivity formulation in Section 4.3. In all considered cases, the interpolation mesh has Q4 elements with characteristic dimension of 0.6 mm. a Mesh type definition in Table 1. detailed structure. The formulation of the compliance sensitivity of Section 4.3, however, allows for the utilization of standard Q8 elements with 9 Gauss points, available in any finite element software. In both cases, sharp boundaries can be observed.
In Figures 6 and 7 and in Tables 2 and 3, the results of the MBB-beam problem 31 are reported for other body mesh types and the two compliance sensitivity formulations. In some cases, especially when linear or bilinear body meshes are considered, a certain amount of islanding and layering has been observed. In these cases, to remove these numerical problems, a sensitivity filter of the kind of the one described in Reference 2 has been applied. The filter radius required to

MBB problem with hole
In this example, a hole is inserted in the cantilever of the MBB-beam problem. 31 The resulting structure is reported in Figure 8. An example of the mesh used for the solution of this problem is given in Figure 9. In the interpolation mesh, the hole is not considered but the same simple rectangular domain employed to solve the MBB-beam problem of Sect. 5.1 is used. The problem is solved for a volume fraction of 30% and the material is steel (elastic modulus 210,000 MPa and Poisson coefficient 0.33). The problem is solved considering a body mesh of T6 elements of characteristic dimension of 1  Tables 4 and 5 for the two formulations.
The results show that the proposed approach can be applied to arbitrary complex body domain without requiring complex mesh transformation for the interpolation mesh. Actually, the interpolation mesh can be realized on simple and regular domain. In this way, the body can be meshed with the most convenient mesh for the displacement computation while a simple mesh for density computation can be used. In fact, the density of the nodes of the interpolation mesh located in the hole and not influencing the density of any Gauss point of the body mesh are set by the optimization algorithm to the minimum value of density, as their sensitivity to compliance is always zero (see Equations 18 and 24).
The compliance of the optimized structures is monotonically decreasing with the reduction of the element dimension of the interpolation mesh. As the interpolation mesh is refined, thinner and more defined features appear in the solutions. The compliance sensitivity formulation of Section 4.2 is able to get more effective structures than the compliance sensitivity formulation of Section 4.3. In particular, a similar structure is obtained by the two formulations with an interpolation mesh size of 0.6 mm for the former and 0.15 mm for the latter. The alternative formulation of Section 4.3 acts like a filter on the sensitivities and, in fact, no filter is required to avoid islanding even for interpolation mesh sizes for which the formulation of Section 4.2 requires a filter. Referring to the filter on sensitivities, the required filter size is much smaller than the body mesh element size. Moreover, the filter does not prevent the computation of more efficient structures with smaller features. The filter seems to have a limited effect on the obtainable topology. Finally, as the interpolation mesh is refined, the number of Gauss points required for the interpolation of the density field increases.

CONCLUSION
In this article, a novel density interpolation scheme suitable for nodal variable based structural optimization has been presented. The proposed scheme decouples the mesh required for the displacements calculation from the mesh describing the material distribution. The scheme is based on the interpolation of the nodal densities by a dedicated mesh of linear, bilinear (2D), or trilinear (3D) elements. The choice of these elements guaranties bounded values of interpolated densities between 0 TA B L E 5 Results of the MBB problem with hole problem of Figure 8  Note: Compliance sensitivity formulation in Section 4.3. In all considered cases, the body mesh has T6 elements with characteristic dimension of 1 mm and the interpolation mesh has Q4 elements. and 1 assuring a physical meaningful value to the interpolated values. The interpolation scheme allows explicit analytical formulations for the computation of the sensitivities that can be used in the optimization algorithm. In this article, two formulations have been proposed for the computation of the compliance sensitivity with respect to the nodal density design variables. The first formulation computes the sensitivity at each body element Gauss point. In this way, a more detailed sensitivity field can be obtained and more effective optimized topologies can be found. However, a very large number of Gauss points have to be used for the numerical integration of the stiffness matrix. As a consequence, special finite elements have to be used and an over-stiff stiffness matrix could be obtained. The second formulation considers a mean density for each element. Standard finite elements can be used, but a filtering effect on the sensitivity is found which reduces the freedom of the algorithm to find small features in the optimized structures. Numerical examples have been solved with the proposed interpolation scheme. The approach has shown to be rather stable, for the same interpolation mesh, to different kinds of elements mesh. As the interpolation mesh is refined, for the same body mesh, more effective optimized structures can be found. Numerical tests have shown a certain tendency of the method to islanding instabilities. The solution is more prone to numerical instabilities if linear elements are used for the body mesh and as the dimension of the interpolation mesh is reduced with respect to the body mesh dimension. A sensitivity filter can be effectively used to prevent such instabilities. The required radius of the filter is small, of the order of about half of the dimension of the body mesh elements. Also, the filter seems to have a very limited effect on the optimized structure topology.
The proposed approach can be easily applied to any body geometry. Also, the interpolation mesh can cover a larger domain with respect to the body mesh. In this way, very regular interpolation meshes can be obtained. Further developments can refer to the inclusion of adaptive body mesh in the optimization algorithm as the proposed approach has very interesting properties of volume and material volume conservation.