A new algorithm for dense ellipse packing and polygonal structures generation in context of FEM or DEM

. A new constructive ellipse packing algorithm is presented. It allows to respect the imposed area, shape and spatial orientation distribution (i


" #"
The problem of packing of objects with different geometrical shape attracted considerable amount of attention during the last few decades due to its application in modelling of granular materials, porous media, foams, polycrystalline structures and so on.Most studies were carried out for the packing of uniform and arbitrary sized spheres [1], considering the problem of optimal packing [2][3][4][5], (and its dual, the sphere cutting problem [6,7]) or space-filling [8][9][10][11][12][13][14].
In the context of space-filling problem, the sphere packing can be obtained by the dynamic methods, in which particles usually change their position and/or their size during the filling process to improve the local and global densities, or by the constructive methods, generating particles sequentially and calculating their positions thanks to geometrical rules.
The dynamic methods were widely used for the ellipsoid packing and high densities were reported in published works [15][16][17][18].For the methods in which ellipsoidal particles are generated by choosing randomly initial position and orientation and letting them fall in an open box under the gravity, an anisotropic behaviour of the packing system was reported [18,19].In the case when the spatial orientation of ellipses is fixed, it was found that the packing by dynamic methods can result in low densities [20].
The constructive method of axisymmetric elliptical particles approximated by four connected arcs was proposed in [21] and it was extended to three dimensional (3D) axisymmetric ellipsoidal forms in [22].In [23] the authors developed an advancing front packing algorithm.
Recently, an overlapping detection algorithm for the sequential packing of elliptical particles was proposed and generalized in [24].Relatively high packing densities were reported, but the generation of packing with an imposed spatial orientation distribution was not considered.It should be noted that due to numerous rejections of randomly generated positions of particles during the filling process in the sequential inhibition method, the computation time for achieving high and even moderate packing densities is typically long.
A constructive ellipse packing method is proposed in the present work.Each elliptical particle is approximated by a set of circles.Recently, this method was examined for applying in discrete element simulations [25].In the present study, we focus on the packing of particles which can be applied not only to consider modelling of granular materials, but also for the generation of polygonal structures (which can be used, for example, to describe polycrystalline materials statistically).In this case, higher packing densities are necessary to provide good accuracy in respecting the desired polygon area, shape and spatial orientation distributions.In the proposed packing method, the list of accommodated in domain ellipses is built in which they are sorted according to the position in the domain (height).To reach high packing densities, each new particle is placed at the lowest possible position by testing pairs of ellipses from the "Height list".
The paper is organized as follows.At first, the discretization of ellipses on multi-circles is described.Then, the algorithm of packing is presented with introducing the Height list used to optimize packing density.Some numerical tests are reported for imposed area, shape and spatial orientation distributions.Finally, densities and computation time provided by the proposed algorithm are compared with the Optimized Dropping and Rolling Method (ODRM) [9] in particular case of disk packing.

! "
In this section, the approach used to approximate an ellipse by a set of circles is described in the first part.In the second part, the algorithm of ellipse packing is detailed.

#" !
Assuming that the principal axes coincide with the axes of local Cartesian coordinates Oxy, the corresponding ellipse equation is: where a and b are the major and minor half-lengths of principal axes.
A random ellipse can be modelled by a set of circles as follows [25]: -All the centers of circles are on the major axis and, hence, the circle equation is , where x i is the abscissa of the center and R i is the radius of the circle with i = 0 for the central circle and i > 0 for the right part of the ellipse.
-All the circles are tangent to the ellipse.If (x it ,y it ) is the point of tangency for the circle i with center (x i ,0), after calculating the derivative with respect to x for the left and right parts of Eq.(1) and Eq.(2) and substituting x = x it , we obtain: , which gives: , with λ = a / b the ellipse shape ratio.Consequently, the radius of the circle R i can be found from Eq.(1) and Eq.( 2) by eliminating y 2 and substituting x = x it where x it is given by Eq.( 5): -The center of the largest circle (0,0) is in the center of the ellipse and its radius is equal to the semi-minor axis b.
-The center of the circle which is the nearest neighbour of the largest (central) circle, (x 1 ,0), is imposed by the user (the closer it is to the central circle, the better the final approximation of the ellipse is obtained).
If Eq.( 7) leads to x i+1 > x N , we impose x i+1 = x N and the iteration is terminated.After reflecting the circles {C i ((x i ,0),R i ), i=1..N} over the y-axis, the description of the ellipse by a set of 2N+1 circles is obtained.The total number of circles approximating an ellipse depends on the ellipse shape ratio, λ, and the precision imposed by the user, x 1 .

"
The proposed ellipse packing method, called Advancing Layer algorithm (ALA), is described in this subsection.The flow-chart of the algorithm for placing the first row (ellipses are in contact with the bottom of the domain) is summarized in Fig. 2. When the first row is completed, the bulk of the domain is filled as described in the flowchart detailed in Fig. 3.
Each new ellipse is generated according to the imposed area, A, shape (aspect ratio), λ, distribution.Then, the ellipse is decomposed on circles as it is described in the previous subsection with the precision imposed by the user and rotated by the angle α (see Fig. 4) generated according to the imposed spatial orientation distribution.
To place the ellipse E i in contact with the bottom edge, the lowest circle should be determined and ycoordinate of its centre is imposed to be equal to its radius.Then, y-coordinates of other circles of the particle are computed.When the particle is placed in the left     bottom corner of the domain, x-coordinates are defined analogically.
After accommodating the first elliptical particle in the corner, other ellipses are generated and placed as described in the flow-chart for filling the first row (Fig. 2).Fig. 5 illustrates an example of two particles placed in contact with bottom: after placing E 1 in the corner, the position of E 2 is calculated by searching one of its circles to be in contact with the circles of E 1 with no overlapping of other circles.
When the first row of ellipses is filled, the index number of each particle is recorded in the Height list at the position which corresponds to the height of the particle (its maximal y-coordinate).The ellipses in the Height list are sorted in ascending order, i.e. from the lowest ellipse (the first position in the list) to the highest one (the last element in the list).
For completing the bulk of the domain, the flow-chart described in Fig. 3 is used for any tested ellipse.After its generation (following imposed distributions), its decomposition and rotation, a new ellipse should be placed in contact with two other particles.Pairs of previously placed ellipses are tested starting from the first element to the last one recorded in the Height list.Fig. 6 illustrates the placement of the ellipse E 3 on the pair of E 1 (with the height H E1 ) and E 2 (with the height H E2 ).The contact position is determined from two conditions of the contact for two pairs of circles (C q E1 with C g E3 , C p E2 with C f E3 in Fig. 6) without any intersection with other circles.Two equalities of the sum of radii to the distance between centres of the corresponding circles from the chosen pairs give two equations with two unknown variables.The system can have two, one or zero solutions.Two solutions correspond to two possible contact positions for the new ellipse -"beneath" and "above".When the new ellipse is placed in contact with the two others, the intersection with other previously placed particles should be checked.In the case of no intersection, the position is accepted and the index number of the new ellipse is recorded in the Height list according to its computed height.In Fig. 6 the ellipse E 3 is placed above E 1 and E 2 (since the position beneath E 1 and E 2 cannot be accepted due to the intersection with other particles).When the overlapping with other particles is detected, another pair of elements is chosen from the Height list.If after testing all possible pairs of ellipses for placing the new particle the contact position was not accepted, the packing is stopped.
In the presented form, the packing algorithm is quite time consuming and generates packings with inhomogeneous densities due to the fact that the pairs of particles to place the new ellipse are tested from the bottom to the top.This strategy is synonymous of higher packing density in the bottom part and in longer computation time since a lot of unsuccessful tests (with overlapping) should be performed before to reach the good position for a given new ellipse.To avoid it, the definition of the Top layer was introduced.It is defined as follows.Let H max be the height (maximal y-coordinate) of the highest ellipse placed in the domain (the last element in the Height list) and H Ei be the height of the ellipse E i (see Fig. 6).The Top layer of the packing is

composed of the particles E i for which
where a max is the maximal semi-major axis of the ellipses and h TL is a coefficient imposed by the user.The Top layer contains all the ellipses recorded in the Height list at last positions.The number of particles in the Top layer depends on the value of h TL .In the present work, the values from 1.1 to 2 were tested to optimize the computational time and the global packing density.In Fig. 7 the new ellipse E 3 is placed on two previously accommodated ellipses: the first ellipse E 1 is chosen from the Top layer and the second ellipse is chosen from the rectangular domain (the biggest rectangle in Fig. 7. called further candidate domain) defined by the rectangular cells of E 1 and E 3 (two small rectangles in Fig. 7).Fig. 7 shows the case when the ellipse E 4 is not tested for placing E 3 since it does not intersect the candidate domain, the ellipse E 2 overlaps with the candidate domain and, hence, can be considered for placing E 3 .
In addition, when the tests of intersection with the previously placed ellipses are performed, first the overlapping of the rectangular cells of ellipses is checked.If the cells overlap, the intersection test for the ellipses is then performed.
In the described form, the ALA was applied for dimensionless ellipse and disk packings and the results are reported in the next section.

!#"! !
Most of the works dealing with ellipse packing consider the case when ellipse orientation (ellipsoids in 3D case) is not conserved during the filling process.According to NUMIFORM 2016 [20], when particles are allowed to rotate, the ellipse packing density increases when the ellipticity increases until reaching a maximum value and then drops below the density of disk packing.Higher densities compared to the packing of disks can then be obtained but without respecting data concerning α.When no rotation is allowed, it was also shown that the packing density decreases as the ellipticity increases.
In this 2D work, we focus on the problem of ellipse packing with spatial orientations conserved during the filling.A uniform ellipse area distribution law with normal and uniform α distributions on ]-π/2, π/2] range for the angle of orientation of the semi-major axis packing of ellipses have been used to generate packing by using the ALA.Depending on the aspect ratio, 4600 to 5100 ellipses were generated (a high number of particles is used for better respecting the imposed distributions).Evolution of the packing density as a function of aspect ratio is described in Fig. 8.These results are in agreement with [20]: density decreases with increasing aspect ratio.The obtained values are lower in the case of the uniform α distribution than the ones with normal α distribution.One may conclude that for a material with a preferential direction of ellipses elongation, packing with higher densities can be generated.
Two examples of ellipse packing with imposed uniform A and α distributions are shown in Fig. 9 and Fig. 10 for the aspect ratios λ = 1.5 and λ = 4, respectively.One can see that there is no preferential direction in the orientation of the ellipses.More voids are observed in the case of higher ratio which is clearly indicated by the obtained packing densities: 0.7284 in the case of λ = 4 and 0.8133 for λ = 1.5.
Finally, to demonstrate the performance of the ALA, it was tested for a more complex case when a bimodal distribution is imposed for all three parameters, A, λ and α.The size of the considered domain is 30×30.A varies from 0.0001 up to 0.181 with two modes of 0.0386 and 0.1271.λ is in the range of 1.5 to 4.4375 with two modes of 2.125 and 3.5625.Finally, α varies from -π/2 to π/2 with two modes of -0.287π and 0.2π.The generated packing of 10696 ellipses with a density of 0.7979 is shown in Fig. 11.

! $" "' "
To demonstrate the efficiency of the ALA for disk packing it was compared with the ODRM [9].The packings of disks by ODR and ALA were generated with an imposed log-normal radius distribution.The distribution width was fixed to R max /R min = 48, where R max and R min are, respectively, the maximal and the minimal disk radius.To define the top layer in the ALA, the parameter h TL was chosen to be equal to 1.1.The computation results are presented in Fig. 12 and Fig. 13.One can see that the ALA allows to place disks in the domain more homogeneously with higher densities compared to ODR.Packing densities obtained by ODR and ALA are 0.8083 and 0.8493, respectively.The CPU time needed for the ALA computation (1.9 s for the packing of 5220 disks) is about two times longer than the   time required for disk packing by ODR method (about 1 s for the packing of 5139 disks).The computation time for the ODR (t ODR ) and ALA (t ALA ) was compared for filling domains of different sizes (from 2000 to 110000 disks) with the same (log-normal type) size distribution.In all the cases, the time required for the AL method is about two times longer.However, the analysis shows that the ratio t ALA /t ODR decreases when the size polydispersity R max /R min decreases.This dependence is illustrated in Fig. 14 (red curve).The ALA is about two times faster than ODR algorithm when the size polydispersity is equal to 4 and about two times slower when R max /R min = 50.This property is explained by the fact that, when the variability in disk size decreases, the number of particles tested in the top layer for placing a new one decreases  also (since less small particles are generated).The same analysis was performed for other types of distributions.
The results are shown in Fig. 14 as well.In contrast to the case of log-normal type A distribution, the ALA is faster than the ODR method for disk packing with imposed bimodal and normal A distribution laws.The results shown in Fig. 14 were obtained for packings of about 50000 disks (the corresponding domain size is different for all the distributions and polydispersities), and a similar dependence of the CPU time on the polydispersity was found for packings of about 10000 disks for all three distributions.Packing densities obtained with the ALA for all the distributions (see Fig. 15) and polydispersities were found to be higher than ones provided by the ODR algorithm.Fig. 15 illustrates that in the case of the AL method the densities increase with increasing the polydispersity and can achieve very high values (0.88) whereas they are typically in the range of 0.8 and 0.81 for packings generated by the ODR method.In addition, the study of the packing density as a function of the domain size was conducted.Fig. 16 illustrates the evolution of density as a function of the domain size for the imposed log-normal size distribution law.Therefore, the proposed ALA is more effective and provides higher packing densities than the classical or Optimized Dropping and Rolling methods in the considered 2D context.

#!!
The results reported in this paper demonstrate that the Advancing Layer Algorithm for packing of disks and ellipses provide high packing densities which are achieved thanks to the constructed Height list.The comparison with recently developed Optimized Dropping and Rolling method shows that the proposed packing algorithm has lower computational cost in 2D case for low and moderate variabilities in disks' size.In addition, the generated packing has higher densities than the ones obtained by the ODR method.
The dense packing of ellipses represented by multicircles was constructed for different imposed A, λ and α distributions.The analysis exhibits that ellipse packing densities are lower than the ones obtained for disks and decrease as the ellipticity increases.Nevertheless, density of 0.8 for ellipse packing can be achieved and provides good agreement with imposed area, shape and orientation distribution.
The packing method developed in the present work is based on the approximation of ellipses by a set of circles.Firstly, this facilitates the calculation of contact position and overlapping detection.Secondly, since the obtained packing is in fact a disk packing, the Laguerre-Voronoï Tessellation Method can be easily applied for generating polygonal structures [9].Indeed, the high packing densities achieved thanks to the Height list used in the ALA provide good accuracy in respecting the desired inertia matrix distributions imposed for generating statistically equivalent cell structures (see Fig. [17][18]. The extension of ALA to 3D packing and polygonal structure generation is also currently developed.In this case, triplets of ellipsoids (instead of pairs of ellipses in 2D) are tested for placing each new particle.Some additional optimization techniques need to be implemented to improve the efficiency of the algorithm and to enable parallel processing.

Figure 1 .
Figure 1.Approximation of the ellipse by a set of circles.

Figure 2 .
Figure 2. Flow-chart of the filling algorithm for the first row.

Figure 3 .
Figure 3. Flow-chart of the algorithm for filling the bulk.

Figure 4 .
Figure 4.The first ellipse in the bottom corner of the domain.The spatial orientation angle, α, is shown.

Figure 5 .
Figure 5. Placement of the second ellipse in contact with the first one (the first row).

Figure 6 .
Figure 6.Placement of the ellipse E3 on E1 and E2.The corresponding pairs of circles C q E1 -C g E3 , C p E2 -C f E3 in contact and ellipses heights (H E1 , H E2 and H max ) are shown.

Figure 7 .
Figure 7. Ellipse cells for E 1 and E 3 and candidate domain (the biggest rectangle, see the text for details).

Figure 8 .
Figure 8. Ellipse packing density as a function of aspect ratio with imposed uniform A distribution law and uniform (shown in red) and normal (shown in blue) α distributions.

Figure 10 .
Figure 10.4596 ellipses with uniform A and α distributions and λ = 4.The packing density is 0.7284.

Figure 13 .
Figure 13.Packing of disks with log-normal type size distribution by ALA (5220 disks, density is 0.8493, CPU time is 1.9 s).

Figure 14 .
Figure 14.Comparison between the ALA and the ODR CPU times for different disk size distribution width.

Figure 15 .
Figure 15.Dependence of the packing density provided by ALA on the disk size distribution width for packing about 50000 disks with normal, log-normal and bimodal size distributions.

Figure 16 .
Figure 16.Dependence of packing density obtained by ODR and ALA on the domain size for the imposed log-normal type distribution.