Level set topology optimization with nodally integrated reproducing kernel particle method

A level set topology optimization (LSTO) using the stabilized nodally integrated reproducing kernel particle method (RKPM) to solve the governing equations is introduced in this paper. This methodology allows for an exact geometry description of a structure at each iteration without remeshing and without any interpolation scheme. Moreover, useful characteristics of the RKPM such as the easily controlled order of continuity and the ability to freely place particles in a design domain wherever needed are illustrated through stress based and design-dependent surface loading examples. The numerical results illustrate the effectiveness and robustness of the methodology with good optimization convergence behavior and ability to handle large topological changes. Furthermore, it is shown that different particle distributions can be used to increase efficiency without additional complexity.


Introduction
The vast majority of topology optimization methods are based on the finite element method (FEM).In recent years, authors started exploring meshfree methods in topology optimization and take an advantage of their useful features such as mesh independency, flexible control of approximation and smoothness and convenience in performing adaptive local refinement.Indeed, in the last decade meshfree methods have been shown to be useful in a wide range of problems.Examples include nonlinear fracture mechanics [1], large deformation problems [2] and contact mechanics [3].Over the years, a wide variety of meshfree methods have been studied and they can generally be categorized into collocation and Galerkin meshfree methods.Collocation methods are based on the strong form of partial differential equations (PDEs).Galerkin meshfree methods on the other hand, are based on the weak form of PDEs.The element free Galerkin method (EFG) [1] and the reproducing kernel particle method (RKPM) [4,5] are typical Galerkin meshfree methods.A recent review on meshfree methods can be found in [6].
A significant part of the literature on meshfree topology optimization combines meshfree methods with densitybased approaches such as the Solid Isotropic Material with Penalization (SIMP) method.Especially the use of Galerkin based methods has been shown to alleviate the need for sensitivity filtering due to the higher order smoothness of the meshfree shape functions and avoid mesh-dependence phenomena [7,8].Commonly, densitybased meshfree methods are based on point-wise density interpolation schemes, with either the densities at the Gauss points considered directly as design variables or with the nodal densities defined as design variables and then used to interpolate the density field at the computational points based on the meshfree shape functions as explained in the next paragraph.
Du et al. [9] used the EFG method integrated into topology optimization with the SIMP method for the design of thermomechanical compliant mechanisms with geometrical nonlinearities.To eliminate the appearance of discontinuous scattered points in the topology, sensitivity filtering is used to smear out the numerical instability.For the integration of the weak form Gaussian quadrature on rectangular background cells is used, and the artificial densities at the Gauss points are considered as the design variables.The EFG method with Gaussian integration on a background mesh in combination with the SIMP method was subsequently used by other authors to study a variety of problems.Gong et al. [7] presented a study with multiple loading conditions and compliance minimization under stress constraints.The densities at the particles are chosen as the design variables and are used to interpolate the densities at the computational Gauss points using the meshfree shape functions.This exploits the smoothness of the meshfree shape functions to improve the smoothness of the density field and eliminate the checkerboard problem.Zheng et al. [10] studied free vibrating continuum structures with the objective of maximizing the fundamental eigenvalue and considering the nodal density as design variable and [11] applied the EFG-SIMP method to the modal topology optimization method for maximizing the first-order natural frequency.Luo et al. [8] proposed an approach in which the Shepard function is applied to construct a physically meaningful dual-level density approximant, due to its non-negative and range-restricted properties.The density at any computational point is interpolated by the Shephard function and the density values of nodes inside the influence domain of the current point.This density approximant acts as a heuristic filtering scheme to enhance the smoothness of the density field.Geometrically nonlinear structures were considered by He et al. [12].The approach is based on a similar point-wise material density field and Shepard interpolation used in [8].The major difference in this work is that the density design variable points are freely positioned independently of the displacement nodes used in the displacement analysis.However for ensuring well-posedness of the optimization problem, the density points should be distributed in a reasonable manner according to the displacement node arrangement.Gong et al. [13] proposed a SIMP based topology optimization-EFG method based on moving particles rather than "soft delete" or "hard delete".The moving of particles is controlled based on the density values change in each optimization step, by moving particles with densities lower than a threshold value to positions with higher density values and keeping high density particles still.Zhang et al. [14] used the EFG method combined with the rational approximation of material properties (RAMP) model for topology optimization of heat conduction in both isotropic and anisotropic materials.The densities of the EFG nodes are chosen as design variables and the densities at the Gauss points are computed by the interpolation of the nodal densities with the moving least squares (MLS) shape function in the domain of influence.
There are also implementations based on different schemes in terms of the meshfree method and integration scheme.Cho and Kwak [15] used the reproducing kernel particle method in combination with a density-based approach for the topology optimization of geometrically non-linear structures.The bulk density at the integration points of the background mesh is interpolated using the nodal densities with support domains covering each integration point using the RK shape function.Du et al. [16] used the EFG method for multi-material topology optimization based on the alternating active-phase algorithm and within the SIMP framework.The nodal relative density which is taken as the design variable is obtained through the Shepard interpolation in combination with the MLS shape function.In contrast with most meshfree topology optimization methods presented in the literature, this work employs nodal integration instead of Gauss quadrature which increases the computational efficiency.Due to the order smoothness of the MLS shape functions sensitivity filtering is not used.A direct coupling between the FEM and EFG was presented by Zhang et al. [17] to reduce the computational cost of the EFG method.A constraint centroidal Voronoi tessellation (CCVT) algorithm linked with the density variable of SIMP is used to generate the point set for discretizing the EFG domain.A dual-level density approximant is adopted to formulate a continuous material distribution.The density field at the Gauss points is interpolated by the density field at the FE and EFG nodes in the influence domain using two different Shepard functions.The domain integration is performed using the Gauss quadrature on the background cells whereas multilevel Gauss quadrature points are used to further reduce the computational cost.The authors later extended the approach to topology optimization of hyperelastic structures [18].Zhou and Zou [19] presented a meshless topology optimization method based on the implicit topology description and the reproducing kernel particle method using nodal design variables and the smoothed Heaviside function with a regular background mesh for the integration of the weak form in the domain.
Some works exist in the literature working with the smoothed particle hydrodynamics (SPH) method.Lin et al. [20] introduced the corrective smoothed particle method and total Lagrangian formulation to eliminate intrinsic problems in the (SPH) such as inconsistency and instability.Each SPH particle is assigned a density design variable and topology optimization is performed using a modified SIMP approach.In contrast to the Galerkin meshfree topology optimization methods, this approach requires a sensitivity filter to avoid the formation of checker-board patterns.Minimization of compliance was considered in the examples.The corrective smoothed particle method (CSPM) was also used by Li et al. [21] for optimizing linear structures under multiple load cases.Particle's density is chosen as the design variable and it is approximated using the Shepard interpolation scheme to make it smooth and range-restricted in [0, 1].The modified SIMP method is applied to determine the Young's modulus, which is related to the particle densities.A sensitivity filtering technique is also introduced to enhance the numerical stability and prevent the checker-board patterns.
In the context of level set methods, Luo et al. [22] proposed a meshless Galerkin parametric level set method based on compactly supported basis functions (CSRBFs).CSRBFs are used to both parameterize the level set function and to construct the shape functions for the meshfree approximation.For the domain integration a fixed background mesh is employed with 4 × 4 Gaussian quadrature.The particles were placed at the nodal positions of the background mesh and their positions remain unchanged throughout the optimization process.To account for the boundary discontinuity, the CSRBFs were used to interpolate the level set function at quadrature points based on the discrete level set function nodal values.Quadrature points with a positive or zero level set value were considered as solid whereas a weak material was assigned to points with negative level set values.Compliance minimization for linear elastic structures under constant point loads was solved in this work and the authors later extended methodology to the optimization of compliance multiphysics actuators [23].Subsequent level set works based on meshless methods also used interpolation schemes and Gaussian integrations on fixed background meshes.For example the parametric level set method based on CSRBFs was also used by Ai and Gao [24] to optimize twodimensional microarchitected periodic metamaterials for the cases of single and multiple materials.Similarly to Luo et al. [22], the material discontinuity is treated by an interpolation scheme in which the level set function value at a Gaussian quadrature point is determined through the interpolation of the field nodes in the support domain using the meshfree shape functions.Khan et al. [25] combined the EFG method with the level set topology optimization using the Hamilton-Jacobi equation with a topological derivative term.Two-dimensional linear elastic problems under single and multiple point loads were examined considering compliance minimization as the objective function.The particles were at fixed locations throughout the whole optimization process.Recently, Neofytou et al. [26] employed the reproducing kernel particle method in combination with level set topology optimization to solve minimum compliance problems with design-dependent pressure loads.The capability of RKPM to freely place particles along the structural boundary is exploited in this work to apply the moving loads directly onto the structural boundary, avoiding the need to use any load interpolation scheme.The particle distribution coincides with the nodal positions of the background mesh plus extra particles placed along the boundary defined by the level set method.An interpolation scheme based on the RK shape functions and level set values at the particles is used to treat the boundary discontinuity during the domain integration.
The bi-directional evolutionary structural optimization (BESO) method was also used in combination with the EFG method by Shobeiri [27] employing uniformly distributed nodes and Gaussian quadrature in the examples for minimum compliance of linear elastic structures.Zhao [28] incorporated an improved meshless density variable approximation into the BESO method and used the Shepard function to create a physically meaningful dual-level density approximation based on the previous work by Luo et al. [22].The shape functions of the meshless Galerkin method are constructed using CSRBFs and domain integration is performed on a regular background cell structure using high order Gauss quadrature.
As seen from the literature, previous works are based on interpolation schemes based on nodal densities.Even level set methods use such interpolation schemes thus not transferring the clear boundary representation generated by the level set method on the computational domain.Conventionally, in level set topology optimization fixed grids and interpolation schemes are used to avoid the cumbersome task of remeshing the structure at each iteration.As an alternative, immersed methods such as the extended finite element method (XFEM) [29] and the finite cell method (FCM) [30] can be used to maintain the crisp boundary.The XFEM achieves this by adding enrichment functions and additional degrees of freedom to nodes around the boundary discontinuity.However, several difficulties can occur due to the cut elements.For example, noise and ill-conditioning of the discretization may occur in case of small intersections of finite elements (i.e., cut elements with small solid fractions) [31].Another potential source of error is the situation of non-physical coupling when the gap between solid components is less than the size of an element [32].The FCM can be seen as an extension of the classical finite element method that uses higherorder non-boundary-fitted meshes.An adaptive quadrature scheme is used by decomposing the cut elements into sub cells which may lead to an increased number of quadrature points, since full Gauss quadrature is employed in each sub-cell [33].On the other hand, works that employ LSTO with remeshing have also been presented in the literature.For example, the work by Xia et al. [34] uses such a remeshing approach for stress-based topology optimization.In this work, a conforming triangular mesh is created at each iteration.This is done by modifying an initial background mesh such that intersected elements are first clipped to conform to the boundary, and subsequently adjusted to improve the quality of the mesh.Then only solid elements are considered for the analysis.Later the authors extended this work by converting the triangular mesh to a quadrilateral body-fitted mesh for improved accuracy [35].Allaire et al. [36] also presented such an approach based on mesh evolution and demonstrated its applicability in a range of different physics such as coupled thermal fluid-structure interactions [37].Nevertheless, a series of operations have to be performed at each optimization step to ensure mesh quality.One of the main benefits of the meshfree methods is the ability to freely place particles in the design domain without relying on a mesh and without the problems caused by mesh dependency.This can allow for an exact representation of shapes without remeshing and the requirement of element quality operations.However, meshfree topology optimization approaches that use interpolation schemes do not fully exploit the advantages of the meshfree methods and lose the potential of an exact representation of the structural geometry.Furthermore, most Galerkin meshfree methods used in the literature use Gauss integration on a background mesh with higher order quadrature.This is necessary to ensure accuracy, as it was identified in [1] as a requirement due to misalignment of the support with the background cells.
In this work we combine a level set topology optimization method with the reproducing kernel particle method with nodal integration for the domain.The RKPM is used here to solve the governing equations, then the level set function is used for topology optimization.The aim is to have an exact geometry description of the structure without remeshing.This is achieved by placing particles within the structure and on the structural boundaries.This way, the clear boundary generated by the level set method is maintained on the computational mesh with an exact description of shapes without any interpolation scheme, enjoying the full advantages of the meshfree method.Moreover, several useful features of RKPM are investigated such as the ability to have different orders of continuity by choosing different kernel functions, the freedom of having different particle distributions, and the efficiency due to the use of nodal integration compared to Gauss quadrature approaches.
The remainder of the paper is organized as follows: Section 2 outlines the level set topology optimization method.Section 3 describes RKPM and its main ingredients such as construction of RK shape functions, Galerkin formulation and the naturally stabilized nodal integration.Section 4 describes the implementation of the methodology.In Section 5 stress-based and design-depended numerical examples are presented.Section 6 offers some concluding remarks.

Level set topology optimization method
This section briefly summarizes the level set topology optimization method used in this study.More details of the method can be found in [38] and [39].In the level set topology optimization, the structural boundary is defined as the zero level set of an implicit function: where φ is the level set function, Ω is the structural domain and Γ is the structural boundary.Conventionally, the implicit function is initialized as a signed distance function [40,41].
The structural boundary is optimized by iteratively solving the following Hamilton-Jacobi equation ∂φ(x, t) ∂t where t is a fictitious time domain for the level set evolution and V n is the normal velocity.
The level set function at each node is updated by solving the following discretized Hamilton-Jacobi equation using up-wind differential scheme, where r is a discrete node in the design domain, V nr is the normal velocity at node r , k is the iteration number and ⏐ ⏐ ∇φ k r ⏐ ⏐ is computed for each node using the Hamilton-Jacobi weighted essentially non-oscillatory method (HJ-WENO).To improve the computational efficiency, the level set update is restricted to nodes within a narrow band close to the boundary.This means that φ r is given by a signed distance to the boundary only within the narrow band.To correct this effect, φ r is periodically reinitialized to a signed distance function.For the reinitialization and velocity extension the fast-marching method is used [42].The velocities required for the level set update are obtained by solving the linearized optimization problem, where f is the objective function, g m is the mth inequality constraint function, ∆Ω k is the update for the design domain Ω and g k m is the change in the mth constraint.Shape derivatives that provide information about how the objective and constraint functions change with respect to the boundary movement, typically take the form of boundary integrals [41]: where s f and s gm are the shape sensitivity functions for the objective and the mth constraint, respectively.The integrals in Eqs. ( 5) and ( 6) can be estimated as, where j is a discrete boundary point, V n j , s f, j and s gm, j are the normal velocity and sensitivities for the objective and mth constraint functions, respectively, at point j.l j is the length of the local boundary around the boundary point j, C f and C g m are vectors containing the product of boundary lengths and shape sensitivities and V n is the vector of normal velocities.
The linearized optimization problem is solved using the Simplex method [43].The obtained optimum boundary point velocities are then used in Eq. ( 3) to update the level set function and this process is repeated until convergence is obtained.This method is implemented in the object oriented C++ code [44] and is available as opensource at http://m2do.ucsd.edu/software/.Further explanatory details on the implementation of the LSTO method described here can also be found in [45].

Reproducing kernel particle method (RKPM)
This section outlines the reproducing kernel particle method (RKPM) and the main ingredients of the methodology used here.For more details the readers are referred to [46].

Reproducing kernel approximation
To construct the RK approximation for a finite dimensional solution of the PDEs, the domain Ω is discretized by a set of nodes {x 1 , x 2 , . . ., x N P }, where x I is the position vector of node I , and N P is the total number of particles, then the RK approximation of a function u(x), denoted by u h (x), is expressed as, where Ψ I (x) is the RK shape function at node I , and u I is the corresponding nodal coefficients to be determined.Then the RK shape function can be expressed as where the basis vector H(x − x I ) is defined as where n is the specified order of completeness and the moment matrix M(x) is The kernel function Φ a (x − x I ) centered at x I with compact support size a controls the smoothness (continuity) and locality of the approximation function.The order of continuity in this approximation can be introduced without complexity.For example, the box function gives C −1 continuity, the hat function leads to C 0 continuity, the quadratic spline function results in C 1 continuity, the cubic B-spline function yields C 2 continuity and the quartic and quintic functions give C 3 and C 4 continuity, respectively.This is illustrated in Fig. 1 where the C 0 linear B-spline (tent) kernel is compared with the C 2 cubic B-spline kernel function: Linear B-spline (tent), Cubic B-spline, where z I is defined as: The support size a I is typically defined as: where c is the normalized support size, chosen between 1.5 and 2.0 in practice, and h I is the nodal spacing associated with point x I [46].The scheme for defining the nodal spacing in this work is explained in Section 4.2.

Galerkin formulation
In this work we consider the case for linear elasticity where u i is the displacement, σ i j is the Cauchy stress, n j is the surface normal, b i is the body force and h i and g i are the prescribed traction and displacement on the Neumann and Dirichlet boundaries, Γ h i and Γ g i , respectively.The RK shape functions lack the Kronecker delta property and thus applying the essential boundary conditions is not as straightforward as in FEM.Several approaches have been proposed for this purpose.In this paper, Nitsche's method [47] is used which results in the following matrix equations where and for the two-dimensional case where K is the stiffness matrix, B is the gradient matrix, C is the elasticity tensor, β is the penalty taken as β = β nor E h with β nor the normalized penalty parameter, E the Young's modulus and S is used to impose each component of the boundary displacement by setting s i = 0 or 1.

Naturally stabilized nodal integration (NSNI)
Galerkin meshfree methods such as RKPM are based on the weak form of PDEs.Thus, domain integration is required to evaluate the integrals in the weak form.This can be done either by performing Gaussian integration with background cells as the quadrature domains or by using nodal integration employing nodal representative domains.In this work we employ the naturally stabilized nodal integration (NSNI) technique proposed in [48] which presents good accuracy and stability while maintaining the efficiency of nodal integration.An outline of the method is given in this section whereas for a more in-depth discussion the reader is referred to [46].
To employ the NSNI, a smoothed gradient ΨI,i of the shape functions at each nodal point is computed using the divergence theorem as follows, where A N is the area of the nodal representative domain Ω N associated with node N , Γ N is the boundary of the nodal representative domain, and n i denotes the ith component of the outward unit normal vector to the representative domain.The nodal representative domains can be constructed by Delaunay triangulation or Voronoi diagram as shown in Fig. 2. For the evaluation of ΨI,i at the nodal points using Eq. ( 23), boundary integration of the nodal representative domain is needed.This is done in this work using a one-point Gauss integration rule [49] where S N is the number of midpoints for each boundary segment of the Voronoi cell associated with node x N and L k is the length of the k th segment.
Finally the smoothed gradient matrix BI for the two-dimensional case can be expressed as The smoothed gradient matrix B can now replace matrix B in Eqs.(19) and (20).It is important to emphasize that although the strain smoothing involves evaluation of shape functions at the edges of the representative nodal domain, it is only needed for shape functions.The stiffness is integrated entirely at the particles.Also, since the smoothed gradient matrices are consistent with the weights in the nodal integration, the distribution of weights, i.e. the areas of the cells, have marginal effect on the solution accuracy.
The NSNI introduces an implicit gradient expansion of the strain field based on the first order Taylor expansion of strain around x N , to avoid instabilities associated with conventional nodal integration techniques.The implicit gradient is expressed as where for a linear basis the vector H i takes the following forms The stiffness matrix can be expressed as where B ∇ I 1 (x N ) and B ∇ I 2 (x N ) are defined as with Ψ ∇ I i, j (x N ) obtained by direct differentiation of Ψ ∇ I i in Eq. ( 26) with respect to x j [48].The terms M 1 (x N ) and M 2 (x N ) are the second moments of inertia in each nodal integration domain given by Stabilization with Eq. ( 28) is termed naturally stabilized nodal integration (NSNI) since the constants M 1 (x N ) and M 2 (x N ) associated with the additional terms occur completely naturally, and thus no tuning of any parameters is required, which is in contrast to other stabilized methods [50,51].

Level set method and nodally integrated RKPM
The nodally integrated RKPM is used in this work within the context of level set topology optimization to achieve an exact description of structural geometry without re-meshing.The level set boundary and signed distance values are used for the particle distribution and the construction of the Voronoi diagram of the nodally integrated RKPM method as explained in this section.In turn, the LSTO benefits from the useful properties of the RKPM method such as placing particles on the discretized boundary points to compute sensitivities, and the higher order basis functions to improve accuracy.This section describes several ingredients of the methodology which are further demonstrated through the examples in Section 5.

Particle positions and boundary sensitivity
As explained in Section 2, to obtain the optimum velocities required for the update of the level set function and thus the topology, shape sensitivities need to be computed at the discretized boundary points i.e. the points at which the boundary intersects the elements of the level set mesh.The positioning of particles at these boundary points are shown in Fig. 3, allows for direct computation of sensitivities.Thus, at each optimization iteration, particles are always placed at the boundary points first.These boundary particles are stored at the beginning of the list of all particles to separate them from the particles in the interior of the structure.The particles in the interior of the structure can be freely positioned.For example in Fig. 4(a) the particles are regularly positioned at the nodal positions of the level set mesh with the respective Voronoi diagram for the NSNI shown in Fig. 4(b).In Fig. 4(c) the particles are randomly positioned and the resulting voronoi diagram is shown in Fig. 4(d).Based on the above discussion on boundary sensitivity, accuracy is more critical close to the boundary.In Section 5.1.1 a particle placement scheme where more particles are strategically placed within a region close to the boundary is discussed and illustrated through a stress minimization example.

Construction of voronoi diagram
For the construction of the nodal representative domains used in the NSNI, the Voronoi diagram is employed here using the open source software library voro++ [52].For the construction of the diagram, the particle coordinates are given as an input which results in an unclipped Voronoi diagram.An example is shown in Fig. 5. Fig. 5(a) shows the unclipped diagram for a rectangular plate with a hole, with the cells along the boundary extending outside of the structure and get intersected by the boundary segments.In order to make the diagram conforming to the structure, which is employed for the NSNI method, the intersected cells need to be clipped as shown in Fig. 5(b).For such clipping process, the signed distance values at the vertices of each intersected cell are employed as described in Fig. 6.The intersections of level set boundary segments with the edges of the Voronoi cells are first computed.For the most part, the cells of boundary particles will be intersected by the boundary segments, however it can happen that interior cells are intersected as well.If an intersected cell is also a boundary cell, the intersections along with the particle to which the cell belongs split the cell into two polygons as shown in Fig. 6(a).When interior cells are intersected, only the intersection points split the cell into polygons as shown in Fig. 6(c).The polygon with the biggest sum of the signed distance values at its vertices is kept as the new clipped cell, while the vertices of the remaining polygon are being removed.This is illustrated in Figs. 6 (b) and (c) for boundary and interior cells, respectively.The remaining polygon after clipping, for example polygon A in Figs.6(b) and (c), is then used as the integral domain to calculate the stiffness matrix.Since this clipping operation only has to be done for the cells along the boundary and maybe a few neighboring interior cells, the process is not computationally expensive.

Size of support domain
The voronoi diagram information can also be used to define the support domain size for each particle in a way that it is consistent with the density of the particle distribution around it, as explained in [53].The neighbors of each particle are identified as those that share an edge with the particle's cell.For example in Fig. 7 particles 2-11 are neighbors of particle 1.The support size a I of particle 1 is then defined as a = c • d max (31) where d max is the distance to the furthest neighbor as shown in Fig. 7 in which the distance between particles 1 and 11 is the maximum amongst all neighbors.The normalized support size c is chosen as 2.0, which is twice the order of the linear basis used in this work.As can be seen in the figure, in areas with denser particle distribution support sizes are smaller than areas with less dense distribution.

Examples
In the literature of Galerkin meshfree methods, stress-based computations have been used to validate against the FEM in terms of accuracy and convergence properties [1,54,55] for linear elastic structures.As a simple illustration, the example of a linear elastic plate with a hole is presented here as shown in Fig. 8(a).Specifically, the FEA is carried out via the partial differential equation toolbox in MATLAB, using linear, triangular elements.The FEA mesh is shown in Fig. 8(b).For the RKPM analysis, the mesh generated in MATLAB is used to place the particles at the nodal positions of the elements in order to create a comparable discretization.The particles along with the voronoi diagram are shown in Fig. 8(c).Symmetry conditions are applied at the left and bottom edges.The plate is loaded uniaxially and analyzed under plane stress conditions.The values used for Young's modulus and Poisson's ratio are 1 and 0.3, respectively.The von Mises stress is computed by the two methods at the nodes on the circular boundary to compare the accuracy of each approach against the analytical solution [1].As can be seen in Fig. 8(d), the solution obtained by RKPM is in good agreement with the analytical solution for this specific discretization, which consists of 230 degrees of freedom (DOF).Of course, when a finer FEA mesh is used as shown in Fig. 8(d) for 5114 DoF, FEA also provides an accurate solution as expected.In the context of topology optimization, stress-based problems are challenging and remain an active area of research [39,[56][57][58].For the aforementioned reasons, stress-based problems are thus chosen here to illustrate the effectiveness of the proposed RKPM-LSTO methodology and to investigate several interesting features such as different particle distributions, kernel functions and efficiency.A design-dependent problem with hydrostatic pressure loads is also used to show how the crisp representation of the level set boundary in combination with particle placement at boundary points can straightforwardly handle the design dependency without interpolation for either the loads or the domain integration.For the examples Young's modulus value is set to 1, Poisson's ratio is 0.3, the dimensionless support size c is 2.0 and the normalized penalty parameter β nor for Nitsche's method is chosen as 1000.Convergence of the objective is checked during 5 consecutive iterations and the tolerance is 0.001.

Stress-based examples
The p-norm function is used here which approximates the maximum stress in the structure.This is a global stress measure and can be used either as an objective for stress minimization, or as an inequality constraint.The p-norm functional defined in the domain Ω can be written as follows where p is the p-norm parameter, σ P N denotes the p-norm stress and σ υm is the von Mises stress.
For the solution of the stress-based problems the shape sensitivity, G ′ (Ω ), of the p-norm functional G(Ω ) is used, the derivation of which can be found in [39], where C is the elasticity tensor, ε represents the mechanical strain and u and λ denote the state and adjoint variable vectors respectively.Eq. ( 33) indicates how the p-norm stress functional changes when a structural boundary point moves in its normal direction.The shape sensitivity terms are computed directly at the boundary points where RKPM particles are placed.

Example 1: Stress minimization
Stress minimization under a volume constraint for the L-bracket shown in Fig. 9(a) is considered as a first example.The optimization problem is expressed as follows where V is the limit for the volume constraint set as 70% and the p-norm value used is 6.Starting from a design without any holes as shown in Fig. 9(b), the optimum hook-like structure of Fig. 9(c) is obtained as expected.For this result the cubic kernel spline was used for the construction of RK shape functions.In the following paragraphs we make use of this example to investigate the effects of different features such as particle distribution and kernel function choice.

Regular particle distribution
Having a regular particle distribution is useful when trying to reduce the computational cost during the optimization procedure.This is because with such a distribution, new particles are only inserted along the boundary at each iteration.The supports of these newly inserted particles cover the particles in the region around them and   thus, shape functions of the particles only in that region need to be corrected to ensure polynomial reproductivity.The particles further away from the boundary, i.e. greater than a certain value of the signed distance function, do not see any change and thus their information can be stored before the optimization starts and then re-used throughout the optimization process.This is illustrated in Fig. 10.In Fig. 10(a) the information for particle "P" including support domain and stiffness matrix are computed on the regular grid.In Fig. 10(b) after new particles have been inserted on the boundary created by the circular hole, the information for particles within the highlighted region has to be computed at this step, whereas since particle "P" is unaffected its information can be read from the previous step.Support sizes vary depending on neighbor distances.However, since particles are irregular only up to the neighbors of boundary points, this variation is limited compared to support sizes of the regularly distributed particles in the interior region.This means that the size of the zone within which the shape functions need to be re-computed does not have to be much larger than the size of the average support size.A size of twice the size of a support of an interior particle is a sufficiently large value.This size makes sure that prestored information is only used for particles that are not covered by new particles.The scheme described above is shown for the stress minimization of the L-Bracket example in Fig. 11.The dark highlighted zone indicates particles for which the information has to be re-computed.The light highlighted regions indicate particles for which the prestored information can be used.In terms of computational efficiency, we compare the overall time for the whole optimization process with and without the prestoring scheme for the same number of iterations.The prestoring scheme turns out to be 1.4 times faster than the normal computation ( without pr estor e with pr estor e ≈ 1.4).The overall assembly time for the stiffness matrix for the entire  optimization process is two times faster for the prestoring scheme.It is of course expected that this difference will be more significant as the number of particles increase.

Different choices of kernel functions
RKPM provides controllable orders of continuity and completeness, independent from one another, without additional complexity just by choosing the kernel function Φ a and the basis H T (x), respectively, in Eq. (10).We experiment here with three different kernel functions to see the effect on the convergence behavior.Specifically we solve the L-Bracket example with the linear (C 0 continuity) and cubic (C 2 continuity) kernel functions given in Section 3.1 and the quintic kernel (C 4 continuity) given as, The optimization is left to run for 200 iterations.Figs.12(a), (b) and (c) show the history p−norm stress for the linear, cubic and quintic kernel functions respectively.Figs.12(d), (e) and (f) present the final topologies for the three cases.For the 1st and the 3rd order kernels although the final structure is reached around iteration 140, the 0.001 tolerance is not achieved before 200 iterations whereas with the 5th order kernel the tolerance is achieved at 148 iterations as shown in Fig. 12.The final hook-like optimum solutions are obtained in all three cases, although boundary oscillations are visible in the case of linear kernel function as shown in Fig. 12(d).As the continuity order of the kernel function increases, the smoothness of the RK shape function increases accordingly.As can be seen from these results increasing the continuity order of the kernel can make the optimization process smoother and speed up the convergence.

Quadtree particle distribution
As explained in Section 2, the level set function is updated based on shape sensitivities computed at the boundary points.Good accuracy is thus more important when close to the boundary than in the interior of the structure.The flexibility that RKPM offers in adaptivity can be used to create a denser particle distribution within a narrow band close to the boundary whereas keeping the distribution coarser away from the boundary at each iteration.Such a scheme is shown in Fig. 13 for the L-Bracket example where a quadtree structure is used to create the particle distribution.Specifically the particles away from the boundary are kept regular and coarse, while within a narrow band close to the boundary the distribution gets denser according to the quadtree data structure.As can be seen in the figure, the same end result is obtained as with the fully dense particle arrangement in Fig. 9, only here much less particles are used.In the case of the fully dense particle distribution, the optimization starts with 6601 particles for the first iteration and ends with 5090 particles for the last iteration, whereas the quadtree distribution starts with 1838 particles and ends with 1542 particles.Regarding computational efficiency, the overall time for the whole optimization process is 1.5 times faster for the quadtree distribution compared to the regular dense distribution (i.e.r egular distribution time quadtr ee distribution time ≈ 1.5) which again will scale up as the mesh size increases.The true benefits of this scheme are expected to be more significant for larger scale and 3-dimensional structures and this will be the aim of our future work.Moreover, for such larger scale problems the quadtree distribution can be combined with the prestoring scheme for the regularly spaced interior particles to achieve even better efficiency.The time spent for the construction of the Voronoi diagram per iteration is shown in Fig. 14 for both regular and quadtree particle distributions.The time comparison is done first at the initial iteration when the number of particles is maximum and thus the process will take the most time.The comparison is repeated for iteration 120 at which point the particle number is decreased due to the removal of material.As can be seen, the construction of the voronoi diagram takes about 11% of the total time per iteration for the regular distribution at the first iteration.For the quadtree distribution the process is faster due to the smaller number of particles, taking about 4.5% of the total time.As expected, the process becomes faster for less particles at iteration 120 in both cases.Alternative implementations for the construction of the voronoi diagram can of course further improve the computational efficiency.For example, since the positions of the particles in the interior of the domain do not change as the topology evolves, only partial reconstruction of the voronoi diagram close to the boundary is required.Such an efficient implementation will be further examined in the future.It is important however, to note that the main difference between the voronoi diagram used here and remeshing, is that the RKPM with NSNI is not sensitive to the shape of the polygons, as long as the polygons are conforming.Once the diagram is constructed, no additional operations are required to modify the shape of the polygons, whereas for FEA with remeshing additional operations are performed to improve the mesh quality.

Example 2: Volume minimization and stress constraint
In this next example volume minimization under stress constraint is considered.As explained in [39] the p-norm stress is always greater than the actual maximum.Thus, an adaptive scaling scheme is used, where σ is the stress constraint limit and c is defined as where η ∈ (0, 1] controls the variations between c k and c k−1 .The scaling uses information from the previous optimization iteration to normalize the constraint as c k σ P N so it approximates the actual maximum stress σ max as the design converges. The optimization problem is formulated as, minimize ∫ Ω dΩ subject to c k σ P N ≤ σ (38) The cubic kernel spline is used in this example and the evolution of the structure throughout optimization is shown in Fig. 15 Starting from an initial design with holes, optimization converges at 178 iterations.The convergence history for σ P N , σ max and the volume is shown in Fig. 16.This example also aims to illustrate the robustness of the methodology in the presence of holes and large topological changes.pressure load is applied and p is the pressure load.The pressure load is assumed to be constant although the method can be generalized for varying pressure.
where p 0 is the constant magnitude of the pressure load, and n is the surface normal.Shape sensitivity for the structural compliance function when the surface load is a pressure load was derived by [41] as, where p 0 is the pressure load and V n is the normal velocity on the boundary.The result obtained here as shown in Fig. 18 is the same as the one in [26] and the solution is achieved in about the same number of iterations.However now the process is 20 times faster.This is because in the case of the Gaussian integration approach almost 3×10 5 Gauss points were used to achieve the sufficient accuracy (16 Gauss points per cell), whereas the nodally integrated RKPM achieves the same level of accuracy with only 1800 particles.Fig. 18 shows the evolution of the structure whereas Fig. 19 shows the Voronoi diagram at the initial and final iterations.As it can be seen the methodology has no problem handling the large topological changes of this example such as holes merging and islands disappearing.

Conclusion
In this paper a combination of the nodally integrated RKPM with the LSTO method has been proposed for an exact geometry description of structures during optimization without any interpolation schemes and without remeshing.The effectiveness and robustness of the methodology has been illustrated through stress-based examples and a design-dependent problem with hydrostatic pressure loads.As shown, the method presents useful features such as the ability to change the order of continuity simply by changing the continuity order of the kernel function without adding computational complexity.Increasing the order of continuity in the RK shape function has been shown to speed up the convergence.Moreover, the method is efficient in handling different particles distributions which can be used to increase efficiency by increasing the particle density only around the boundaries where the sensitivity accuracy is more crucial.Furthermore, large topological changes such as holes merging and disappearing can be handled naturally.As meshfree methods have been shown to perform well in problems where the FEM struggles, such as cases with large structural deformations [6] and contact mechanics [3], we plan to explore these classes of problems for topology optimization.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 3 .
Fig. 3. Particles are first placed at the discretized boundary points.

Fig. 6 .
Fig. 6.Clipping process for intersected cells for: (a) Boundary intersected cells and (b) interior intersected cells.First, an intersected cell is split into two polygons based on intersection points with the boundary segments.The next step is clipping the cell based on signed distance values at the vertices.The sum φ 1 + φ 2 + φ 5 is larger than φ 3 + φ 4 and thus polygon B is removed.

Fig. 7 .
Fig. 7. Support domain size definition based on Voronoi diagram information.

Fig. 8 .
Fig. 8. Plate with a hole example: (a) Problem setup, (b) FEA mesh, (c) RKPM particles and Voronoi diagram and (d) von Mises stress computation at the nodes along the circular hole by FEA and RKPM, and comparison with the analytical solution.

Fig. 9 .
Fig. 9. L-bracket example: (a) problem definition, (b) von Mises stress field for the initial design and (c) von Mises stress field for the optimum solution.

Fig. 10 .
Fig. 10.Prestoring scheme: (a) Supports and local stiffness matrices are computed and stored for a regular particle arrangement before optimization starts and (b) Once new boundary particles are created, information only has to be computed for particles in the highlighted zone.Particles such as particle P do not see any change and thus the information can be used from the initial regular distribution.

Fig. 11 .
Fig. 11.Prestoring scheme for the L-Bracket example: The darkly highlighted zone indicates particles for which the information has to be re-computed.Lightly highlighted regions indicate particles for which prestored information can be used.

Fig. 12 .
Fig. 12. Convergence history of the minimum stress L-bracket design for different kernel functions: (a) Linear B-spline (C 0 continuity), (b) Cubic B-spline (C 2 continuity) and (c) Quintic B-spline (C 4 continuity), (d) optimum solution for linear kernel function, (e) optimum solution for cubic kernel function and (f) optimum solution for quintic kernel function.

Fig. 13 .
Fig. 13.Stress minimization for the L-bracket example with quadtree particle distribution: (a) Initial design with von Mises stress field, (b) optimum solution with von Mises stress field, (c) Voronoi diagram for initial design and (b) Voronoi diagram for optimum solution.

Fig. 14 .
Fig. 14.Average time per iteration spent for voronoi diagram construction for (a) regular particle distribution and (b) Quadtree particle distribution.

Fig. 15 .
Fig. 15.Volume minimization with stress constraint for the L-bracket example: Evolution of the structure with von Mises stress field.

Fig. 18 .
Fig. 18.Snapshots of the piston-head solution with pressure loads.

Fig. 19 .
Fig. 19.Voronoi diagram for the piston-head example: (a) Initial diagram and (b) Diagram at the optimum solution.