Proportional topology optimisation with maximum entropy-based meshless method for minimum compliance and stress constrained problems

In this paper, proportional topology optimisation (PTO) with maximum entropy (maxent)-based meshless method is presented for two-dimensional linear elastic structures for both minimum compliance (PTOc) and stress constraint (PTOs) problems. The computation of maxent basis functions is efficient as compared to the standard moving least square (MLS) and possesses a weak Kronecker delta property leading to straightforward imposition of Dirichlet boundary conditions. The PTO is a simple, non-gradient, accurate, and efficient method compared to the standard topology optimisation methods. A detailed and efficient implementation of the computational algorithms for both PTOc and PTOs is presented. The maxent basis functions are calculated only once at the start of simulation and used in each optimisation iteration. Young’s modulus for each background cells is calculated using the modified solid isotropic material with penalisation (SIMP) method. A parametric study is also conducted on the degree of proportionality and history dependence of both PTOc and PTOs algorithms. A variety of numerical examples with simple and complex geometries, and structured and unstructured discretisations are presented to show the accuracy, efficiency, and robustness of the developed computational algorithms. Both PTOc and PTOs algorithms can handle large topological changes, and provide excellent optimisation convergence characteristics.


Introduction
Topology optimisation (TO) is an effective tool for the optimal solutions of structural engineering problems with a variety of objective functions, e.g., minimum compliance (maximum stiffness), minimum stress, minimum weight, minimum thermal potential energy, etc. The required objective function is achieved through an iterative process which distributes the given amount of material within the design domain imposed through a volume/area constraint [1]. The increasingly high demand of lightweight and cost-efficient structural designs in many industrial applications, e.g., automotive, aerospace, wind turbines, buildings, etc., can rightly be fulfilled with TO [2,3]. After the pioneering work of Bendsoe and Kikuchi [4], a variety of TO methods have been developed, such as solid isotropic material with penalisation (SIMP), evolutionary structural optimisation (ESO) method, bi-directional evolutionary structural optimisation (BESO) method, and level set-based optimisation method [4][5][6][7].
Majority of the existing TO techniques used finite-element method (FEM) for the solution of underlying structural mechanics or thermal problems; for details, see [8,9]. However, in the literature, meshless methods are also used as alternative to FEM [10][11][12][13]. As compared to FEM, only a set of nodes is required in meshless methods. In addition, meshless methods provide flexible control of approximation and ease in performing local adaptive refinements. Meshless methods have been used for the solution of a large number of engineering problems such as large deformation, fracture mechanics, and contact mechanics problems [14][15][16]. The reproducing kernel particle method (RKPM) coupled with variable density techniques was proposed for the TO of geometrically non-linear structures in [12]. The meshless corrected smoothed particle method (CSPM) was proposed for the TO of plane structure in [17]. Optimisation is based on a modified SIMP approach with a sensitivity filter to avoid the checker-board pattern. The meshless smoothed particle hydrodynamics (SPH) method coupled with SIMP method with a filtering technique was used for the optimisation of linear structures under multiple load cases in [18].
The computational cost of the element-free Galerkin method (EFGM)-based TO algorithm was improved using coupled EFGM and FEM in [19]. The EFGM in combination with SIMP was used for the TO of thermo-mechanical compliant mechanisms with geometrical non-linearities in [13]. Compliance minimization under stress constraint with multiple load condition was studied in [11]. A meshless Galerkin level set method (LSM) was proposed for structural shape and TO problems in [20], where compactly supported radial basis functions (CSRBFs) were used to construct the basis functions for meshless approximations and for parametrisation of the level set function. The EFGM with moving least square (MLS) basis functions with dual-level point-wise density was implemented in [20]. In [21], the EFGM was coupled with LSM for two-dimensional TO problems with single and multiple load cases. For hole insertion inside the structure during the optimisation process, the topological derivative term was added in Hamilton-Jacobi equation. TO of continuum structures using BESO-based EFGM was implemented in [22]. The EFGM for the TO of heat conduction problems was proposed in [23,24].
Stress-based TO is considered as one of the important areas of structural design. If stress is not taken into account during the conceptual design stage of the structural component, the design will require post-processing or rework and will lead to unexpected costs [25]. For the stress-based structural shape optimisation and TO was proposed in [26]. A mixed finite-element method for stress constraint TO was proposed in [27]. Both displacements and stresses were used as variables in the formulation. The ESO method for stiffness maximisation and stress minimisation was developed in [28]. The proposed method was applied to both single and multiple load cases.
Most of the optimisation methods presented in the literature are based on the evaluation of gradients of objective function and the imposed constraints. In the case of minimum compliance-based TO, design is updated using information from the associated sensitivities obtained through the gradient information of objective function. In the compliance-based TO, a self-adjoint option makes the gradients evaluation as the evaluation of the displacement field only. However, in the case of a minimum stress-based objective function approach, the evaluation of the gradients is based on complex analytical expressions, which are computationally expensive as well as may pose implementation issues [29]. The evaluation of gradient information enhances the convergence speed of an optimisation approach; however, there exist non-gradient techniques with lower implementational complexities and comparable efficiencies [30]. A proportional topology optimisation (PTO) method is presented in [29], which is a non-gradient-based approach for the solution of minimum compliance and stress constraint-based TO problems. The heuristic nature of the proposed approach allows simplicity in its implementation and understanding in addition to obtaining efficient and accurate optimal solution. An improved PTO method for solving the minimum compliance problems was presented in [31], with some advantages in terms of certain performance aspects as compared to the algorithms proposed in [29,32].
In this paper, PTO with maximum entropy (maxent)-based meshless method is presented for two-dimensional linear elastic structures for both minimum compliance (PTOc) and stress constraint (PTOs) problems. The maxent basis functions possess a weak Kronecker delta property leading to straightforward imposition of Dirichlet boundary conditions as in the case of finite-element analysis [33][34][35][36][37][38]. Furthermore, the calculation of maxent basis functions and corresponding spacial derivatives are more efficient as compared to the standard moving least square (MLS) counterpart [35]. The PTO is a simple, non-gradient, accurate, and efficient method as compared to the standard TO methods [29]. In the PTO method, constraint is imposed on the global structure. Due to its global nature, PTO is considered as heuristic method and search for optimal solution globally. This heuristic nature of PTO is leading to its straightforward implementation in computer program as compared to alternative TO methods. In addition, the PTO does not require the calculation of sensitivities. A detailed and efficient implementation of the computational algorithms for both PTOc and PTOs is presented. The maxent basis functions and its spacial derivative are calculated only once at the start of simulation and used in each optimisation iteration for the calculations of stiffness matrix, stresses, and compliance. Young's modulus for each background cells is calculated using the SIMP method and is then subsequently used for the calculation of stiffness matrix for structural analysis. A parametric study is also conducted on degree of proportionality and history dependence of the both PTOc and PTOs algorithms. Open-source pre-processor, CUBIT, is used with the developed optimisation algorithms, which allow to construct complex geometries and use robust meshing algorithms [39]. A variety of numerical examples with simple and complex geometries with structured and unstructured discretisation are presented to show the accuracy, efficiency, and robustness of the developed computational algorithms. Both algorithms can 1 3 handle large topological changes and provide excellent optimisation convergence characteristics. This paper is organised as follows. A computational framework is given in Sect. 2. The computational framework includes an overview of the maxent-based meshless method which is described in Sect. 2

Computational framework
The developed computational framework consists of maxentbased meshless method for elasto-statics problems and associated PTO method. Mathematical formulations and computational details of these are provided in the following.

Maximum entropy meshless method
An overview of the maxent-based meshless method is given in this section, full details of which is available in [33][34][35][36][37][38]. A two-dimensional formulation is provided, which is straightforward to modify for one-and three-dimensional problems. The final discretised system of equations is written as where For the plane strain and plane stress cases, matrix is written as where is the Poisson's ratio and E is the Young's modulus, and f and are body force and traction, respectively. To perform the integrations in Eqs. (2) and (3) numerically, the problem domain Ω and traction boundary Γ t are divided into a number of non-overlapping background cells. In Eq.
(1), the vector contains fictitious nodal values or nodal parameters, and and are global stiffness matrix and force vector, respectively. The maxent basis functions possess a weak Kronecker delta property leading to straightforward imposition of Dirichlet boundary conditions. The original maxent shape functions derived for the polygonal interpolation in [33] are not ideal for meshless methods as they extend to the global problem domain. The non-local and non-interpolating characteristics of these shape functions are highlighted in [34] where the (more useful) local maxent formulations are introduced and incorporated into meshless modelling of linear and non-linear elasticity. The weak Kronecker delta property of these shape functions was also observed, and its correspondence with the MLS was also highlighted. Compact support shape functions are derived using Gaussian weight functions (or priors) in [34], work which is extended in [40] for any weight function (or generalized prior). In this paper, the local maxent formulation based on [40] is used. A detailed derivation of both global and local maxent basis function is also given in [35]. The i is a matrix of the local maxent basis functions for node i at a point , that is The maxent shape functions is defined as where in which w i is the weight function associated with node i, evaluated at point = (x, y) T , x i = x i − x and ỹ i = y i − y are shifted coordinates. n is the number of nodes in support of Gauss point (nodes, the influence domains of whose encompass the Gauss point) and λ 1 and λ 2 are Lagrange multipliers which can be found from (6b) = for plane stress, 1− for plane strain,

3
F is a convex function and Newton's method is used to solve (10) to find the Lagrange multipliers which can then be used in the expressions for the shape functions [33,34,40]. In this paper, cubic spline weight function is used. The cubic spline weight function in one dimension is given as follows: the distance between the node i and point of interest x. The d mi is the domain of influence of a node and is a very important concept in meshless methods, as it determines the region in which it has influence. The size of domain of influence for For static analysis, its value ranges from 2.0 to 4.0 and c i is determined by searching for enough neighbour nodes.
The shape function derivatives follow as: where where is the Hessian matrix and the dyadic product ⊗ of two vector and is a second order tensor, i.e., ⊗ defined as ⊗ = T .

Proportional topology optimisation
Mathematical background and computational algorithms for both PTOs and PTOc and their coupling with maxentbased meshless method are given in the following sections. Both PTOc and PTOs algorithms are implemented in MATLAB for two-dimensional problems.

Stress constraint topology optimisation
The PTOs search for the optimal design with minimisation of available material or volume fraction for constraining the stress to user-defined limit. Mathematically, PTOs is written as [29]: Here, n BC is the number of background cells used for the discretisation of design domain, and is the design variable or density. i and v i are density and volume/area of background cell i. i = 1 means solid material and i = 0 refer to as hole. max = 1 and min = 0.001 are the maximum and minimum values of the design variable. min = 0.001 is used instead of min = 0 to avoid singularity of the stiffness matrix. The , , and are standard global stiffness matrix, nodal parameters vector, and external force vector, respectively. i is a stress measure calculated at the centre of each background cell, while l is the user defined stress limiting value. von Misses stress measure is used in this paper.
The computational algorithm for the PTOs is shown in Fig. 1. The algorithm starts with setting up the problem. This step consists of either creating discretisation information for simple geometries directly within the MATLAB program or importing these from the input files created using CUBIT for complex geometries. The discretisation information consists of data for background cells, nodal position, and Dirichlet and traction boundary conditions. Domains of influence for each node within the design domain are subsequently calculated. The maxent basis functions and their spacial derivatives are calculated and saved to be used in each optimisation iteration. External force vector is also calculated only once at the start of simulation. In each optimisation iteration, Young's modulus (E) for each background cells is calculated from its density ( ) using the SIMP method [41]. The calculated E for each background cell is then used for the calculation of its stiffness matrix. The modified SIMP approach is given as Here, E min is minimum value of Young's modulus, typically 10 −9 , and is assigned to void background cells to avoid singularity of global stiffness matrix . After the solution of final system of equations, i.e., Eq. (1), stresses and compliance are calculated. In the case of two-dimensional problem, stresses = [ x y xy ] T are calculated at the centre ( ) of each background cell from the nodal parameters i = [u xi u yi ] T as follows: where x and y are normal stresses in x and y direction, respectively, and xy are shear stress. Equations for matrices and i are given in Sect. 2.1 and n are the number of nodes in the support of point . The von Misses stress is subsequently calculated at the centre of each background cell using The maximum vM is compared with user-defined allowable stress limit l . If the convergence criterion satisfies, i.e., vM − l = , then the analysis is terminated, where is a small number. In this paper, = 0.01 is used. If the convergence criteria are not satisfied, the simulation continues to add or remove material from the design domain. If vM > l material is added to the design domain; otherwise, material is removed, that is where TM is the new target material for the design domain, CM is the current material, and MM is the move material. In this paper, MM = (0.001 × number of background cells ) is used. The calculated TM is then iteratively distributed to the design domain within the second internal loop in each optimisation iteration. Material is proportionally assigned to each background cell according to its von Misses stress where RM = TM − AM is the remaining material, AM is the actual material, and q is the degree of proportionality. Filter is applied subsequently to the density field. A cone-type filter is used in this paper [29]. For each background cell, density ( i ) is calculated using the following equation: Due to the use of complex geometries and unstructured discretisation, Eq. (21) is implemented using influence domain associated with each background cell. Domain of influence d max = r 0 c i is calculated for each background cell, where c i is found by searching for enough neighbouring background cells and r 0 is a user-defined scaling parameter. r 0 = 1.2 is used to include enough neighbouring background cell in the calculation. At the same time, this smaller value of r 0 will ensure local characteristics. c i = max(r 1 , r2) where r 1 and r 2 are two diagonal of background cell. w ij is the weigh function shown in Eq. (11), n is the number of neighbouring background cells in the support of cell i.
, where i and j are spacial coordinates of geometric centres of background cells i and j, respectively. This density filtering is local averaging and preserve volume of the design domain.
The RM distribute into background cells without regard to the density limits of 0 and 1. Subsequently, density limits are applied and actual material is calculated using The calculated AM is then used to calculate the remaining material using RM = TM − AM . Due to this difference of TM and AM, iterative procedure is adopted for the distribution of material to the design domain. The convergence criteria for the second internal loop are RM < , where = 0.0001 is used here. Finally, the calculated where new i is the new density of the background cell i, and is a coefficient for blending densities between two optimisation iterations. = 0 represents no contribution from the previous optimisation iteration and = 0.5 represents half of the contribution from the previous optimisation iteration.

Minimum compliance topology optimisation
In the case of PTOc, optimisation algorithm search for optimal design with minimum compliance subjected to volume constraint such that where C is the compliance of the whole structure and can be calculated from the summation of compliance of each background cell ( ∑ n BC i=1 T i i ) or directly for the global problem T . The computational algorithm for the PTOc is shown in Fig. 2.
In the PTOc, the target material TM is fixed during the whole optimisation process and determined a priori using a user-defined volume fraction, that is where v f is the user-defined volume fraction. TM is distributed proportional to the compliance in each background cell, that is where RM is the remaining material, opt i is the optimised density, and C i is the compliance of background cell i and q is the degree of proportionality. At each optimisation iteration, compliance is calculated for each background cell, and then, the TM is iteratively distributed to each background cells proportional to their compliance. The convergence criteria in the case of PTOc are given as (22)   the MATLAB program. In the final numerical example, complex geometries and unstructured discretisations are considered including serpentine beam, hook problem, and wrench problem. CUBIT is used to create these geometries and input data for their discretisations.

MBB beam
The first numerical example is MBB beam. Geometry, boundary conditions, and loading for this problem are shown in Fig. 3. For this problem, the input parameters used are  Table 1 can provide a reasonable optimised geometry, as shown in Fig. 4a. The optimised geometries for = 0.5 and q = (0.5, ⋯ 3.0) are shown in Fig. 4a. For q = (0.75, 1.0, 1.25, 1.5) , it gives expected optimised geometries, but for q > 1.5 , PTOc generates unreasonable optimised geometries with very high compliance values as given in Table 1. For = (0.4, 0.6, 0.7, 0.8, 0.9) and q = 3.0 , optimised geometries are also given in Fig. 4a. For = (0.7, 0.8, 0.9) and q=3.0 it provides reasonable optimised geometries, but for = (0.4, 0.6) and q=3.0 it leads to unreasonable optimised geometries with very high values of compliance.
Similarly, obtained compliance values and volume fractions for the PTOs algorithm are given in Tables 2 and 3, respectively. The minimum value of compliance in Table 2 is 249. 38 and compliance values within 10% range are shown as bold text. Similar to the PTOc, for any combination of and q shown as bold text in Table 2, PTOs provide reasonable optimised geometries. In the case of PTOs, volume fractions of the optimised geometries for each versus q are shown in Fig. 4c. Lower volume fractions are obtained with higher values of , but this leads to increase in compliance as given in Table 2. The optimised geometries for = 0 and q = (0.5, ⋯ 3.0) are shown in Fig. 4b. Similarly, optimised geometries for = (0.1, 0.2, 0.3) and q = 3.0 are also shown in Fig. 4b. With few exceptions, all these are reasonable optimised geometries.
For the MBB problem, discretisation are uniformly refined to assess its effect on the optimised geometries and     [29]. Although, the overall optimised geometries are very similar for all discretisations but with increasing refinement, extra details appear for discretisation with n BC of (50 × 150) as shown in Fig. 5e. For all these five discretisations, compliance and max( vM ) versus number of iterations are shown in Fig. 6a and b, respectively. For each discretisation, the final compliance value, max( vM ) and number of iteration required for convergence are also given in Table 4. As expected, compliance is decreasing and max( vM ) is increasing with increasing refinement level. For the PTOs, in addition to the previously mentioned parameters, = 0 and q = 2 are used. For all five discretisations, optimised geometries are shown in Fig. 7. For each discretisation, l is used from the corresponding converged value of max( vM ) of the PTOc algorithm. All discretisations provide similar results to the one reported in the literature [29] and the one obtained with the PTOc algorithm shown in Fig. 5. The only exception is Fig. 7a, which is due to the use of very coarse discretisation with n BC of (10 × 30) . In this case, compliance and max( vM ) versus iterations is shown in Fig. 8. In addition, final compliance values and the number iterations required for convergence for all the discretisations are given in Table 4. As compared to the PTOc, the compliance value for the PTOs is lower for discretisations with n BC of (10 × 30) and (20 × 60) , but for the remaining discretisation with n BC of (30 × 90) , (40 × 120) , and (50 × 150) , the compliance values are higher for the PTOs. In the case of PTOs, except for the discretisation with n BC of (10 × 30) , less number of iterations are required for convergence.
Local maxent basis functions are used in this paper. Figure 9a and b shows compliance and max( vM ) versus number of iterations, respectively, for the PTOc using meshless with different d max and finite-element analysis. Only 100 iterations are considered. In addition, the final values of compliance, max( vM ) and simulation time are also given in Table 5. The compliance versus iteration plot in Fig. 9a is very similar for meshless with different d max and finite element. The max( vM ) versus iteration plot is very similar for meshless with d max = 1.5 and the FEA, as shown in Fig. 9b.
With increasing values of d max from 1.5 to 3.0, max( vM ) versus number of iterations plot converges. For d max = 3.5 , max( vM ) plot starts to diverge and the stress values are even higher than for meshless with d max = 1.5 and the FEA. From Table 5, it can be seen that simulation time for meshless with d max = 1.5 and FEA is very similar. In the case of the FEA, shape functions and their derivatives are calculated   Similarly, compliance and max( vM ) versus iterations for the PTOs using meshless with different d max and the FEA are shown in Fig. 10. The final values of compliance, max( vM ) , and simulation time are also given in Table 6. Similar to Fig. 9, compliance values are very similar for meshless with different values of d max and the FEA. In the case of max( vM ) , stress values converge for meshless with increasing d max up to 3.0 and start to diverge for d max = 3.5.

Cantilever beam
The second numerical example is cantilever beam subjected to point load in the middle of the free edge. Geometry, boundary conditions, and loading for this problem are shown in Fig. 11. Input parameters used for this problem are the same as used in the MBB problem for both PTOc and PTOs algorithm. For the PTOc algorithms, = 0.5 and q = 1 are used. Four discretisation with n BC of (30 × 60) , (40 × 80) ,  versus iterations are shown in Fig. 12a and b respectively. The final optimised geometries are also shown in Fig. 12b. The final compliance, max( vM ) , and number of iterations   are also given in Table 7. The final optimised geometries are in a very good agreement with the one achieved in the literature for the same problem [29]. The problem converges within reasonable number of iterations. As expected, with increasing refinement level, compliance values decreases, while max( vM ) increases. For the PTOs, the same four discretisations are also considered with = 0 and q = 2 . For each discretisation, l for the PTOs is used from the corresponding max( vM ) of the PTOc algorithm. Compliance and max( vM ) versus iterations are shown in Fig. 13a and b, respectively. The final optimised geometries are also shown in Fig. 13b. In addition, the final compliance values and number of iterations required for convergence are also given in Table 7. For = 0 and q = 2 , the PTOs algorithm struggles to converge and required large number of iterations to achieve the same stress level as in the case of PTOc.The final optimised geometries are not in a very good agreement with the one given in literature for the same problem [29] and corresponding results obtained from PTOc shown in Fig. 12b. This problem is subsequently solved with PTOs with = 0 and q = 1 , for which compliance and max( vM ) versus iterations are shown in Fig. 14a and b, respectively. The final optimised geometries are also shown in Fig. 14b. In addition, the final compliance values  and number of iterations required for convergence are also given in Table 7. All these optimised geometries are in a very good agreement with the ones given in the literature [29] and with those obtained using the PTOc shown in Fig. 12b. The convergence in this case is very quick and it only takes very few iterations as compared to the PTOc and PTOs with = 0 and q = 2 . As compared to the compliance and max( vM ) versus iteration plots shown in Fig. 13, plots in Fig. 14

L-shaped bracket
The third numerical example is L-shaped bracket for which geometry, boundary conditions, and loading are shown in Fig. 15. It is fully fixed at the top edge and a point load is applied at the right corner. Input parameters used for this problem are also the same as used in the MBB problem for both PTOc and PTOs algorithms. In this case, four discretisations with n BC of 567, 1600, 6400, 14400 and number of nodes of 637, 1701, 6601, and 14701 are considered. For = 0.5 and q = 1 , compliance and max( vM ) versus the number of iterations are shown in Fig. 16a and  b, respectively. The optimised geometries are also shown in Fig. 16b. The optimised geometries are in a very good agreement with the one given in the literature [29]. The problem converges within reasonable number of iterations. Similar to the previous problems, with increasing refinement, compliance values decreases, while max( vM ) increases. Additional details can also be observed in the links on the bottom side for the fine discretisation with n BC of 14,400 and number of nodes of 14,701.
In the case of PTOs, for = 0, q = 2 and = 0, q = 1 , compliance versus iterations are given in Figs. 17a and 18a, respectively. For the same parameters, max( vM ) versus iterations and final optimised geometries are shown in Figs. 17b and 18b, respectively. The PTOs with = 0, q = 1 converge very quickly as compared to PTOc and PTOs with = 0, q = 2 and leading to optimised geometries which are very similar to the one given in literature and ones obtained using the PTOc.

Problems with complex geometries
The first three numerical examples are with relatively simple geometries and structured discretisations. In this numerical examples, both PTOc and PTOs algorithms are tested for complex geometries with unstructured discretisations to show the robustness of the developed computational framework. These geometries includes serpentine beam, hook problem, and wrench problem [42][43][44][45]. Geometries, boundary conditions, and loadings for serpentine beam, hook problem, and wrench problem are shown in Fig. 19a-c, respectively. Background cells and nodes for these problems are shown in Fig. 20. The discretisation of serpentine beam, hook problem, and wrench problem are (n BC = 5660, nodes = 5842) , (n BC = 5802, nodes = 6455) , and (n BC = 4448, nodes = 4655) respectively. In the case of PTOc, = 0.5 and q = 1 are used, while in the case of PTOs, = 0 and q = (0.5, 1, 2.0) are considered. For the hook problem, plots for compliance and max( vM ) versus iterations are shown in Fig. 21a and b, respectively. The final optimised geometries in the case of PTOc, PTOs ( q = 2 ), PTOs ( q = 1 ) and PTOs ( q = 0.5 ) are shown in Figs. 22a-d, respectively. These optimised geometries are in a good agreement with the one obtained in literature for similar problem [42][43][44][45]. Compliance and max( vM ) plots for PTOc are very smooth and converge in    PTOc, but material is inappropriately distributed and there are many background cells with density between 0 and 1.
Similarly, serpentine beam and wrench problem is solved with the same parameters. For the serpentine beam, plots for compliance and max( vM ) versus iterations are shown in Fig. 23a and b, respectively. The optimised geometries in the case of PTOc, PTOs ( q = 2 ), PTOs ( q = 1 ), and PTOs ( q = 0.5 ) are shown in Fig. 24a- Fig. 25a and b, respectively. The final optimised geometries in the case of PTOc, PTOs ( q = 2 ), PTOs ( q = 1 ), and PTOs ( q = 0.5 ) are shown in Fig. 26a-d, respectively. These optimised geometries are in a very good agreement with the one obtained in the literature [42][43][44][45]. Similar to the hook problem plots for PTOc and PTOs with q = 1, 0.5 are very smooth but results for PTOs with q = 2 fluctuates.

Conclusion
Maximum entropy-based proportional topology optimisation (PTO) is presented in this paper for two-dimensional linear elastic structures for both minimum compliance (PTOc) and stress constraint (PTOs) problems. The use of maximum entropy basis functions provides efficient calculation of basis functions and its spacial derivatives and leading to straightforward imposition of Dirichlet boundary conditions as in the case of finite-element analysis. The non-gradient nature of the PTO method makes it convenient to implement in a computer program and at the same time provides very accurate results efficiently as compared to the standardised TO methods. The calculation of maxent basis functions and its spacial derivative only once at the start of simulation improve the computational efficiency. In addition, background cell-based approach, i.e., the calculation of cell-based density, compliance, and stress, add to the computational efficiency. A variety of numerical examples with both simple and complex geometries, and structured and unstructured discretisations are presented to show the accuracy, efficiency, and robustness of both PTOc and PTOs algorithms. It was found that both  algorithms can handle large topological changes and provide excellent optimisation convergence characteristics in terms of compliance and max( vM ) versus iterations curves. The optimised geometries are found to be in a very good agreement with the one found in the literatures. It was shown that PTOs with appropriate proportionality and history dependence parameters provide excellent results compared to PTOc in terms of compliance values and number of iteration required for convergence. It was also shown that for PTOc algorithms compliance decreases and stress increases with increasing refinement level. In addition, it was shown that final optimised geometries are very similar for different refinement level. Due to the use of coupled meshless method and PTO, the developed computational framework provides a solid foundation to be extended for the solution of three-dimensional and nonlinear problems (both material and geometrical) and will be pursued in the future. Moreover, the current computational framework for isotropic material behaviour will  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Fig. 26
Optimised geometries for the wrench problem