A Minkowski difference‐based advancing front packing technique for generating convex noncircular particles in complex domains

In this work, a Minkowski difference‐based advancing front approach is proposed to generate convex and non‐circular particles in a predefined computational domain. Two specific algorithms are developed to handle the contact conformity of generated particles with the boundaries of the computational domain. The first, called the open form, is used to handle the smooth contact of generated particles with (external) boundaries, while the other, called the closed form, is proposed to handle the internal boundaries of a computational domain with a complex cavity. The Gilbert‐Johnson‐Keerthi (GJK) method is used to efficiently solve the contact detection between the newly generated particle at the front and existing particles. Furthermore, the problem of one‐sided particle lifting, which can cause some defects in the packing structure in existing advancing front methods during packing generation, is highlighted and an effective solution is developed. Several examples of increasing complexity are used to demonstrate the efficiency and applicability of the proposed packing generation approach. The numerical results show that the generated packing is not only more uniform, but also achieves a higher packing density than existing advancing front methods.


INTRODUCTION
In general, three different types of methods have been developed: the dynamical approach, [11][12][13] the constructive approach, [14][15][16][17][18][19][20] and the coupled constructive-dynamic approach. [21][22][23] In the dynamical approach, particles are first generated in a domain and then their positions and/or sizes are dynamically changed under the applied compression/gravity force within the standard DEM framework. However, a large number of numerical steps are often required to obtain a satisfactory packing density, making this approach very time consuming. Additionally, there is a certain amount of overlap between particles due to external/internal forces during packing, resulting in a physically unrealistic pre-stress within the particle assembly.
In the constructive approach for generating particle packings, particles are added sequentially to the domain in accordance with certain geometric rules that are based on the existing particles already generated. This method ensures that there is no overlap between particles. The constructive approach is very efficient, typically several orders of magnitude faster than the dynamical approach. This is because the position of a particle is determined by only a few known geometric constraints. In the coupled constructive-dynamical approach, the packing initially generated by the constructive approach can be further densified by dynamic techniques to achieve a higher density. The efficiency of this coupled approach lies between that of the dynamic approach and that of the constructive approach.
Due to the high efficiency, several constructive techniques have been developed in the last two decades: regular arrangements, sequential inhibition, sedimentation techniques, mesh-based approach and advancing front approach. 24 Among them, the advancing front approach (AFA), which is originally proposed in the seminal work by Feng et al. 16 to pack discs in 2D domains, is very efficient with a linear complexity. For example, on a PC with a single 1 GHz processor, it takes only 3.77 s to pack 1 million discs. Another important advantage of AFA is that the particle size distribution can be user defined, which is very versatile.
There are two versions of the AFA: the closed form and the open form. In the closed form, the initial front is formed first by placing three touching discs in the centre of the domain. Then discs are generated in contact with the other two previous discs in an outward direction. The final particle assembly is obtained by deleting those discs generated outside the external boundary (EB) of the domain. Thus, in the closed form, there may be some gaps between the assembly and the boundaries. To overcome this problem, the open form is developed to include the EB in the generation of the discs. However, the gaps at the top of the domains are not addressed. Alternatively, Bagi 14 built the initial front near the EB and then generated discs inwards. This is called the inward-AFA based on the closed form. However, the central region cannot be completely filled. To reduce the gap near the EB, Dong et al. 25 used the Newton-Raphson iteration method to re-determine the radius and position of the new disc near the EB, which can guarantee that the disc is tangent to the boundary and the other two discs simultaneously. Recently, Xu and Xia 26 have extended the method to the open form to eliminate the gaps at the upper boundary. Overviews of AFA can be found in Löhner and Oñate 18 and Morfa et al. 27 Due to the simple geometric description of circular particles, it has been widely used in particle-based numerical modelling. However, particle shape has some profound effects on the macro-and micromechanical behaviour of granular materials, 28,29 and particles are non-circular in reality. Therefore, how to efficiently generate convex non-circular particle assemblies has been a forefront topic in recent years. While the packing of convex non-circular particles has received considerable attention, the ellipse, being the simplest such particle, has been the focus of much research. However, non-elliptical particle packing remains challenging due to the complex geometry involved, making it difficult to achieve efficient handling. Extending previous packing approaches from circular to non-circular particles is not straightforward. Many existing approaches are based on either the dynamic approach or the coupled constructive-dynamic approach.
As mentioned above, these approaches are very time consuming. Therefore, some researchers have mainly used the constructive approach to generate non-circular particles. For example, Feng et al. 30 extended their closed-form AFA to generate convex polygons and ellipses based on Minkowski Sum. Subsequently, the AFA has been modified and extended to generate convex polygons. 27,31 However, there are three shortcomings: (1) relatively large gaps can be generated near EB, which reduces the quality of particle assemblies; (2) both interior boundary (IB) and EB cannot be well considered; and (3) the one-sided lifting problem during packing (as will be detailed in Section 4.1) using the open form is not reported.
In this work, a Minkowski difference-based advancing front approach is developed to overcome the aforementioned shortcomings of AFA in generating convex non-circular particle packings. Furthermore, careful consideration is given to IB/EB to solve the gap problem between particles and domain boundaries. The fundamentals of the Minkowski difference and its application are described in detail in Section 2. Then, the closed form and the open form are proposed in Sections 3 and 4, respectively. In particular, Section 4 highlights the one-sided lifting problem during packing using the open form, and the corresponding treatment is proposed. Section 5 presents several examples to illustrate the performance and effectiveness of the proposed approach. Finally, conclusions are drawn in Section 6.
It is worth noting that a convex polygon with a certain number of vertices can be used to approximate convex non-circular particles of any shape. Therefore, this work focuses on the packing of convex polygons. In the case of analytically defined convex shapes, such as ellipses and super-quadrics, they can be first discretised into convex polygons, and then the associated packing is reduced to the packing of convex polygons, which can be effectively handled by the proposed procedures.

Minkowski difference and contact state
A geometric shape (2D or 3D) is considered to be the set of all individual points contained in the shape. 4 The Minkowski difference of A and B, denoted as A ⊖ B, is defined here to be the Minkowski sum of A and (−B): where (−B) can be viewed as the symmetric block of B with respect to the origin: The boundaries of A ⊕ B and A ⊖ B are denoted as (A ⊕ B) and (A ⊖ B), respectively. Consider a hexagon A and a quadrilateral B as an example, and define the origin o to coincide with the centroid of B (i.e., c b ). Figure 1A shows the corresponding Minkowski sum and difference of these two blocks. We will focus on the IB and EB in the present study. As a 2D domain can be considered as being formed by a number of connected lines, the Minkowski sum/difference of a line W and the quadrilateral B is also illustrated in Figure 1B  geometric meaning of (A ⊖ B) and (W ⊖ B) are the locus of the centroid of block B (i.e., c b ) sliding along the boundary of block A and line W, respectively, while the boundary of B is kept in contact with the boundary of A or W. The procedures to construct the Minkowski sum of two convex polygons or a polygon and a line are summarised in Algorithms 1 and 2, respectively. Note that in Algorithms 1 and 2, the vertices of a polygon or a line are ordered anti-clockwise with the first vertex having the lowest y coordinate.

Algorithm 2. Minkowski sum of a line and a convex polygon
When two blocks A and B are in contact (in touch or overlap), there must exist at least one point p that is shared by the two objects, p ∈ A; p ∈ B. Hence, the origin o must be enclosed in the Minkowski difference A ⊖ B: Consequently, the contact state of two blocks A and B is equivalent to the following statement: Contact state between two blocks ∶ Figure 2 demonstrates that the relationship between the Minkowski difference of a fixed hexagon A and a moving quadrilateral B, and the origin of the Minkowski difference of the two shapes corresponds to the contact state of the two shapes. 4,33 Consequently, the contact state between the two blocks can be determined by whether or not the origin is enclosed in the corresponding Minkowski difference of the two blocks.
The Minkowski difference of a line and a polygon can also be utilised to evaluate their contact state. However, determining the contact state between a polygon and a line can be done in a simpler manner. For example, if we consider a convex polygon A with k vertices/edges and a line segment W with two vertices (as shown in Figure 3), there are three contact states that can occur between them: separation, touch, and overlap.
The distance d(A, W) between A and W is defined as the minimum distance from the boundary of A to the line W: where v i a is the i th vertex of the polygon A; v s W is the start point of the line W (it can be any point in the line W); n W is the unit vector normal to the line W and points inward towards the packing domain as shown in Figure 3.
Consequently, the contact state between block A and line W is simplified to the following statement: Contact state between block and line ∶ If the EB is convex, the contact state between block A and line W can be determined using Equations (6) and (7). However, if the EB is concave, the polygon should locate in the active area associated with the line before establishing the contact state between block A and line W. As shown in Figure 3B, the active area is defined using two lines (W l and W r ) associated with the line W. Both W l and W r are vertical to W.
The distance d(A, W l ) between A and W l is defined as the minimum distance from the boundary of A to the left line W l : where n W l is the unit vector normal to the line W l formed by rotating n W 90 • anticlockwise as shown in Figure 3B. The distance d(A, W r ) between A and W r is defined as the minimum distance from the boundary of A to the right line W r : where v e W is the end point of the line W; n W l is the unit vector normal to the line W l formed by rotating n W 90 • clockwise as shown in Figure 3B.
The polygon A is located in the active area associated with the line W if the following condition is fulfilled: If the polygon A is located in the active area, the contact state between block A and line W can be determined using Equations (6) and (7). If a polygon is completely on the opposite side of the line, then this polygon is in overlap with the line. In this case, this polygon should be excluded as it lies outside the packing domain. Thus, this situation is considered in the present approach.

Placing a polygon in contact with other polygons/lines based on Minkowski difference
Unlike discs, determining the position of a new polygon in contact with other polygons/lines is not straightforward. As illustrated in Figure  (a) A convex polygon in contact with two convex polygons In Figure 4A, P d and P e are two polygons and the position of a new polygon P n needs to be determined so that it touches polygons P d and P e simultaneously. The two loci of the polygon P n with polygons P d and P e are (P d ⊖P n ) and (P e ⊖P n ), respectively. The two possible centroids of polygon P n are the intersection points of (P d ⊖P n ) and (P e ⊖P n ). This case is used in determining convex particles during outwards-direction packing as discussed in Section 3.3 and associated with the closed form.
(b) A convex polygon in contact with a convex polygon and a line In Figure 4B, P d and W are a polygon and a line respectively, and the position of a new polygon P n needs to be found so that it touches polygon P d and line W simultaneously. The two loci of the polygon P n with P d and W are (P d ⊖P n ) and (W ⊖ P n ), respectively. The intersection points of (P d ⊖P n ) and (W ⊖ P n ) are the potential centroids of polygon P n . This case is often used for determining the second convex particle during the first layer packing or convex particles near IB/EB as discussed in Section 4 and associated with an open form.
(c) A convex polygon in contact with two lines In Figure 4C, W 1 and W 2 are two lines and the position of a new polygon P n needs to be found so that it touches lines W 1 and W 2 simultaneously. The two loci of the polygon P n with W 1 and W 2 are (W 1 ⊖P n ) and (W 2 ⊖P n ), respectively. The intersection points of (W 1 ⊖P n ) and (W 2 ⊖P n ) are the possible centroids of polygon P n . This case is often used for the determination of the first convex particle during the first layer packing in Section 4 and associated with the open form.
Note that convex polygons/lines (i.e., P d , P e , W 1 and W 2 ) are fixed, while the new polygon P n is movable (but only its centroid is unknown). After the two loci of the polygon P n with other polygons/lines are obtained, a linear algorithm proposed by Han et al. 34 is performed to search for convex polygon intersections. The determination of the centroid of a new polygon is given in Algorithm 3.

Algorithm 3. Determination of the position of a new polygon that touches two polygons/lines
It should be noted that Du et al. 31 constructed the Minkowski difference between a polygon and a line or two lines in a more complex way. They first introduced an auxiliary point associated with the line to convert the line to a convex polygon, and then constructed the corresponding Minkowski difference. However, as stated in Section 2.1, the Minkowski difference between a polygon and a line or two lines can be easily and directly constructed. For this reason, the present approach is much more efficient to place new polygons in contact with lines.

Geometric description of convex polygons
Two approaches can be used to define the template or reference shape for generating convex polygons in a geometric domain. The first approach involves random generation, while the second approach involves using realistic irregular particle images. In the latter approach, the boundary points are obtained from imaging and then made convex by extracting the convex hull. Alternatively, a predefined reference polygon can be translated, scaled, and rotated to generate all convex polygons, which can then be sequentially placed in a 2D domain. Suppose that the reference polygon P ref has k vertices/edges. The coordinates of the ith vertex of the nth generated polygon P n (i.e., v i n ) (n = 1, … , NP max ) are computed by: where R n is the rotation matrix for polygon P n ; n is the rotation angle in the interval [0, 2 ]; U n is the scaling matrix and n1 and n2 are two scaling factors; NP max is the maximum number of polygons need to be generated in the domain; c ref is the centroid of the reference polygon P ref ; and v i ref is the coordinates of the ith vertex. Initially, the centroids of all polygons P n (n = 1, … , NP max ) generated from the reference shape are at the origin for two reasons: (1) (−P n ) can be easily obtained for determining the corresponding Minkowski difference; (2) It is convenient for implementation and programming.
After the centroid of polygon P n (i.e., c n ) is determined using the approach proposed in Section 2.2, all vertices of polygon P n (i.e., v i n ) can be updated as follows: F I G U R E 5 New polygon generation based on a reference polygon.
In this study, the two scaling coefficients for the length and width of the generated polygons are set to be the same within the interval [ min , max ], that is, n1 = n2 = n ∈ [ min , max ]. As a result, all polygons have an identical aspect ratio. Also assume that n and n can be determined from a desired probability distribution function (PDF), such as the uniform, Gaussian or lognormal distribution. 19 The process of generating a new polygon is illustrated using a triangle in Figure 5, where a circumscribed disc is also assigned to the polygon (i.e., triangle). The radii of the circumscribed discs for the reference polygon P ref and the generated polygon P n are r ref and n r ref , respectively.
Equations (11)- (13) show that a polygon P n , denoted as P n ( and two shape descriptors ( n , n ). As c ref and v i ref are pre-defined for the reference P ref , and the rotation angle ( n ) and scaling factor ( n ) of the polygon P n also follow a pre-defined PDF, only c n is unknown. Therefore the polygon can be simply denoted as P n (c n ).

3.2
The initial front generation The first step of AFA is to form an initial front, which is a group of polygons or lines in the 2D domain. Then, an active segment is selected, and a new polygon is placed to make contact with the polygons or lines associated with the active front, aiming to achieve a high packing density. After placing the new polygon, the front is advanced, and the process is repeated by selecting a new active front. This procedure continues until the entire 2D domain is filled with non-overlapping polygons. 16 In the present study, both IB and EB are considered (Figure 6), and the corresponding initial fronts are formed in different ways. The domain boundary is first decomposed into straight lines composed of one or more segments and arranged in an anti-clockwise sense. When the shape of IB or EB is circular, elliptical or super-quadric, it can be approximated by a convex polygon through an adaptive sampling algorithm proposed by Han et al. 34 For an IB with super-quadric shape, its surface is expressed as: where a, b and m are three positive numbers, which determine the shape and size of the super-quadric. An ellipse is obtained when m = 2. Many other shapes can be obtained by varying m. When m < 1, the super-quadric is concave, which is not considered here.

Without consideration of IB
When IB is not considered, the first three polygons, denoted as P 1 , P 2 and P 3 , are placed tangent to each other in the centre of the 2D domain, as shown in Figure 7A. The position of the third polygon P 3 can be determined using the method shown in Figure 4A as explained in Section 2.2. However, polygons P 1 , P 2 and P 3 should form an anti-clockwise cycle, as in the case with discs in Feng et al. 16 Then P 1 → P 2 → P 3 → P 1 form the initial front ( Figure 7A).

With consideration of IB
When IB is considered, the initial front is more complex. In this situation, the IB can be considered as a convex polygon. If the IB is formed from a circular solid region as shown in Figure 8A, the region is first approximated as a convex polygon.
To initiate the AFA process with IB, we designate the IB and the first polygon as P 0 and P 1 , respectively. The position of the second polygon P 2 is determined using the method outlined in Section 2.2 to establish an initial front P 0 → P 1 → P 2 → P 0 . Subsequently, the last segment P 2 → P 0 of the front is selected as the first active segment, as depicted in Figure 8A, to determine the third polygon P 3 . It is worth noting that the third polygon P 3 must be positioned on the right-hand side of the segment, as shown in Figure 8B. As more polygons are generated, the front may be extended to P 0 → P 1 → P 2 → · · · → P 0 , as shown in Figure 8C. The initial front is considered complete when a polygon, such as P 22 , touches the first polygon P 1 . At this point, the front becomes P 1 → P 2 → P 3 → · · · → P 22 → P 1 , and the IB P 0 is removed from the front, as shown in Figure 8C.
Because the generated polygons forming the initial front conform to the IB, the problem of tangency between the generated convex non-circular particles and the IB is naturally resolved. This eliminates both the non-smoothness of particle packings and the relatively large gaps that would have been present near the IB as shown in Figure 8C.

New polygon generation and front update
Once the initial front has been formed, the domain can be filled by generating new polygons and incrementally advancing the front outward. The directions of the frontal segments are defined to ensure that any newly generated polygon is placed on the right-hand side when moving along the positive directions of the segments. This procedure of front update for convex noncircular particles is similar to the method used for discs in Feng et al. 16 For completeness, a brief summary of the front update procedure is provided below. The initial front is updated by incrementally advancing it in the outward direction to generate new polygons and fill the domain. The first segment P 1 → P 2 in the initial front without IB is chosen as the current active front. Then, a new polygon P 4 is generated and placed in contact with both P 1 and P 2 , while lying on the right-hand side of the segment P 1 → P 2 , as described in Section 2.2. The initial front is updated by deleting the segment P 1 → P 2 and insert two new segments P 1 → P 4 and P 4 → P 2 to form a new front P 1 → P 4 → P 2 → P 3 → P 1 . The process is then repeated by selecting the segment P 2 → P 3 as the next active front, and generating a new polygon P 5 , as shown in Figure 7B,C.
However, it is possible that a newly generated polygon from an active segment may overlap with other existing polygons, resulting in four different cases as illustrated in Figure 9. To address this issue, a set of front updating rules for discs was proposed by Feng et al. 16 A similar methodology is adopted here for packing convex particles. More details can be found in Feng et al. 16

Treatment of polygons near EB
During the outwards-direction packing, some generated polygons may overlap with the EB. Discarding such polygons can result in large gaps near the EB, leading to non-smooth packing. To address this, the position and size of an overlapping polygon should be adjusted so that it becomes tangent to the EB. Figure 10A illustrates the distance d n between the centroid of the new polygon P n (c n , n ) and the boundary line B t . If d n < u r ref , then the new polygon may overlap with the boundary. To ensure tangency, a new polygon P n (c, ) is generated while maintaining the rotation angle n , but with adjusted size and position c. The adjustment is achieved by solving the following problem: o ∈ (P e ⊖ P n (c, )) ≡ P n touch P e o ∈ (P e ⊖ P n (c, )) ≡ P n touch P d ,  However, there are two extremely cases that need to be considered: (1) If the largest polygon (its scaling coefficient is u ) is in separation with the boundary B t , that is, d (P n (c, u ) , B t ) > 0, another same sized polygon P n (c, u ) is generated and placed near the boundary B t ; and (2) If the smallest polygon (its scaling coefficient is l ) overlaps with the boundary B t , that is, d (P n (c, l ) , B t ) < 0, the polygon is rejected and no polygon is placed at this place.
Except for the above two extremely cases, the contact state between a newly generated polygon P n (with scaling coefficient n ) and the EB line B t can be used to further narrow the range for the lower and upper bounds as follows: To solve Equations (15) and (16), a bisection method for determining the polygon size and position near the EB are implemented in this work, and the algorithm is presented in Algorithm 4. When generating polygons near the boundary EB, it is important to check if the new polygon overlaps with other boundary lines, especially at the corners. As shown in Figure 10B, the new polygon P n ( c ′ , ′ ) may be closest to one boundary (e.g., B l ) before treatment. After determining the size and position of the polygon, it may be tangent to that boundary but overlap with other boundaries (e.g., B t ), as shown in Figure 10B. In this case, a new polygon P n (c n , n ) should be re-generated based on the overlapping boundary. This check for overlapping boundaries should be repeated until the new polygon does not overlap with any other boundaries.
After generating each new polygon, it is necessary to check whether it is near the EB or not. However, since most of the polygons (maybe 90% or more) are not near the EB, a significant amount of computational time can be wasted in this process. To address this issue, a simple yet effective scheme is proposed as shown in Figure 11. The scheme involves introducing two circles: the inscribed circle of EB and the circumscribed circle of IB. If there is no IB, the point c 0 is selected as the geometric centre of the initial triangle ΔP 1 P 2 P 3 (front). Otherwise, the point is chosen as the centre of the circumscribed circle of the IB. The radius of the inscribed circle of EB (i.e., r EB ) is defined as the minimum distance from the point to the exterior boundary W: where v s W i is the start point of the boundary line W i (it also can be any point in the line W i ); n W i is the unit vector normal to the boundary line W i pointing to inner domain; k is the total number of the boundary lines associated with EB.
The treatment for a newly generated polygon near the EB is active if the following criterion is satisfied: where |c 0 c n | is the distance between the centroid of the new polygon P n (i.e., c n ) and the point c 0 . Equation (18) reveals that a new polygon fully contained within the inscribed circle of the EB requires no further checks for proximity to the EB. With this approach, the algorithm can rapidly determine if a polygon is near the EB or not, without having to perform individual polygon checks. Consequently, this strategy significantly enhances the computational efficiency of the algorithm.

Overlap checking of two polygons
It is worth noting that checking for overlaps between polygons is a computationally intensive process, as it must be conducted for each new polygon with respect to all the existing polygons on the front. Therefore, the overlap checking algorithm must be as efficient as possible. In this study, we have utilised four different algorithms for this purpose: 1. The simplest algorithm is the direct search algorithm (DSA), which sequentially checks the intersection of each edge of one polygon with each edge of the other polygon.

Algorithm 5. The GJK algorithm
In order to compare the performance of these four overlap-checking algorithms, a large number of packing cases with up to 5 million convex polygons of different shapes are employed to examine their computational costs and the results are depicted in Figure 12. The efficiency of these four algorithms is in this order: GJK > DSA > LA > MDA when the number of vertices is lower than 5, and GJK > LA > DSA > MDA when the number is greater than 5. Obviously, GJK is the most efficient, while MDA, which is independent of polygonal shapes, is the worst. Most importantly, the packing speed using GJK is about 2 times faster than that of MDA. Thus, GJK is used to speed up the overlapping checking for new generated polygons with the existing polygons on the front. Also note that the actual computational efficiency of each algorithm is sensitive to the number of vertices of polygons.

F I G U R E 13
Rough contact detection between convex polygons using circumscribed disk.
Furthermore, prior to any overlap checking using GJK, a broader contact detection can be used to exclude those newly generated polygons that cannot be in contact with the existing polygons, thereby further improving the performance of the overlap check. As described in Section 3.1, the circumscribed disc is assigned to each polygon. So if the distance between the centroid of two polygons is greater than the sum of the two radii of their circumscribed discs (i.e., |c 1 c 2 | > r 1 + r 2 ), they can be excluded from further overlap checking, as shown in Figure 13

3.5.2
The value of l and u l and u are the lower and upper bounds of that variable that controls the size of polygons near the EB. Thus, the value of u should be moderate. On the other hand, to avoid unwanted small polygons generated near the EB, if the scaling coefficient of a new polygon is less than l , the polygon should be abandoned. In the present study, l = min , u = max . In this situation, all generated polygons in the domain, including those generated near the EB, are within the range [ min , max ]. Other lower and upper bounds can also be chosen by users.
The complete procedure for generating convex noncircular particles in a complex domain with IB/EB using the closed form is given in Algorithm 6. 1: Compute all polygons Pn(o, n) (n=1,…, NPmax) based on the reference polygon Pref using Eqs. (11) and (12) 2: Discretise IB/EB into a set of connected segments as shown in Fig.6 3: Form the initial front if (IB is not considered) then Generate three contacting polygons in the centre of the domain as shown in Fig.7(a) else Generate a array of smoothly placed polygons which conforms to the IB shape using Algorithm 3 as shown in Fig. 8 As B l and B b are two lines, the first polygon P 1 can be generated using the approach shown in Figure 4C in Section 2.2. After P 1 is generated, the front becomes B l → P 1 → B b → B r . Then P 1 → B b is chosen as the next active front. The second polygon P 2 , which is tangent both to the first polygon P l and the boundary line B b , can be generated using the approach shown in Figure 4B in Section 2.2. Accordingly, the front is now updated as The last polygon in the first layer, for instance, polygon P i as marked in brown in Figure 14A, should be in contact with a previous polygon P n and the bottom wall B b , but is in overlap with the right wall B r . In this case, the polygon P i should be placed to touch P n and the left wall B r , as the one marked in yellow in Figure 14A, to complete the first layer of polygons. The corresponding front becomes B l → P 1 → P 2 → · · · → P i → B r and the bottom wall B b has been removed from the front.
However, the first layer generated in this way can cause the one-sided lifting problem as shown in Figures 14A  and 15, with the following two issues: (1) the particles near the right boundary will continue to grow in the subsequent packing steps, as shown in Figure 9B, which further increases the computing cost and decreases the packing density; and (2) a potential discontinuity from the bottom corner of one side to the top corner of the other side may be encountered in the packing of mono-sized polygons without rotation. A similar phenomenon has been reported in disc packings. 26 To overcome the problem, the following scheme is proposed. If a polygon P n (c n , n ) is near the right boundary B r (see Figure 14B), its size and position c should be re-determined by satisfying the following conditions: Equation (19) is similar to Equation (15), except that the polygon P d is replaced by the right wall B r . Then Algorithm 4 can be used to determine the size and position of the polygon near the right boundary. After this treatment, the right-sided lifting problem can be resolved as shown in Figures 14B and 15.

Packing subsequent layers of polygons
The second and third layers can be built upon the first layer following a similar procedure, as shown in Figure 15B,C. However, each newly generated polygon must be checked for possible overlap with all existing polygons on the front. The same four cases as described in Section 3.3 can occur, and the same actions must be taken, as shown in Figure 9. Figure 15B,C also demonstrate that the number of polygons on the front is nearly the same as when using the proposed algorithm to solve the one-sided lifting issue. Successive layers can be placed in a similar fashion. It should be noted that this study takes the top and right boundaries into account. When a polygon is near either boundary, its size and position need to be re-determined using Algorithm 4. Specifically, after treatment of polygons near the right boundary, the packing densities with the final assembly ( Figure 15D) increases from 0.72 to 0.74. This indicates that the packing density can be further improved through treatment. The one-sided lifting issue can create a 'void' or 'hole' in the domain, which reduces the packing density.
Algorithm 7 outlines the complete procedure for generating convex noncircular particles in complex domains using the open form with EB.

EXAMPLES
To demonstrate the effectiveness and versatility of the proposed approach, various particle assemblies have been generated using the closed form and the open form on an Intel Core i5-8300H CPU 2.3 GHz laptop. The approach is programmed using Fortran95 in an in-house code. Figures 16 and 17 illustrate the packing examples for different polygons generated using the closed form and open form respectively. The proposed approach is capable of generating various convex non-circular particle shapes including tri-, quad-, pen-, hex-, hept-agonal and mixed shapes. The packing assemblies show smoothness and minimal gaps near the EB, indicating that the proposed approach is robust and able to handle any type of convex non-circular particles. The IB includes triangular, rectangular, horseshoe and circular shapes, while the EB includes rectangular, circular and  super-elliptical shapes, demonstrating the versatility of the proposed approach to handle different packing domains with either convex or concave interior cavity.
To assess the computational efficiency of the proposed approach, we measured the time required to pack 1 million convex non-circular particles of different shapes into a rectangular domain, taking into account EB. In these simulations, we limited the number of front segments checked during polygon generation to a maximum value of 10 for both forward and backward checks, and the scaling coefficient was uniformly distributed between 1.0 and 1.2, while the particle rotation ranged from 0 to 2π.
The results, presented in Figure 18A, show that although the open form algorithm is generally more complex to programme, it is more efficient than the closed form algorithm. For instance, the open form approach generated 1 million convex rectangles, pentagons, hexagons, heptagons, octagons, and mixed polygons in 3.89, 5.14, 5.59, 6.86, 7.11, and 7.42 s, respectively, using an Intel Core i5-8300H CPU 2.3 GHz laptop. In contrast, Du et al. 31 required 1 s to generate 1101 rectangles with an aspect ratio of 1 using MATLAB on an Intel Core i7-6700k CPU with 16 GB of RAM. Remarkably, our approach generated 1 million rectangles with an aspect ratio of 1 in only 3.89 s, highlighting its superior performance.
We also found that the computational cost is linearly proportional to the number of edges/vertices in a planar domain, provided the number of polygons to be generated is fixed. Notably, the open form algorithm only required 7.11 s to generate 1 million convex octagons, demonstrating its effectiveness. Furthermore, our approach produced a higher packing density between 0.72 and 0.80, as shown in Figure 18B. As the surface of a convex particle becomes more circular, a higher density is achieved. Figure 19 shows the packings of rectangular particles and contact statistics in four square domains. The coordination number is the average number of contacts per particle. It can be seen that the generated packing assembles have almost the same average coordinate number around 4.07, indicating the homogeneity of the packing. There is also a similar particle orientation distribution.
The proposed algorithms can also handle particles with a wide range of size ratios. Figure 20 shows the packings of polygons generated in a square by the closed form algorithm for two different scaling coefficients: 1.2 and 5.0. Tests also indicate that the computation cost increased for a greater size ratio is negligible.

CONCLUSIONS
This work presents a Minkowski difference-based advancing front approach, with both a closed form and an open form, for generating convex noncircular particles with arbitrary shapes in complex domains. The approach places new polygons in contact with existing polygons/lines using the Minkowski difference, which considers both interior and exterior boundaries. To ensure that polygons generated near exterior boundaries are strictly tangent to those boundaries, the bisection method is introduced.
The efficiency of overlap checking between new and existing polygons using four different algorithms is compared, with the Gilbert-Johnson-Keerthi (GJK) algorithm being the most efficient and the explicit construct Minkowski difference being the slowest. The GJK is about two times faster than the explicit construct Minkowski difference. During the packing process using the open form, the one-sided lifting problem is identified and addressed after treating polygons near the right boundary, which further improves computational efficiency and packing density.
Several numerical examples with various interior and exterior boundaries demonstrate the effectiveness of the proposed approach, which can generate 1 million rectangles and octagons in just 3.89 and 7.11 s, respectively, on a laptop with a 2.3 GHz processor. The results show that the proposed approach is highly efficient and robust and can achieve a packing density between 0.72 and 0.80.
It should be noted that the proposed computational framework is limited to generating arbitrary convex particles. However, it can be extended to packings of concave particles with two modifications: (1) When concave particles can be decomposed into convex components, their Minkowski difference can be obtained as the union of individual differences of the convex components; (2) The GJK algorithm is no longer valid for detecting contact states for concave particles, so the DSA contact detection algorithm can be used instead. Additionally, this work is limited to two dimensions, and future work will explore the extension to 3D.