- Email: [email protected]
Generation of triangular mesh with specified size by circle packing
Advances in Engineering Software 38 (2007) 133–142 www.elsevier.com/locate/advengsoft
Generation of triangular mesh with specified size by circle packing W.X. Wang a
a,*
, C.Y. Ming b, S.H. Lo
c
Intelligence Engineering Lab, Institute of Software, Chinese Academy of Sciences, P.O. Box 8718, Beijing 100080, China b Applied Science School, University of Science and Technology Beijing, Beijing 100083, China c Department of Civil Engineering, The University of HongKong, HongKong, China Received 22 February 2005; accepted 11 April 2006 Available online 7 September 2006
Abstract This paper describes an algorithm for the generation of a finite element mesh with a specified element size over an unbound 2D domain using the advancing front circle packing technique. Unlike the conventional frontal method, the procedure does not start from the object boundary but starts from a convenient point within the open domain. As soon as a circle is added to the generation front, triangular elements are directly generated by properly connecting frontal segments with the centre of the new circle. Circles are packed closely and in contact with the existing circles by an iterative procedure according to the specified size control function. In contrast to other mesh generation schemes, the domain boundary is not considered in the process of circle packing, this reduces a lot of geometrical checks for intersection between frontal segments. If the mesh generation of a physical object is required, the object boundary can be introduced. The boundary recovery procedure is fast and robust by tracing neighbours of triangular elements. The finite element mesh generated by circle packing can also be used through a mapping process to produce parametric surface meshes of the required characteristics. The sizes of circles in the pack are controlled by the principal surface curvatures. Five examples are given to show the effectiveness and robustness of mesh generation and the application of circle packing to mesh generation over curved surfaces. Ó 2006 Elsevier Ltd. All rights reserved. Keywords: Mesh generation; Circle packing; Advancing front; Specified size; Boundary recovery; Parametric surface
1. Introduction Finite element mesh generation has evolved rapidly over the last two decades and meshing techniques seem recently to have reached a level of maturity. Various methods have been developed for unstructured mesh generation, among which, the quadtree technique, the Delaunay triangulation and the advancing front approach have been admittedly thought of as three classical methods widely employed by the meshing community today [1–4]. By these methods, the domain boundaries have to be considered first. The circles/spheres packing technique [5–8] has recently attracted great attention in computational geometry as a tool for the placement of nodes within a bound domain. The packing circles provide valid nodal points for a quality Delaunay *
Corresponding author. E-mail address: [email protected] (W.X. Wang).
0965-9978/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.advengsoft.2006.04.006
mesh. For example, Bern et al. [5] and Eppstein [6] triangulated a n-vertex polygonal region (with holes) by packing circles of different radii in it so that no element has angle larger than p/2. Li et al. [9] made use of the advancing front approach to construct quality circle packing, with the first circle found on the boundary, the packing grows gradually towards the interior of the domain. Feng et al. [10] extended the idea of advancing front approach to construct packing of circles of different radii with considerations for application to discrete element methods (DEM). However, no discussion has been given as to how a finite element mesh is generated from the circle packing. It is important to generate meshes with strict limits on element size. Muylle et al. [11] presented a new point creation scheme for generating meshes of an exact element size using the Delaunay triangulation method. In this paper, a new scheme for the generation of finite element meshes with variable sizes over an unbound planar
134
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
domain will be presented. The idea of mesh generation is to connect the centres of packing circles efficiently and robustly by the advancing front method. However, unlike the conventional frontal method, the procedure does not start from the boundary but starts from any convenient point within an open domain. The sequence of construction of the packing circles is determined by the shortest distance from an arbitrary chosen point in such a way that the generation front will be more or less a circle with occasional minor concave parts due to element size variation. As soon as a circle is added to the generation front, finite elements are directly generated by proper connection of the frontal segments with the centre of the new circle. Unlike other mesh generation schemes, the domain boundary is not considered in the process of circle packing, this reduces a lot of geometrical checks for intersection between frontal segments. For the mesh generation over an unbound domain controlled by the distance from a fictitious centre, the structure of the generation front is a simple closed loop circular in shape. The circular generation front will grow in the radial direction when more and more circles are packed together. The sizes of the triangular elements are controlled by packing circles of different sizes as specified by the control function closely together such that as far as possible circles are tangent to one another to minimize gap spaces. The mesh generation can be terminated either by specifying the diameter of the generation front or by the number of packing circles constructed. The advancing front circle packing approach seems to be a very effective tool for such constructions in terms of mesh quality, efficiency and robustness. Finally, if the mesh generation of a physical object is required, the object boundary can be introduced by any known well-established procedure [12,13]. As the neighbouring relationship of elements has already been established in the circle packing process, a boundary recovery procedure by the method of neighbour tracing is very efficient and robust. Details on such a boundary insertion procedure by neighbour tracing will be discussed in Section 3. In case mesh generation over a curved surface is needed, the finite element mesh generated by circle packing can also be used through a mapping process to produce parametric surface meshes of the required characteristics. The sizes of circles in the pack are controlled by the principal surface curvatures. Examples of various characteristics and practical meshing problems are given to illustrate the combined use of the unbound meshing scheme, the boundary recovery procedure by neighbour tracing, and the application of circle packing to mesh generation over curved surfaces.
within an open domain. Usually the convenient point is the centre of the object. The initial front consists of three tangent circles, which will expand towards the exterior and has to be updated when a new circle is added. Elements are subsequently generated when the centres of the circles are connected to form triangles. The data structure and the rules of circle packing are described in the following sections. 2.1. Data structure The structure of the generation front is very simple and is a closed loop of line segments (circles). As a result, in the mesh generation process, only those circles on the generation front need to be considered. In the implementation, the generation front is represented by a double list of circles, for each circle the centre, the radius and the distance from the origin are stored. The process of front updating is just to apply the operations of inserting (deleting) a point to (from) the generation front. It is interesting to note that the circles along the frontal loop are closely packed and nearly tangent to each other. 2.2. Three criteria for circle packing (i) Densest: The circles are to be packed as close to one another as possible. It is best that all circles are packed closely together so that the gaps between them are minimized. However, it is difficult to pack them densely because the specified radius for a given circle do not always fit its neighbours. Hence, not only the centre of the inserting circle but also the size have to be adjusted several times to obtain a dense packing. (ii) No overlapping: There is no or little overlapping between any two circles. It is best that any new circle added on the generation front does not overlap with existing circles. It is helpful to control element size. The overlapping of packing circles leads to perhaps not only bad mesh quality but also faulty mesh connection such as mesh intersection, overlapping or holes. (iii) Nearest: A new circle is always generated at the location closest to the origin. The front nearest to the origin is often the most concave part. Packing new circles at the concave part can fill gaps, which reduce the chance of forming large holes/voids. By packing circles nearest to the origin, the shape of the generation front is basically convex like a circle with minor concave parts.
2. Circle packing algorithm 2.3. Circle packing and mesh generation The idea of this mesh generation procedure is to connect the centres of packing circles with specified sizes by the advancing front approach. Unlike the conventional advancing front method, the procedure does not start from the object boundary but starts from any convenient point
2.3.1. The initial pack and the coefficient b To simplify, the first three circles, denoted by C1, C2 and C3, are placed tangent to each other around the origin O(0, 0) as shown in Fig. 1a. The sizes of the three circles
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
C1
C1
Ci O
O C3
C2
C3
C2
Fig. 1. Initial front and progression of the generation front: (a) initial front; (b) update front after a new circle is added.
are determined according to a given size control function. Triangle C1C2C3 is the first element of the mesh. The initial front {C1, C2, C3} including information of the circle centre, the radius and the distance to the origin is recorded in a double directed linked list {C1 M C2 M C3 M C1}. The parameter bC1 C2 described below shows how close two circles are 1 bC1 C2 ¼ min k; ; ð1Þ k with k¼
r1 þ r2 ; d
where r1 and r2 are the radii of the circles C1 and C2 respectively, and d is the distance between the two centres of circles C1 and C2. The bigger bC1 C2 ð0 6 bC1 C2 6 1Þ, the closer are the two circles. In case bC1 C2 ¼ 1, i.e. d = r1 + r2, the two circles are touching each other. 2.3.2. Fitting a circle to the existing pack Suppose that C1 and C2 are two neighbouring circles on the generation front. P is an arbitrary point left of C1C2 as shown in Fig. 2. A circle starting from point P can be packed to C1 and C2 following the steps: (i) Calculate the radius r at point P according to the size control function.
a
b
P
P
r
ri
r2 C2
Pi
r1 C1
C2
C1
Fig. 2. Move P to closer C1C2: (a) circle P is far away C1C2; (b) P moves to C1C2.
135
(ii) Calculate b ðb ¼ bC1 P bC2 P Þ at point P using Eq. (1). (iii) If b > k (k 0.95), the circle at P can be fitted to C1C2, end; otherwise. (iv) Calculate bx at point (P + Dx), by at point (P + Dy), Dbx (Dbx = bx b) and Dby (Dby = by b) for arbitrarily small Dx and Dy. Db x (v) Move P to P* along the direction Db ; y for a Dx Dy distance (r1 + r2)d%, where d takes progressively descending values of 10, 8, 6, 4, 2 and 1 in the iterative process. (vi) Update P with P* and go to (i). According to the above procedure, the circle P moves towards C1 and C2 gradually as shown in Fig. 2, and takes up position Pi as the final position when the closeness criteria 2.3.2 (iii) is satisfied. The circle Pi is in contact with circles C1 and C2 to within the specified tolerance k. In the iterative process, the size of circle P is specified by control functions. As the position of initial circle P would significantly affect the rate of convergence in fitting P to C1 and C2, it is important to have an accurate estimate of the initial position of point P. The following gives a brief description of how the initial position of point P is determined. Suppose rp = (r1 + r2)/2, where rp denotes the radius of circle P, r1 and r2 are the radii of the circles C1 and C2, respectively. Then, the centre of circle P is calculated based on rp, C1, r1, C2, r2 and the tangent conditions to circles C1 and C2. 2.3.3. Checking intersections and mesh generation Let Cm be the nearest circle to the origin on the front and Cm+1 be the circle on the left of Cm. Following the procedure as discussed in Section 2.3.2, Ci is the circle to be inserted into the pack at the site CmCm+1. The coefficient fij = dij (ri + rj) indicates if Ci intersects with the generation front, where dij is the distance between the centres of Ci and Cj, ri and rj are the radii of the two circles Ci and Cj respectively, and Cj is a circle on the front with j 5 m and j 5 m + 1. If circles Ci and Cj are overlapped, fij is less than zero; and if circles Ci and Cj are far apart, fij is greater than zero. Hence, coefficient fij measures the relative positions of Ci and Cj. Usually, fij > 0, that is, Ci is not overlapping with any circles Cj on the front, a new element CmCiCm+1 is directly generated and the generation front has to be updated to include Ci such that { . . . Cm M Cm+1 . . . } becomes { . . . Cm M Ci M Cm+1 . . . }. As shown in Fig. 1b, a new element C1CiC2 is generated at the site C1C2. However, when a large circle is required at relatively small insertion site, and as a result Ci penetrates deeply into some neighbouring circles as shown in Fig. 3a. To remedy this situation, we have to push Ci away from the line segment CmCm+1 until Ci is tangent with the two outermost circles and no overlapping with other circles on the front as shown in Fig. 3b. This pushing-out operation is discussed in detail as follows:
136
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
C m+ p C m−q
C m+ p ' ri Ci
ρi
C m +1
Ci
Cm − q '
Cm
C m +1
b C m+ p
Cm
Ci C m−q
b C m+ p C m−q
C m+ p '
C m +1 Cm Fig. 3. Serious overlapping of frontal circles: (a) insertion of large circle Ci at a concave site; (b) push Ci away from the segment CmCm+1.
(i) Calculate Ci touching Cm and Cm+1 according to the procedure 2.3.2. (ii) A search on the neighbouring circles of CmCm+1 is conducted to find out if any circle near CmCm+1 intersects with Ci. (iii) If Ci does not intersect the generation front, element CmCm+1Ci is directly generated, and packing circle Ci to circles Cm and Cm+1 terminates. (iv) Suppose Ci penetrates into Cj, such that fij < 0, the process of pushing-out Ci is invoked. (v) If the circle Cj intersecting with Ci is on the left-hand side of CmCm+1, then j = m + p with p P 2. On the other hand, if the circle Cj is on the right-hand side, then j = m q with q P 1. (vi) Updating the supporting line segment CmCm+1 for the construction of Ci, such that Ci is now generated from {Cm+p, Cmq}. It is noted that either p or q will be updated within each cycle of this pushing-out operation. (vii) Go to step (i). Here, the problem under consideration is the triangulation of polygon {Cm, Cm+1, . . . , Cm+p, Ci, Cmq . . . }. From the construction process of the polygon, it is known that no circle of radius greater than ri can be filled at the interior of the polygon without intersection with other circles. However, if the radius is not required to be in strict, then a circle with radius qi smaller than ri can be constructed within the polygon to produce triangular elements with better quality. From the polygon, take any three points
Ci
Cm − q '
C m +1 Cm Fig. 4. Adjustment of Ci and triangulation of n-polygon: (a) fitting Ci in a n-polygon; (b) triangulation of n-polygon.
to construct the new Ci. Ci is defined as the circle which is tangent to the circles associated with the three chosen points. The new Ci is selected as the largest circle formed tangent to three circles on the polygon. This new Ci will be accepted if its radius qi does not deviate much from the required value of ri. In this case, a new circle Ci can be formed, the polygon has to be updated to fC m ; C mþ1 ; . . . ; C mþp0 ; C i ; C mq0 ; . . .g in which 1 6 p 0 6 p and 0 6 q 0 6 q as shown in Fig. 4a and b. The updated polygon can now be triangulated by any standard procedure [3]. From our experience, even in the most extreme cases, n (the maximum of p and q) is less than or equal to 9; and in a majority of cases n is only 5 or 6. As for the generation front, it could be updated by including two new segments C mq0 C i and C i C mþp0 , and at the same time segments from C mq0 to Cm, CmCm+1 and from Cm+1 to C mþp0 have to be deleted from the front. 2.4. Termination procedure There are several ways to control the termination of the advancing front mesh generation procedure. One is the number of packing circles. The other is the shortest distance from the front to the origin. If the number of circles
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
packed is more than a given number or if the shortest distance to the origin is larger than a specified value, the procedure terminates. We also terminate a section of the front once a generated circle lies seriously outside the boundary of a domain. This would make sense when generating a mesh for domains with an extreme aspect ratio. 2.5. Efficiency and time complexity
(seconds)
There are three main operations in the circle packing process. The first operation is to determine the circle insertion site nearest to the origin. The second operation is to determine the location, size of the circle to be packed at the insertion site. The third operation is to check if there are any intersections between the inserted circle and the existing circles on the generation front. As the distances of the circles from the origin are stored when they are created, the searching for the nearest circle is only a simple comparison. The computational cost grows as the number of circles on the front increases. However, the number of circles on the front increases varies slowly, because some circles on the front will be deleted when a new circle is added. Further investigation into the searching process, we found that the searching time could be much reduced by restricting the distance check to a few circles near the current insertion site. In the present implementation, only five circles to the left and five to the right of the current insertion site are checked to locate the nearest point from the origin for the next insertion. The meshes generated are essentially the same compared to the scheme in which all the circles on the front are checked for the absolute minimum distance point. The result, however, is a great reduction of mesh generation time, and a linear relationship between the CPU time and the number of elements generated can be found as shown in Fig. 5. The time of iteration to determine the position and the size of the inserting circle is controlled by two major factors. The first factor is the initial position of the circle where iteration starts, and the second is the precision adopted for the final converged position of the inserting circle. A scheme has been proposed in Section 2.3.2 to estimate the initial position of the circle to be inserted. Following this simple procedure, the number of iterations has been greatly reduced, and the average number of iterations
16 14 12 10 8 6 4 2 0
137
Table 1 CPU time for packing circles with random radius from 1 to 100 Circles
Elements
Seconds
10,000 20,000 30,000 40,000 50,000 60,000 70,000 80,000 90,000 100,000
19,610 39,440 59,326 79,212 99,119 119,024 138,930 158,873 178,818 198,720
1.496 3.063 4.516 6.033 7.553 9.133 10.675 12.152 13.695 15.112
for various size control functions used in the examples is given in Table 2. A tighter control on the condition of touching between the inserting circle and the existing circles would mean more iteration. However, a tolerance of 0.95 or above adopted in the examples is usually sufficient, and circles are virtually touching each other in the pack as will be shown in the case examples in Section 4. The third operation is to check for overlapping between the proposed circle and the existing circles on the front. If all the circles on the front are checked, the cost could be very expensive. However, because of the very simple ring structure of the front, the number of frontal circles involved in the overlapping check could be reduced to a fixed number M, resulting in a linear complexity in the checking process. This number M is related to the maximum ratio ri/rj, where, ri and rj are the radii of two neighbouring circles. From the experience of some artificial extreme cases, in which the ratio of radii of neighbouring circles can be as large as 1000, we found that 7 circles away from the inserting site had to be checked to ensure no circles overlapping. In all the examples, we set M = 10, that is checking of overlapping was conducted for 10 circles on the left and 10 circles on the right of the insertion site (minimum distance from the origin). Indeed, 90% of the insertion cases, there was no overlapping or overlapping with only one neighbouring circle. The justification for the use of a smaller number M = 10 is that the front is more or less in a convex shape of a circular ring with minor concave parts, hence the new circle does not overlap with circles far away from the insertion point. Table 1 shows the CPU times of a series of examples of circle packing and mesh generation with random radius from 1 to 100. The meshes were generated on a PC with CPU speed of 2.4 GHz and 512 MB RAM using a VC++6.0 compiler. From Fig. 5, the CPU time increases slowly with the number of circles packed, and a quasi-linear relationship can be observed. 3. Boundary insertion
1
2
3
4
5 6 7 (×10000) circles
8
Fig. 5. The CPU time plot for circle packing.
9
10
On a two-dimensional domain, the boundary of an object can be defined as a series of line segments, which could be introduced to the unbound mesh one by one. The elements outside the object boundary or, in case of
138
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
openings, inside the interior boundaries can then be deleted to recover the meshed object. The boundary recovery procedure can be done in five steps as follows: (1) Identify the boundary segments to be inserted. (2) Insert the boundary segments into the background mesh and update the mesh topology locally. (3) Recover the boundary edge segments by element swap. (4) Improve the mesh quality locally along the boundary. (5) Remove the elements not belonging to the object outside of boundaries. As a typical object boundary is a loop of ordered line segments, the most effective way of introducing the loop is to insert the line segment one by one following its natural order. As the neighbouring relationship of the triangular mesh has already been established in the generation process, the triangular elements intersected by the boundary segments can be determined through a neighbour tracing process [14,15] without searching. Whenever a node is inserted into the background mesh, mesh topology has to be modified to include the new node. If the node is at the centre of the element, the element can be divided into three triangles. If the node is close to a vertex of the triangle, the vertex of the triangle can be moved to the position of the node. And if the node falls on an edge sharing by two triangles, simply divide this edge into two to include the newly inserted node. After node insertion, the two ending nodes are mesh vertices and all the elements cut through by the line segment are known through the neighbour tracing process. By means of the neighbouring relationship, an efficient edge recovery scheme can be devised through a systematic use of edge swaps [15]. The entire boundary loop can be recovered by recovering line segments in a sequential manner until we come back to the starting node. After the introduction of the first node, the boundary recovery is fairly automatic as no search is needed throughout the whole process, and neighbouring relationship has been used extensively to facilitate the determination of intersecting elements and edge swaps. Mesh optimization along the object boundary could be applied to improve the element quality as elements are disturbed in the boundary recovery process. Here again, the neighbouring relationship can be used to speed up many basic smoothing operations, such as edge swapping, edge collapsing, node deleting and node re-positioning. Finally, the meshed object can be obtained by removing all the elements not belonging to the object. This can be easily done with the neighbouring relationship. We can start with any element inside the object. This seed element can grow into a larger patch by including all its neighbouring elements. However, in the process, neighbouring elements separated by a boundary segment need not to be included. The meshed object is given by the patch of elements bound by the boundary segments of the object. With
the help of the neighbouring relationships, the recovery process is fast and robust, and no search or intersection check is required. 4. Examples and applications In this section, we provide five examples of mesh generation of specified element size over a 2D unbound domain by circle packing. In the first of the examples, the element size is set randomly. In the next two examples, the element size is controlled by a distance function. In the last two examples, the element size of 2D mesh is controlled by surface curvatures and surface mesh is generated through a mapping process. Table 2 shows the statistics of the example meshes. N and M are respectively the number of nodes and the number of elements in the mesh. Max1 is the maximum deviation of element size from the required size as specified by the size control function. Avg1 is the average of the deviation of element size from the required size. Max2 and Avg2 are respectively the maximum and the average of the ratio of neighbouring element size. Min3 and Avg3 are respectively the minimum and the average of element shape quality as defined in Ref. [3]. It is noted that the element shape quality is calculated from the raw meshes without mesh optimization. Example 1 shows a rather difficult case of generating finite element mesh of random size with a range up to 1– 1000. In other words, the difference in size between neighbouring elements can be as large as 1000. Under this extreme condition, the deviation from the required size is about 15%(Avg1), the maximum size ratio between neighTable 2 The statistics of the five examples Random size
‘‘a’’ Shape
Structure
Surface 1
Surface 2
N M Max1 Avg1 Max2 Avg2 Min3 Avg3 Iteration CPU (s)
1000 2552 8864 8138 40,000 1898 5058 17,697 16,200 79,391 0.998917 0.998713 0.999211 0.989928 0.999281 0.151276 0.062647 0.056686 0.031940 0.024702 327.3385 3.918444 2.961564 4.055701 4.112251 4.222539 1.757723 1.395665 1.145573 1.105738 0.107286 0.401702 0.486297 0.501469 0.618883 0.823691 0.878948 0.921719 0.946467 0.951014 5 3 4 2 3 0.152 0.379 1.327 1.216 6.082 Max1: max qri ; qrii 1; ði ¼ 1; 2; . . . ; N Þ, i i P h Avg1: N1 Ni¼1 max qri ; qrii 1 , i r Max2: max rrji ; rji ði; j ¼ 1; 2; . . . ; N e ; i; j are neighbouring circlesÞ, P r Avg2: N1e Ni¼1 max rrji ; rji , pffiffiffi jABACji Min3: minðai Þ; ai ¼ 2 3 jABj2 þjBCj2 þjCAj 2 ði ¼ 1; 2; . . . ; MÞ, i i i PM 1 Avg3: M i¼1 ðai Þ ði ¼ 1; 2; . . . ; MÞ, where N is the number of nodes, M is the number of elements, Ne is the number of edges, qi is the specified radius of circle i, ri is the actual radius of circle i, rj is the radius of the neighbouring circle of i, A, B and C are three vertices of triangle i.
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
139
b
b
c
Fig. 6. Circle packing and mesh generation for random size control function: (a) packing of circles of random sizes; (b) mesh of finite elements of random sizes.
bouring elements is 327 times (Max2); however, the average shape quality is fairly good at a value of 0.82(Avg3) and the average iterative time is only 5. Fig. 6a and b shows the circle packing and the triangular mesh for this random size control function. Fig. 7a and b shows the circle packing and the corresponding mesh controlled by a curve with shape ‘‘a’’, in which the average deviation of element size from the required size specified by the distance function is about 6%(Avg1). A magnified view of the left of Fig. 7a and b
Fig. 7. Circle packing and mesh generation controlled by a curve: (a) packing of circles controlled by a curve; (b) mesh generation of (a); (c) a magnified view of a portion of (a) and (b).
is shown in Fig. 7c. The rapid size change between two neighbouring circles is 1.76(Avg2). The average iterative time is relatively little only 3. Even without mesh optimization, the elements are of high quality with an average avalue of 0.88(Avg3).
140
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
Fig. 8. Circle packing and mesh generation for a building structure: (a) a transfer plate structure of a tall building; (b) circle packing and mesh generation controlled by boundaries and constraints of (a); (c) insertion of boundary segments; (d) discarding elements outside object boundary.
As shown in Fig. 8a, the boundaries and constraints in a transfer plate structure of a tall building specify the variable size by a distance function. Setting an origin roughly at the geometrical centre of the object, the circle packing technique is applied to generate closely packed circles of specified sizes to cover the entire object as shown in Fig. 8b. A triangular mesh over the unbound domain covering the transfer plate structure is obtained by connecting up the circle centres as shown in Fig. 8b. The domain boundaries and other meshing constraints are introduced and are marked with thicker lines as shown in Fig. 8c. Upon removal of elements outside the object boundary, the meshed object is recovered as shown in Fig. 8d. From Table 2, the average a-value of the mesh is 0.92, and the average deviation from the specified size is about 5%. The next examples are about the application of planar meshing to curved surfaces. In general, the surface can be represented by a bivariate mapping of the form [16,17] T
ðx; y; zÞ ¼ rðu; vÞ
ð2Þ
Fig. 9a and b shows the circle packing and the corresponding triangular mesh in which the size of circles are controlled by the principal curvature of parametric surface. The parametric surface with two distinct peaks is defined by x¼u y¼v
ð3Þ 2
2
z ¼ ðu þ 3v Þe
1u2 v2
In this example, the average deviation of element size from the required size is about 3%(Avg1). The average iterative time is only 2. Even without mesh optimization, the elements are of high quality with an average a-value of 0.95(Avg3). Fig. 9c shows the surface mesh mapped from the mesh of Fig. 9b. Fig. 10a and b shows the circle packing and the corresponding triangular mesh of a wavy curved surface. The parametric surface is defined [18] by
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
Fig. 9. Parametric mesh generation by circles packing: (a) packing circles controlled on curvature of parametric surface; (b) triangular mesh corresponding to (a); (c) surface mesh mapped form (b).
141
Fig. 10. Parametric mesh generation by circle packing: (a) packing circles controlled on curvature of parametric surface; (b) triangular mesh corresponding to (a); (c) surface mesh mapped form (b).
x ¼ u3 þ 10u y ¼ v3 þ 10v z ¼ 100 sin u cos v
ð4Þ
In this last example, the average deviation of element size from the required size is very little, only about 2%(Avg1). The average iterative time is equal to 3. Fig. 10c shows the surface mesh generated from mapping the mesh of
Fig. 10b. In 3D space, the boundary of surface can be conveniently defined through a surface intersection process [15]. 5. Conclusions and discussions In this paper, a procedure for the generation of a triangular mesh of variable element size in compliance with a
142
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
specified function over an unbound domain is presented. Based on the advancing front circle packing technique, with a control of the generation front by the distance from a fictitious centre, the patch of elements will grow in the radial direction as circles are packed together on the generation front. In the mesh generation process, circles of size consistent with the distance function or metric tensor are inserted on the generation front one by one. The proposed iterative procedure allows circles to be packed tightly on the generation front. Very high quality meshes can be obtained if a smooth size control function is used such that the ratio of radii between adjacent circles is not excessive. As in most cases, there is no circle overlapping or only minor adjustment is needed in the insertion process, the quality of the mesh is only governed by the length requirement as specified by the size control function. As usually very good triangular elements are generated by just connecting up the packing circles, in all the examples shown, no optimization on mesh quality has been done on the raw meshes. In case the mesh generation of a physical object is required, what is needed is to introduce the boundary segments into the background mesh and remove those elements not belonging to the object. With illustration by practical examples, an efficient boundary segment insertion process by tracing of neighbours has also been described in Section 3. The packing process is robust as no difficulty is encountered in the generation of meshes of random size with a length between adjacent elements up to 1:1000. In fact, as long as a tight packing of circles is maintained in the insertion process, there is no need to check for intersection as the generation front controlled by the distance from a fictitious centre will always be a closed loop with a simple geometrical shape of a circle. The circle packing by advancing front without checking for intersection is very fast as it is a local process within a few circles at the insertion site on the generation front. The relationship between CPU time and the number of elements (circles) generated is virtually linear. A useful application of the mesh generation is to generate parametric curved surface mesh of various mesh characteristics. The metric tensor required for circle packing can be derived from the surface curvatures and other constraints. After packing circles as specified by the metric tensor, the 2D triangular mesh is mapping to 3D space and a 3D surface mesh is generated. In case the boundary of the meshed surface is needed, it can be conveniently introduce through the process of intersection of discretized surfaces. Finally, the extension of the present scheme for mesh generation by the advancing front sphere packing technique over three-dimensional domain is now under serious investigation.
Acknowledgement The work described in this paper was supported by a grant from National Key Basic Research and Development Program, No. 2002CB312103 and the Research Grants Council of the Hong Kong SAR Government on the project HKU7117/04E ‘‘Analysis of transfer plate structures using high-performance solid 3D hybrid stress hexahedral elements’’. References [1] Watson DF. Computing the n-dimensional Delaunay tessellation with application to voronoi polytopes. Comput J 1981;24(2):167–72. [2] Yerry MA, Shephard MS. A modified quadtree approach to finite element mesh generation. IEEE Comput Graph Appl 1983;3(1): 39–46. [3] Lo SH. A new mesh generation scheme for arbitrary planar domains. Int J Numer Methods Eng 1985;21(8):1403–26. [4] Frey PJ, George PL. Mesh generation: application to finite elements. Oxford, UK: Hermes Science; 2000. [5] Bern M, Mitchell S, Ruppert J. Linear-size nonobtuse triangulation of polygons. Discrete Comput Geomet 1995;14(4):411–28. [6] Eppstein D. Faster circles packing with application to nonobtuse triangulation. Int J Comput Geomet Appl 1997;7(5):485–91. [7] Bern M, Eppstein D. Quadrilateral meshing by circle packing. Int J Comput Geomet Appl 2000;10(4):347–60. [8] Shimada K, Yamada A, Itoh T. Anisotropic triangulation of parametric surfaces via close packing of ellipsoids. Int J Comput Geomet Appl 2000;10(4):417–40. [9] Li XY, Teng SH, Ungor A. Biting: advancing front meets sphere packing. Int J Numer Methods Eng 2000;49(1–4):61–81. [10] Feng YT, Han K, Owen DRJ. Filling domains with disks: an advancing front approach. Int J Numer Methods Eng 2003;56(5): 699–713. [11] Muylle J, Ivanyi P, Topping BHV. A new point creation scheme for uniform Delaunay triangulation. Eng Comput 2002;19(5–6):707–35. [12] Owen SJ, Staten ML, Canann SA, Saigal S. Q-Morph: An indirect approach to advancing front quad meshing. Int J Numer Methods Eng 1999;44(9):1317–40. [13] George PL, Hecht F, Saltel E. Automatic mesh generator with specified boundary. Comput Methods Appl Mech Eng 1991;92(3): 269–88. [14] Lo SH, Wang WX. An algorithm for the intersection of quadrilateral surfaces by tracing of neighbours. Comput Methods Appl Mech Eng 2003;192(20–21):2319–38. [15] Lo SH, Wang WX. Finite element mesh generation over intersecting curved surfaces by tracing of neighbours. Finite Elem Anal Des 2005;41(4):351–70. [16] Lee CK. On curvature element-size control in metric surface mesh generation. Int J Numer Methods Eng 2001;50(4):787–807. [17] Lee CK. Automatic metric 3D surface mesh generation using subdivision surface geometrical model. Part 1: Construction of underlying geometrical model. Int J Numer Methods Eng 2003;56(11):1593–614. [18] Borouchaki H, Laug P, George PL. About parametric surface meshing. In: Second symposium on trends in unstructured mesh generation, USNCCM’99, University of Colorado, Boulder, CO, USA, 1999.
Generation of triangular mesh with specified size by circle packing W.X. Wang a
a,*
, C.Y. Ming b, S.H. Lo
c
Intelligence Engineering Lab, Institute of Software, Chinese Academy of Sciences, P.O. Box 8718, Beijing 100080, China b Applied Science School, University of Science and Technology Beijing, Beijing 100083, China c Department of Civil Engineering, The University of HongKong, HongKong, China Received 22 February 2005; accepted 11 April 2006 Available online 7 September 2006
Abstract This paper describes an algorithm for the generation of a finite element mesh with a specified element size over an unbound 2D domain using the advancing front circle packing technique. Unlike the conventional frontal method, the procedure does not start from the object boundary but starts from a convenient point within the open domain. As soon as a circle is added to the generation front, triangular elements are directly generated by properly connecting frontal segments with the centre of the new circle. Circles are packed closely and in contact with the existing circles by an iterative procedure according to the specified size control function. In contrast to other mesh generation schemes, the domain boundary is not considered in the process of circle packing, this reduces a lot of geometrical checks for intersection between frontal segments. If the mesh generation of a physical object is required, the object boundary can be introduced. The boundary recovery procedure is fast and robust by tracing neighbours of triangular elements. The finite element mesh generated by circle packing can also be used through a mapping process to produce parametric surface meshes of the required characteristics. The sizes of circles in the pack are controlled by the principal surface curvatures. Five examples are given to show the effectiveness and robustness of mesh generation and the application of circle packing to mesh generation over curved surfaces. Ó 2006 Elsevier Ltd. All rights reserved. Keywords: Mesh generation; Circle packing; Advancing front; Specified size; Boundary recovery; Parametric surface
1. Introduction Finite element mesh generation has evolved rapidly over the last two decades and meshing techniques seem recently to have reached a level of maturity. Various methods have been developed for unstructured mesh generation, among which, the quadtree technique, the Delaunay triangulation and the advancing front approach have been admittedly thought of as three classical methods widely employed by the meshing community today [1–4]. By these methods, the domain boundaries have to be considered first. The circles/spheres packing technique [5–8] has recently attracted great attention in computational geometry as a tool for the placement of nodes within a bound domain. The packing circles provide valid nodal points for a quality Delaunay *
Corresponding author. E-mail address: [email protected] (W.X. Wang).
0965-9978/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.advengsoft.2006.04.006
mesh. For example, Bern et al. [5] and Eppstein [6] triangulated a n-vertex polygonal region (with holes) by packing circles of different radii in it so that no element has angle larger than p/2. Li et al. [9] made use of the advancing front approach to construct quality circle packing, with the first circle found on the boundary, the packing grows gradually towards the interior of the domain. Feng et al. [10] extended the idea of advancing front approach to construct packing of circles of different radii with considerations for application to discrete element methods (DEM). However, no discussion has been given as to how a finite element mesh is generated from the circle packing. It is important to generate meshes with strict limits on element size. Muylle et al. [11] presented a new point creation scheme for generating meshes of an exact element size using the Delaunay triangulation method. In this paper, a new scheme for the generation of finite element meshes with variable sizes over an unbound planar
134
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
domain will be presented. The idea of mesh generation is to connect the centres of packing circles efficiently and robustly by the advancing front method. However, unlike the conventional frontal method, the procedure does not start from the boundary but starts from any convenient point within an open domain. The sequence of construction of the packing circles is determined by the shortest distance from an arbitrary chosen point in such a way that the generation front will be more or less a circle with occasional minor concave parts due to element size variation. As soon as a circle is added to the generation front, finite elements are directly generated by proper connection of the frontal segments with the centre of the new circle. Unlike other mesh generation schemes, the domain boundary is not considered in the process of circle packing, this reduces a lot of geometrical checks for intersection between frontal segments. For the mesh generation over an unbound domain controlled by the distance from a fictitious centre, the structure of the generation front is a simple closed loop circular in shape. The circular generation front will grow in the radial direction when more and more circles are packed together. The sizes of the triangular elements are controlled by packing circles of different sizes as specified by the control function closely together such that as far as possible circles are tangent to one another to minimize gap spaces. The mesh generation can be terminated either by specifying the diameter of the generation front or by the number of packing circles constructed. The advancing front circle packing approach seems to be a very effective tool for such constructions in terms of mesh quality, efficiency and robustness. Finally, if the mesh generation of a physical object is required, the object boundary can be introduced by any known well-established procedure [12,13]. As the neighbouring relationship of elements has already been established in the circle packing process, a boundary recovery procedure by the method of neighbour tracing is very efficient and robust. Details on such a boundary insertion procedure by neighbour tracing will be discussed in Section 3. In case mesh generation over a curved surface is needed, the finite element mesh generated by circle packing can also be used through a mapping process to produce parametric surface meshes of the required characteristics. The sizes of circles in the pack are controlled by the principal surface curvatures. Examples of various characteristics and practical meshing problems are given to illustrate the combined use of the unbound meshing scheme, the boundary recovery procedure by neighbour tracing, and the application of circle packing to mesh generation over curved surfaces.
within an open domain. Usually the convenient point is the centre of the object. The initial front consists of three tangent circles, which will expand towards the exterior and has to be updated when a new circle is added. Elements are subsequently generated when the centres of the circles are connected to form triangles. The data structure and the rules of circle packing are described in the following sections. 2.1. Data structure The structure of the generation front is very simple and is a closed loop of line segments (circles). As a result, in the mesh generation process, only those circles on the generation front need to be considered. In the implementation, the generation front is represented by a double list of circles, for each circle the centre, the radius and the distance from the origin are stored. The process of front updating is just to apply the operations of inserting (deleting) a point to (from) the generation front. It is interesting to note that the circles along the frontal loop are closely packed and nearly tangent to each other. 2.2. Three criteria for circle packing (i) Densest: The circles are to be packed as close to one another as possible. It is best that all circles are packed closely together so that the gaps between them are minimized. However, it is difficult to pack them densely because the specified radius for a given circle do not always fit its neighbours. Hence, not only the centre of the inserting circle but also the size have to be adjusted several times to obtain a dense packing. (ii) No overlapping: There is no or little overlapping between any two circles. It is best that any new circle added on the generation front does not overlap with existing circles. It is helpful to control element size. The overlapping of packing circles leads to perhaps not only bad mesh quality but also faulty mesh connection such as mesh intersection, overlapping or holes. (iii) Nearest: A new circle is always generated at the location closest to the origin. The front nearest to the origin is often the most concave part. Packing new circles at the concave part can fill gaps, which reduce the chance of forming large holes/voids. By packing circles nearest to the origin, the shape of the generation front is basically convex like a circle with minor concave parts.
2. Circle packing algorithm 2.3. Circle packing and mesh generation The idea of this mesh generation procedure is to connect the centres of packing circles with specified sizes by the advancing front approach. Unlike the conventional advancing front method, the procedure does not start from the object boundary but starts from any convenient point
2.3.1. The initial pack and the coefficient b To simplify, the first three circles, denoted by C1, C2 and C3, are placed tangent to each other around the origin O(0, 0) as shown in Fig. 1a. The sizes of the three circles
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
C1
C1
Ci O
O C3
C2
C3
C2
Fig. 1. Initial front and progression of the generation front: (a) initial front; (b) update front after a new circle is added.
are determined according to a given size control function. Triangle C1C2C3 is the first element of the mesh. The initial front {C1, C2, C3} including information of the circle centre, the radius and the distance to the origin is recorded in a double directed linked list {C1 M C2 M C3 M C1}. The parameter bC1 C2 described below shows how close two circles are 1 bC1 C2 ¼ min k; ; ð1Þ k with k¼
r1 þ r2 ; d
where r1 and r2 are the radii of the circles C1 and C2 respectively, and d is the distance between the two centres of circles C1 and C2. The bigger bC1 C2 ð0 6 bC1 C2 6 1Þ, the closer are the two circles. In case bC1 C2 ¼ 1, i.e. d = r1 + r2, the two circles are touching each other. 2.3.2. Fitting a circle to the existing pack Suppose that C1 and C2 are two neighbouring circles on the generation front. P is an arbitrary point left of C1C2 as shown in Fig. 2. A circle starting from point P can be packed to C1 and C2 following the steps: (i) Calculate the radius r at point P according to the size control function.
a
b
P
P
r
ri
r2 C2
Pi
r1 C1
C2
C1
Fig. 2. Move P to closer C1C2: (a) circle P is far away C1C2; (b) P moves to C1C2.
135
(ii) Calculate b ðb ¼ bC1 P bC2 P Þ at point P using Eq. (1). (iii) If b > k (k 0.95), the circle at P can be fitted to C1C2, end; otherwise. (iv) Calculate bx at point (P + Dx), by at point (P + Dy), Dbx (Dbx = bx b) and Dby (Dby = by b) for arbitrarily small Dx and Dy. Db x (v) Move P to P* along the direction Db ; y for a Dx Dy distance (r1 + r2)d%, where d takes progressively descending values of 10, 8, 6, 4, 2 and 1 in the iterative process. (vi) Update P with P* and go to (i). According to the above procedure, the circle P moves towards C1 and C2 gradually as shown in Fig. 2, and takes up position Pi as the final position when the closeness criteria 2.3.2 (iii) is satisfied. The circle Pi is in contact with circles C1 and C2 to within the specified tolerance k. In the iterative process, the size of circle P is specified by control functions. As the position of initial circle P would significantly affect the rate of convergence in fitting P to C1 and C2, it is important to have an accurate estimate of the initial position of point P. The following gives a brief description of how the initial position of point P is determined. Suppose rp = (r1 + r2)/2, where rp denotes the radius of circle P, r1 and r2 are the radii of the circles C1 and C2, respectively. Then, the centre of circle P is calculated based on rp, C1, r1, C2, r2 and the tangent conditions to circles C1 and C2. 2.3.3. Checking intersections and mesh generation Let Cm be the nearest circle to the origin on the front and Cm+1 be the circle on the left of Cm. Following the procedure as discussed in Section 2.3.2, Ci is the circle to be inserted into the pack at the site CmCm+1. The coefficient fij = dij (ri + rj) indicates if Ci intersects with the generation front, where dij is the distance between the centres of Ci and Cj, ri and rj are the radii of the two circles Ci and Cj respectively, and Cj is a circle on the front with j 5 m and j 5 m + 1. If circles Ci and Cj are overlapped, fij is less than zero; and if circles Ci and Cj are far apart, fij is greater than zero. Hence, coefficient fij measures the relative positions of Ci and Cj. Usually, fij > 0, that is, Ci is not overlapping with any circles Cj on the front, a new element CmCiCm+1 is directly generated and the generation front has to be updated to include Ci such that { . . . Cm M Cm+1 . . . } becomes { . . . Cm M Ci M Cm+1 . . . }. As shown in Fig. 1b, a new element C1CiC2 is generated at the site C1C2. However, when a large circle is required at relatively small insertion site, and as a result Ci penetrates deeply into some neighbouring circles as shown in Fig. 3a. To remedy this situation, we have to push Ci away from the line segment CmCm+1 until Ci is tangent with the two outermost circles and no overlapping with other circles on the front as shown in Fig. 3b. This pushing-out operation is discussed in detail as follows:
136
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
C m+ p C m−q
C m+ p ' ri Ci
ρi
C m +1
Ci
Cm − q '
Cm
C m +1
b C m+ p
Cm
Ci C m−q
b C m+ p C m−q
C m+ p '
C m +1 Cm Fig. 3. Serious overlapping of frontal circles: (a) insertion of large circle Ci at a concave site; (b) push Ci away from the segment CmCm+1.
(i) Calculate Ci touching Cm and Cm+1 according to the procedure 2.3.2. (ii) A search on the neighbouring circles of CmCm+1 is conducted to find out if any circle near CmCm+1 intersects with Ci. (iii) If Ci does not intersect the generation front, element CmCm+1Ci is directly generated, and packing circle Ci to circles Cm and Cm+1 terminates. (iv) Suppose Ci penetrates into Cj, such that fij < 0, the process of pushing-out Ci is invoked. (v) If the circle Cj intersecting with Ci is on the left-hand side of CmCm+1, then j = m + p with p P 2. On the other hand, if the circle Cj is on the right-hand side, then j = m q with q P 1. (vi) Updating the supporting line segment CmCm+1 for the construction of Ci, such that Ci is now generated from {Cm+p, Cmq}. It is noted that either p or q will be updated within each cycle of this pushing-out operation. (vii) Go to step (i). Here, the problem under consideration is the triangulation of polygon {Cm, Cm+1, . . . , Cm+p, Ci, Cmq . . . }. From the construction process of the polygon, it is known that no circle of radius greater than ri can be filled at the interior of the polygon without intersection with other circles. However, if the radius is not required to be in strict, then a circle with radius qi smaller than ri can be constructed within the polygon to produce triangular elements with better quality. From the polygon, take any three points
Ci
Cm − q '
C m +1 Cm Fig. 4. Adjustment of Ci and triangulation of n-polygon: (a) fitting Ci in a n-polygon; (b) triangulation of n-polygon.
to construct the new Ci. Ci is defined as the circle which is tangent to the circles associated with the three chosen points. The new Ci is selected as the largest circle formed tangent to three circles on the polygon. This new Ci will be accepted if its radius qi does not deviate much from the required value of ri. In this case, a new circle Ci can be formed, the polygon has to be updated to fC m ; C mþ1 ; . . . ; C mþp0 ; C i ; C mq0 ; . . .g in which 1 6 p 0 6 p and 0 6 q 0 6 q as shown in Fig. 4a and b. The updated polygon can now be triangulated by any standard procedure [3]. From our experience, even in the most extreme cases, n (the maximum of p and q) is less than or equal to 9; and in a majority of cases n is only 5 or 6. As for the generation front, it could be updated by including two new segments C mq0 C i and C i C mþp0 , and at the same time segments from C mq0 to Cm, CmCm+1 and from Cm+1 to C mþp0 have to be deleted from the front. 2.4. Termination procedure There are several ways to control the termination of the advancing front mesh generation procedure. One is the number of packing circles. The other is the shortest distance from the front to the origin. If the number of circles
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
packed is more than a given number or if the shortest distance to the origin is larger than a specified value, the procedure terminates. We also terminate a section of the front once a generated circle lies seriously outside the boundary of a domain. This would make sense when generating a mesh for domains with an extreme aspect ratio. 2.5. Efficiency and time complexity
(seconds)
There are three main operations in the circle packing process. The first operation is to determine the circle insertion site nearest to the origin. The second operation is to determine the location, size of the circle to be packed at the insertion site. The third operation is to check if there are any intersections between the inserted circle and the existing circles on the generation front. As the distances of the circles from the origin are stored when they are created, the searching for the nearest circle is only a simple comparison. The computational cost grows as the number of circles on the front increases. However, the number of circles on the front increases varies slowly, because some circles on the front will be deleted when a new circle is added. Further investigation into the searching process, we found that the searching time could be much reduced by restricting the distance check to a few circles near the current insertion site. In the present implementation, only five circles to the left and five to the right of the current insertion site are checked to locate the nearest point from the origin for the next insertion. The meshes generated are essentially the same compared to the scheme in which all the circles on the front are checked for the absolute minimum distance point. The result, however, is a great reduction of mesh generation time, and a linear relationship between the CPU time and the number of elements generated can be found as shown in Fig. 5. The time of iteration to determine the position and the size of the inserting circle is controlled by two major factors. The first factor is the initial position of the circle where iteration starts, and the second is the precision adopted for the final converged position of the inserting circle. A scheme has been proposed in Section 2.3.2 to estimate the initial position of the circle to be inserted. Following this simple procedure, the number of iterations has been greatly reduced, and the average number of iterations
16 14 12 10 8 6 4 2 0
137
Table 1 CPU time for packing circles with random radius from 1 to 100 Circles
Elements
Seconds
10,000 20,000 30,000 40,000 50,000 60,000 70,000 80,000 90,000 100,000
19,610 39,440 59,326 79,212 99,119 119,024 138,930 158,873 178,818 198,720
1.496 3.063 4.516 6.033 7.553 9.133 10.675 12.152 13.695 15.112
for various size control functions used in the examples is given in Table 2. A tighter control on the condition of touching between the inserting circle and the existing circles would mean more iteration. However, a tolerance of 0.95 or above adopted in the examples is usually sufficient, and circles are virtually touching each other in the pack as will be shown in the case examples in Section 4. The third operation is to check for overlapping between the proposed circle and the existing circles on the front. If all the circles on the front are checked, the cost could be very expensive. However, because of the very simple ring structure of the front, the number of frontal circles involved in the overlapping check could be reduced to a fixed number M, resulting in a linear complexity in the checking process. This number M is related to the maximum ratio ri/rj, where, ri and rj are the radii of two neighbouring circles. From the experience of some artificial extreme cases, in which the ratio of radii of neighbouring circles can be as large as 1000, we found that 7 circles away from the inserting site had to be checked to ensure no circles overlapping. In all the examples, we set M = 10, that is checking of overlapping was conducted for 10 circles on the left and 10 circles on the right of the insertion site (minimum distance from the origin). Indeed, 90% of the insertion cases, there was no overlapping or overlapping with only one neighbouring circle. The justification for the use of a smaller number M = 10 is that the front is more or less in a convex shape of a circular ring with minor concave parts, hence the new circle does not overlap with circles far away from the insertion point. Table 1 shows the CPU times of a series of examples of circle packing and mesh generation with random radius from 1 to 100. The meshes were generated on a PC with CPU speed of 2.4 GHz and 512 MB RAM using a VC++6.0 compiler. From Fig. 5, the CPU time increases slowly with the number of circles packed, and a quasi-linear relationship can be observed. 3. Boundary insertion
1
2
3
4
5 6 7 (×10000) circles
8
Fig. 5. The CPU time plot for circle packing.
9
10
On a two-dimensional domain, the boundary of an object can be defined as a series of line segments, which could be introduced to the unbound mesh one by one. The elements outside the object boundary or, in case of
138
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
openings, inside the interior boundaries can then be deleted to recover the meshed object. The boundary recovery procedure can be done in five steps as follows: (1) Identify the boundary segments to be inserted. (2) Insert the boundary segments into the background mesh and update the mesh topology locally. (3) Recover the boundary edge segments by element swap. (4) Improve the mesh quality locally along the boundary. (5) Remove the elements not belonging to the object outside of boundaries. As a typical object boundary is a loop of ordered line segments, the most effective way of introducing the loop is to insert the line segment one by one following its natural order. As the neighbouring relationship of the triangular mesh has already been established in the generation process, the triangular elements intersected by the boundary segments can be determined through a neighbour tracing process [14,15] without searching. Whenever a node is inserted into the background mesh, mesh topology has to be modified to include the new node. If the node is at the centre of the element, the element can be divided into three triangles. If the node is close to a vertex of the triangle, the vertex of the triangle can be moved to the position of the node. And if the node falls on an edge sharing by two triangles, simply divide this edge into two to include the newly inserted node. After node insertion, the two ending nodes are mesh vertices and all the elements cut through by the line segment are known through the neighbour tracing process. By means of the neighbouring relationship, an efficient edge recovery scheme can be devised through a systematic use of edge swaps [15]. The entire boundary loop can be recovered by recovering line segments in a sequential manner until we come back to the starting node. After the introduction of the first node, the boundary recovery is fairly automatic as no search is needed throughout the whole process, and neighbouring relationship has been used extensively to facilitate the determination of intersecting elements and edge swaps. Mesh optimization along the object boundary could be applied to improve the element quality as elements are disturbed in the boundary recovery process. Here again, the neighbouring relationship can be used to speed up many basic smoothing operations, such as edge swapping, edge collapsing, node deleting and node re-positioning. Finally, the meshed object can be obtained by removing all the elements not belonging to the object. This can be easily done with the neighbouring relationship. We can start with any element inside the object. This seed element can grow into a larger patch by including all its neighbouring elements. However, in the process, neighbouring elements separated by a boundary segment need not to be included. The meshed object is given by the patch of elements bound by the boundary segments of the object. With
the help of the neighbouring relationships, the recovery process is fast and robust, and no search or intersection check is required. 4. Examples and applications In this section, we provide five examples of mesh generation of specified element size over a 2D unbound domain by circle packing. In the first of the examples, the element size is set randomly. In the next two examples, the element size is controlled by a distance function. In the last two examples, the element size of 2D mesh is controlled by surface curvatures and surface mesh is generated through a mapping process. Table 2 shows the statistics of the example meshes. N and M are respectively the number of nodes and the number of elements in the mesh. Max1 is the maximum deviation of element size from the required size as specified by the size control function. Avg1 is the average of the deviation of element size from the required size. Max2 and Avg2 are respectively the maximum and the average of the ratio of neighbouring element size. Min3 and Avg3 are respectively the minimum and the average of element shape quality as defined in Ref. [3]. It is noted that the element shape quality is calculated from the raw meshes without mesh optimization. Example 1 shows a rather difficult case of generating finite element mesh of random size with a range up to 1– 1000. In other words, the difference in size between neighbouring elements can be as large as 1000. Under this extreme condition, the deviation from the required size is about 15%(Avg1), the maximum size ratio between neighTable 2 The statistics of the five examples Random size
‘‘a’’ Shape
Structure
Surface 1
Surface 2
N M Max1 Avg1 Max2 Avg2 Min3 Avg3 Iteration CPU (s)
1000 2552 8864 8138 40,000 1898 5058 17,697 16,200 79,391 0.998917 0.998713 0.999211 0.989928 0.999281 0.151276 0.062647 0.056686 0.031940 0.024702 327.3385 3.918444 2.961564 4.055701 4.112251 4.222539 1.757723 1.395665 1.145573 1.105738 0.107286 0.401702 0.486297 0.501469 0.618883 0.823691 0.878948 0.921719 0.946467 0.951014 5 3 4 2 3 0.152 0.379 1.327 1.216 6.082 Max1: max qri ; qrii 1; ði ¼ 1; 2; . . . ; N Þ, i i P h Avg1: N1 Ni¼1 max qri ; qrii 1 , i r Max2: max rrji ; rji ði; j ¼ 1; 2; . . . ; N e ; i; j are neighbouring circlesÞ, P r Avg2: N1e Ni¼1 max rrji ; rji , pffiffiffi jABACji Min3: minðai Þ; ai ¼ 2 3 jABj2 þjBCj2 þjCAj 2 ði ¼ 1; 2; . . . ; MÞ, i i i PM 1 Avg3: M i¼1 ðai Þ ði ¼ 1; 2; . . . ; MÞ, where N is the number of nodes, M is the number of elements, Ne is the number of edges, qi is the specified radius of circle i, ri is the actual radius of circle i, rj is the radius of the neighbouring circle of i, A, B and C are three vertices of triangle i.
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
139
b
b
c
Fig. 6. Circle packing and mesh generation for random size control function: (a) packing of circles of random sizes; (b) mesh of finite elements of random sizes.
bouring elements is 327 times (Max2); however, the average shape quality is fairly good at a value of 0.82(Avg3) and the average iterative time is only 5. Fig. 6a and b shows the circle packing and the triangular mesh for this random size control function. Fig. 7a and b shows the circle packing and the corresponding mesh controlled by a curve with shape ‘‘a’’, in which the average deviation of element size from the required size specified by the distance function is about 6%(Avg1). A magnified view of the left of Fig. 7a and b
Fig. 7. Circle packing and mesh generation controlled by a curve: (a) packing of circles controlled by a curve; (b) mesh generation of (a); (c) a magnified view of a portion of (a) and (b).
is shown in Fig. 7c. The rapid size change between two neighbouring circles is 1.76(Avg2). The average iterative time is relatively little only 3. Even without mesh optimization, the elements are of high quality with an average avalue of 0.88(Avg3).
140
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
Fig. 8. Circle packing and mesh generation for a building structure: (a) a transfer plate structure of a tall building; (b) circle packing and mesh generation controlled by boundaries and constraints of (a); (c) insertion of boundary segments; (d) discarding elements outside object boundary.
As shown in Fig. 8a, the boundaries and constraints in a transfer plate structure of a tall building specify the variable size by a distance function. Setting an origin roughly at the geometrical centre of the object, the circle packing technique is applied to generate closely packed circles of specified sizes to cover the entire object as shown in Fig. 8b. A triangular mesh over the unbound domain covering the transfer plate structure is obtained by connecting up the circle centres as shown in Fig. 8b. The domain boundaries and other meshing constraints are introduced and are marked with thicker lines as shown in Fig. 8c. Upon removal of elements outside the object boundary, the meshed object is recovered as shown in Fig. 8d. From Table 2, the average a-value of the mesh is 0.92, and the average deviation from the specified size is about 5%. The next examples are about the application of planar meshing to curved surfaces. In general, the surface can be represented by a bivariate mapping of the form [16,17] T
ðx; y; zÞ ¼ rðu; vÞ
ð2Þ
Fig. 9a and b shows the circle packing and the corresponding triangular mesh in which the size of circles are controlled by the principal curvature of parametric surface. The parametric surface with two distinct peaks is defined by x¼u y¼v
ð3Þ 2
2
z ¼ ðu þ 3v Þe
1u2 v2
In this example, the average deviation of element size from the required size is about 3%(Avg1). The average iterative time is only 2. Even without mesh optimization, the elements are of high quality with an average a-value of 0.95(Avg3). Fig. 9c shows the surface mesh mapped from the mesh of Fig. 9b. Fig. 10a and b shows the circle packing and the corresponding triangular mesh of a wavy curved surface. The parametric surface is defined [18] by
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
Fig. 9. Parametric mesh generation by circles packing: (a) packing circles controlled on curvature of parametric surface; (b) triangular mesh corresponding to (a); (c) surface mesh mapped form (b).
141
Fig. 10. Parametric mesh generation by circle packing: (a) packing circles controlled on curvature of parametric surface; (b) triangular mesh corresponding to (a); (c) surface mesh mapped form (b).
x ¼ u3 þ 10u y ¼ v3 þ 10v z ¼ 100 sin u cos v
ð4Þ
In this last example, the average deviation of element size from the required size is very little, only about 2%(Avg1). The average iterative time is equal to 3. Fig. 10c shows the surface mesh generated from mapping the mesh of
Fig. 10b. In 3D space, the boundary of surface can be conveniently defined through a surface intersection process [15]. 5. Conclusions and discussions In this paper, a procedure for the generation of a triangular mesh of variable element size in compliance with a
142
W.X. Wang et al. / Advances in Engineering Software 38 (2007) 133–142
specified function over an unbound domain is presented. Based on the advancing front circle packing technique, with a control of the generation front by the distance from a fictitious centre, the patch of elements will grow in the radial direction as circles are packed together on the generation front. In the mesh generation process, circles of size consistent with the distance function or metric tensor are inserted on the generation front one by one. The proposed iterative procedure allows circles to be packed tightly on the generation front. Very high quality meshes can be obtained if a smooth size control function is used such that the ratio of radii between adjacent circles is not excessive. As in most cases, there is no circle overlapping or only minor adjustment is needed in the insertion process, the quality of the mesh is only governed by the length requirement as specified by the size control function. As usually very good triangular elements are generated by just connecting up the packing circles, in all the examples shown, no optimization on mesh quality has been done on the raw meshes. In case the mesh generation of a physical object is required, what is needed is to introduce the boundary segments into the background mesh and remove those elements not belonging to the object. With illustration by practical examples, an efficient boundary segment insertion process by tracing of neighbours has also been described in Section 3. The packing process is robust as no difficulty is encountered in the generation of meshes of random size with a length between adjacent elements up to 1:1000. In fact, as long as a tight packing of circles is maintained in the insertion process, there is no need to check for intersection as the generation front controlled by the distance from a fictitious centre will always be a closed loop with a simple geometrical shape of a circle. The circle packing by advancing front without checking for intersection is very fast as it is a local process within a few circles at the insertion site on the generation front. The relationship between CPU time and the number of elements (circles) generated is virtually linear. A useful application of the mesh generation is to generate parametric curved surface mesh of various mesh characteristics. The metric tensor required for circle packing can be derived from the surface curvatures and other constraints. After packing circles as specified by the metric tensor, the 2D triangular mesh is mapping to 3D space and a 3D surface mesh is generated. In case the boundary of the meshed surface is needed, it can be conveniently introduce through the process of intersection of discretized surfaces. Finally, the extension of the present scheme for mesh generation by the advancing front sphere packing technique over three-dimensional domain is now under serious investigation.
Acknowledgement The work described in this paper was supported by a grant from National Key Basic Research and Development Program, No. 2002CB312103 and the Research Grants Council of the Hong Kong SAR Government on the project HKU7117/04E ‘‘Analysis of transfer plate structures using high-performance solid 3D hybrid stress hexahedral elements’’. References [1] Watson DF. Computing the n-dimensional Delaunay tessellation with application to voronoi polytopes. Comput J 1981;24(2):167–72. [2] Yerry MA, Shephard MS. A modified quadtree approach to finite element mesh generation. IEEE Comput Graph Appl 1983;3(1): 39–46. [3] Lo SH. A new mesh generation scheme for arbitrary planar domains. Int J Numer Methods Eng 1985;21(8):1403–26. [4] Frey PJ, George PL. Mesh generation: application to finite elements. Oxford, UK: Hermes Science; 2000. [5] Bern M, Mitchell S, Ruppert J. Linear-size nonobtuse triangulation of polygons. Discrete Comput Geomet 1995;14(4):411–28. [6] Eppstein D. Faster circles packing with application to nonobtuse triangulation. Int J Comput Geomet Appl 1997;7(5):485–91. [7] Bern M, Eppstein D. Quadrilateral meshing by circle packing. Int J Comput Geomet Appl 2000;10(4):347–60. [8] Shimada K, Yamada A, Itoh T. Anisotropic triangulation of parametric surfaces via close packing of ellipsoids. Int J Comput Geomet Appl 2000;10(4):417–40. [9] Li XY, Teng SH, Ungor A. Biting: advancing front meets sphere packing. Int J Numer Methods Eng 2000;49(1–4):61–81. [10] Feng YT, Han K, Owen DRJ. Filling domains with disks: an advancing front approach. Int J Numer Methods Eng 2003;56(5): 699–713. [11] Muylle J, Ivanyi P, Topping BHV. A new point creation scheme for uniform Delaunay triangulation. Eng Comput 2002;19(5–6):707–35. [12] Owen SJ, Staten ML, Canann SA, Saigal S. Q-Morph: An indirect approach to advancing front quad meshing. Int J Numer Methods Eng 1999;44(9):1317–40. [13] George PL, Hecht F, Saltel E. Automatic mesh generator with specified boundary. Comput Methods Appl Mech Eng 1991;92(3): 269–88. [14] Lo SH, Wang WX. An algorithm for the intersection of quadrilateral surfaces by tracing of neighbours. Comput Methods Appl Mech Eng 2003;192(20–21):2319–38. [15] Lo SH, Wang WX. Finite element mesh generation over intersecting curved surfaces by tracing of neighbours. Finite Elem Anal Des 2005;41(4):351–70. [16] Lee CK. On curvature element-size control in metric surface mesh generation. Int J Numer Methods Eng 2001;50(4):787–807. [17] Lee CK. Automatic metric 3D surface mesh generation using subdivision surface geometrical model. Part 1: Construction of underlying geometrical model. Int J Numer Methods Eng 2003;56(11):1593–614. [18] Borouchaki H, Laug P, George PL. About parametric surface meshing. In: Second symposium on trends in unstructured mesh generation, USNCCM’99, University of Colorado, Boulder, CO, USA, 1999.