An efficient automatic mesh generator for quadrilateral elements implemented using C++

Finite Elements in Analysis and Design 39 (2003) 905 – 930 www.elsevier.com/locate/nel

An e cient automatic mesh generator for quadrilateral elements implemented using C++ M. Bastian, B.Q. Li ∗ Department of Mechanical and Materials Engineering, Washington State University, Pullman, WA 99163, USA Received 17 April 2001; received in revised form 6 March 2002; accepted 14 April 2002

Abstract This paper presents an e cient mesh generation algorithm for creating unstructured quadrilateral meshes for nite element computations. The algorithm is developed by improving upon an existing recursive mesh generation scheme and implemented using the C++ object-oriented language. Exhaustive cases of misshaped elements for which a general remedy is not possible are identied through extensive numerical testing and treated case by case. These improvements make it possible to generate quadrilateral elements for any case without any need for nal mesh resolution by human interference. The implementation of the algorithm using object-oriented methods allows 9exibility in programming and increases the computational e ciency of mesh generation. The use of object data structures, such as line segments and elements, expedites the algorithm implementation. Numerical experience shows that not only is programming simplied when objects are used, but computational speed is also increased. When comparing two nodes for best splitting lines, line segments do not compare nodes on the same segment. This immediately eliminates all nodes on the segment the starting point of the potential splitting line is on. By automatically doing this, a great saving in computing time is achieved. Examples of automatic mesh generation in both simply and multiply connected domains are given and their computational e ciencies are also assessed. ? 2002 Elsevier Science B.V. All rights reserved.

1. Introduction One of the distinctive advantages of the nite element method or other unstructured numerical methods is the generality and ease with which very irregular geometric domains are treated [1,2]. To eAectively implement these numerical methods, automatic mesh generators need to be developed that break a complex domain into various types of elements of desirable shape before computations are performed. Automatic mesh generation algorithms become even more important when adaptive numerical schemes are used for certain classes of problems involving free and moving boundaries ∗

Corresponding author. Tel.: +1-509-335-7386; fax: +1-509-335-4662. E-mail address: [email protected] (B.Q. Li).

0168-874X/02/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 8 7 4 X ( 0 2 ) 0 0 1 3 8 - 5

906

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

[3,4]. For 9uid 9ow and heat transfer calculations using the nite element method, quadrilateral elements are the preferred choice. Good elements will have corners that meet at right angles and are similar in size. In practice, most computational domains do not have simple geometry and require some irregular elements to completely mesh the domains. Meshes generated by automatic generators fall into two categories: structured and unstructured. Structured meshes are produced over a pre-dened block structure and are numbered in a pre-ordered sequence. For a structured mesh, all interior nodes have an equal number of adjacent elements. Unstructured meshes, on the other hand, are often generated over the entire domain and their nodes are numbered arbitrarily. An unstructured mesh allows any number of elements to meet at a single node. While both structured and unstructured meshes can take regular and irregular shapes, the latter often produce element shapes that are more versatile to better ll the computational domains. Several categories of algorithms exist for automatic mesh generation. A recent survey of mesh generation technology has been conducted by Owen [5]. Some algorithms are designed to produce structured meshes, while others are for unstructured mesh generation. Some of the more commonly used algorithm categories include: Paving, Transnite mapping, Delaunay triangulation and other heuristic algorithms. Paving, introduced by Blacker and Stephenson [6], and improved by White and Kinney [7], is a meshing method in which quadrilateral elements are added to a domain, row by row, around the border of the domain until the entire domain is meshed. This technique requires that tests be made for quad intersections and for quads that are su ciently close that a new connecting quad would be very distorted. Even with these tests, some excessively distorted or infeasible elements may still result. This algorithm, like nearly all other quad meshing algorithms, requires clean up and smoothing phases to make the mesh viable for computation. Paving generally produces unstructured meshes. Transnite mapping, also called transport mapping or transformation mapping or transnite interpolation [8], is a method wherein a domain can be transformed from a simply connected geometry such as a circle, to a regular domain, such as square or rectangle. The mesh is generated over these regular domains with appropriately designed grid distributions, and then mapped back to the original geometry. This allows for simple grids to be generated in the transformed domain. This method only works, however, where a transform is available to go from the initial arbitrary domain to a rectangular domain. The application of the method often requires breaking a domain into small larger blocks that can be transformed into regular domains in the transformed space. For example, a circle is often broken ve subdomains, each of which is then transformed into a square or rectangle over which meshes are generated. A multiply connected domain must be sectioned into simply connected regions to warrant the existence of a unique transformation function. The element shapes are not arbitrary and depend on the transformation function used. Transnite mapping produces structured meshes. It is easy to use and program and is still widely used in the numerical analysis industry [9]. Delaunay triangulation is a very popular method for creating triangular meshes in which very good triangles are generally formed [8]. The unique property of a Delaunay triangulation is that the circle formed by any three points in a Delaunay triangle will not encompass any other points, thereby ensuring a well-distributed mesh. There are many methods for constructing this type of meshes. One of the simplest techniques is to connect the neighboring points of a Dirichlet tessellation of the points, as described by Farin [10]. A Delaunay triangulation mesher produces unstructured meshes.

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

907

The eld of mesh generation is still researched very actively, and there are many other new and innovative methods being developed [8]. Some, like the ones above, are general purpose while others are written for specic cases. One algorithm that has been developed recently by Sarrate and Huerta [11] that breaks a domain into subdomains recursively using domain splitting lines, which are determined by various geometric rules of subdivision. This algorithm is general purpose when clean up and smoothing aspects are added. The meshes created are generally unstructured. This paper presents a mesh generation algorithm for quadrilateral nite elements that improves upon the basic concept of recursive domain decomposition as proposed by Sarrate and Huerta [11] and is implemented using object-oriented programming. Direct implementation of Sarrate and Huerta’s algorithm frequently results in poor meshes that in many cases require a make-up stage to repair. In many cases humane interference is needed to repair the nal part of the mesh. These drawbacks can basically render the algorithm practically useless for the purpose of automatic mesh renement and adaptive numerical computations. To overcome these di culties, we have developed a general and more e cient remedy by which all the poorly meshed cases are identied and treated each by a special scheme. These critical improvements, combined with the original scheme as proposed, allow the mesh generation algorithm to become truly free from human interference. The other improvement comes from the use of object-oriented programming (see below), which provides substantial computing and programming advantages. Compared with the existing paving algorithm, the present scheme represents a new and diAerent concept for quadrilateral mesh generation. The paving algorithm generates a mesh by adding quadrilaterals one by one, whereas this present scheme generates the mesh through recursively sub-dividing the domains. It should be stressed that automation of mesh generation completely free from human interaction is essential for practical computations and in particular for adaptive nite element applications. Another important aspect of the improved mesh generation algorithm, as described, is that is implemented using the C++ object-oriented programming language. The discussion of C++ implementation has not yet been made, to our best knowledge, on mesh generation for nite element applications. So far, most mesh generation code has been developed using the traditional FORTRAN language, which is a natural choice from the point view of continuity in downstream data processing [6]. Such a natural choice may not necessarily be the optimal choice. In fact, the use of an object-oriented language is more desirable for the improved mesh scheme, because the paradigm of object-oriented programming allows the diAerent parts constituting the mesh to be described easily and naturally as if they were real world objects. The characteristics of an object-oriented programming language, encapsulation, inheritance, and polymorphism, are particularly useful for mesh generation applications [12]. Encapsulation encompasses the idea of classes. A class is a way of programmatically describing objects in the real world. A class is a data structure that has member variables that describe the object to be encapsulated, as well as member functions that allow the object to do work appropriate to that object. Polymorphism is the ability for a program to have multiple instances of the same function. Each function is distinguishable only by the number and=or type of function arguments. Inheritance is the ability for a class to be constructed from an existing class, receiving by inheritance all of its parent class’s member variables and functions. Classes that are appropriate for a mesher include lines, elements, or nodes. All of these objects are easily visible in the real world and have properties, such as connectivity and length, that are easily addressed in a class structure.

908

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

2. Mesh generation algorithm 2.1. Domain splitting criteria The concept of quadrilateral mesh generation is based on the recursive subdivision, which is described in detail by Sarrate and Huerta [11]. The meshing idea is illustrated in Fig. 1. In the rst step, a domain is divided into two equal subdomains by determining a splitting line that will connect two of the nodes in the domain to be split. The decision as to which splitting line is the best is based on a set of geometric rules that judge what subdomain will likely produce the best elements. Each new line also produces its own nodes. Now that there are two subdomains, the splitting rules are reapplied to subdivide the subdomains further. This process is repeated until only elements of the desired type and size remain. The domain splitting algorithm uses the following four basic rules that are derived from geometric considerations: (1) subdomains of equal area are desirable for area splitting, (2) a shorter splitting line is better than a longer one, (3) a good split line splits the previous domain at a ◦ right angle, and (4) a subdomain should not be split at an angle of ¡ 120 . The nal mesh is generated by applying these four rules to determine the best split line that subdivides the computational domain in question. The details of these rules and the associated arguments are given in [11]. These rules may be brie9y summarized in terms of the following mathematical

Initial Computational Domain

Fig. 1. Schematic illustration of recursive domain decomposition to generate quadrilateral elements.

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

909

equations: =

|a2 − a1 | ; a2 + a 1

(1)

‘split ; ‘bb  if (1 + 2 ¡ =2) and (3 + 4 ¡ =2);  1 if (=2 6 1 + 2 6 2=3) or (=2 6 3 + 4 6 2=3);  = ( 1 ; 2 )   (1 ; 2 ; 3 ; 4 ) if (1 + 2 ¿ 2=3) and (3 + 4 ¿ 2=3);

‘=

( 1 ; 2 ) = (1 − 1 ) + 1 2 (1 ; 2 ; 3 ; 4 ); (1 ; 2 ; 3 ; 4 ) =

|1 − 2 | + |3 − 4 | ; 1 +  2 +  3 +  4

  (1 + 2 ) − 2=3 if (=2 6 1 + 2 6 2=3); =6

1 =  1 otherwise;   (3 + 4 ) − 2=3 if (=2 6  +  6 2=3); 3 4 =6

2 =  1 otherwise;  80 0 6 i 6 2=3;

i + j with i =

= 200 0 2=3 6 i 6 2; c1  + c2 + c3  + c4 ‘ + c5 :

(2)

(3)

(4) (5)

(6)

(7)

(8) (9)

In the above, Eq. (1) denes the criterion for area splitting, where a1 and a2 are the areas of the new potential subdomains. Eq. (2) is used to determine the splitting line, where ‘split and ‘bb are the lengths of the splitting line and the bounding box, as shown in Fig. 2a. Eqs. (3) – (7) are for perpendicular bisector determination, where ’s are the angles the candidate splitting line makes with the segments of the domain it will split, as shown in Fig. 2b. It should be noted that the angles are measured from the base of the polygon to the split domain and 0 6 i 6 2. Eq. (8) measures the uniformity of local structuring, where  is the measurement between adjacent line segments before splitting, as shown in Fig. 2c, and it is most easily calculated as i = 1 + 2 and j = 3 + 4 (see Fig. 2b). In an ideal quadrilateral mesh, four adjacent nodes meet a central node at right angles from each other, as seen in Fig. 4. This implies that the line segments being split by the candidate ◦ split line ideally should be 180 apart. This leads to a penalty being applied if the splitting line tries to split segments that are already at small angles from each other. The last equation, i.e. Eq. (9), is the mesh generation function, which is to be minimized to determine the best splitting line. The constants are determined as follows: c1 = 0:52; c2 = 0:17; c3 = 0:0; c4 = 0:17 and c5 = 0:14 [11].

910

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930 Candidate Split Line

α3 α1

α4

α2

(a)

Domain Bounding Box Diagonal

(b)

Opposite nodes 180° apart θ

Adjacent nodes 90° apart (c)

Fig. 2. Geometric criteria used for generating quadrilateral elements: (a) segments used for line length criteria, (b) angles used for perpendicularity criteria and (c) angles used for structuring criteria.

2.2. Improved algorithm In general, the algorithm described repeats until a 4-node element results and produces a reasonable mesh for regular geometries such as squares or rectangles, when an appropriate smoothing scheme is used. For many irregular geometric domains, however, an unsatisfactory mesh outlay is generated with distorted elements and 6-node subdomains before reaching the 4-node stop criteria. Sarrate and Huerta [11] suggest using a diAerent set of weights for the basis function in this case, setting c4 and c5 to 0 and having non-zero c3 . Unfortunately, these remedies do not always work out well and, in fact, in many cases tested, the algorithm terminates with a 6-node element that cannot be further reduced into 4-node elements, such as the one shown in Fig. 3. This conguration could be split along the dashed line, but the resulting element on the right is illegal, having a node that is collinear with its neighbors. Through extensive testing, we found that the most eAective means to resolve these 6-node cases is to treat them separately, rather than to develop a generalized remedy to be included in the general formulation. Our exhaustive test of conceivable 2-D geometries shows that all 6-node domains can be classied into one of eight categories, as shown in Fig. 4 [13]. The nomenclature for classication consists of a sequence of integers. The number of integers represents the number of sides of the domain.

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

Fig. 3. Typical poorly dened subdomain without improvements.

!

!

4-1-1 3-1-1-1

!

!

3-2-1 or 2-3-1

2-2-2

!

!

2-2-1-1

2-1-1-1-1

2-1-2-1

1-1-1-1-1-1

Fig. 4. Classication of misshaped 6-node subdomains through exhausted testing.

911

912

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

The integer values represent the base 0 count of collinear nodes forming a side. Note that the order of node sequence matters, as a 2-1-2-1 and a 2-2-1-1 are handled diAerently. Resolution of each of these cases is now described below. The 4-1-1 case is a rarity that only occurs when a high nodal gradient is specied during domain formation. In other cases, the program’s node spacing requirements will add nodes to the segments connecting the single node to the segment with ve nodes, so this case does not otherwise occur. Since this type of element has a known cause and results in very poor elements, it is ignored in the program. If this element type does occur, the user should reform the domain with a smoother node spacing so that the situation does not repeat. If explicit treatment is desired, the easiest solution would be to add a line segment consisting of two nodes from the peak of the triangle to the midnode of the collinear nodes. This segment should have two additional points, thus subdividing the triangle into two 3-2-1 domains, which case is handled below. It should be noted that the other cases fall into three general categories—those in which no additional nodes are needed, those where one additional node is needed, and those where two additional nodes are needed. Cases with four collinear nodes are the worst cases; they require that two nodes be added to form acceptable elements. Cases with multiple sides of three adjacent collinear nodes are the second to worst case, requiring a single node be added. The remaining cases are trivial, requiring no node addition. The 3-1-1 case may be resolved by vectorizing sides 1 and 2 of the polygon, as shown in Fig. 5a. The vectors are then reduced in magnitude by an arbitrary amount, generally about 1=3 their original length. These vectors are then added to the two base points of the side with four nodes, as shown. Note that there is a potential problem with the vectors crossing if they are very long. This problem is generally resolved by making the vectors su ciently short. This is followed by mesh smoothing by overrelaxation to completely eliminate the problem. In the case of a 3-2-1 (or 2-3-1) domain, the problem is nearly identical to that of the 3-1-1-1 domain. Since there are fewer sides, the vectors formed by sides 1 and 2, shown in Fig. 5b, can lead to more distorted elements and the crossover problem described above. In practice, the vectors generated are scaled even further to solve this problem. Additionally, like the 3-2-1 case, this problem is easily resolved when smoothing is performed. Let us now look at cases where only a single central element is added, the 2-2-2 and 2-2-1-1 cases. First, the 2-2-2 case is resolved easily by adding a node at the centroid of the triangle formed by the nodes. The new center node then connects to the triangle midsides, forming three quads. This process is shown in Fig. 5c. The 2-2-1-1 case, shown in Fig. 5d, is nearly identical to the 2-2-2 case. A node is added at the centroid. Elements are then constructed by joining the center node with every other node in the 6-node loop beginning with a center node of a collinear trio. Other cases are relatively simple and the problems can be xed as illustrated in Figs. 5e–g. Experience has shown that, of all these cases described here, the most common case is the 2-1-2-1, followed by the 2-2-1-1 and 3-1-1-1. 2.3. Mesh smoothing The smoothing algorithm described here is a Laplacian technique, named because it is solved by methods identical to the numerical solution of the Laplace equation, with nodal position as a Dirichlet boundary condition. The Laplace equation describes elliptic partial diAerential

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

913

equations for which the solution is felt at all points of the domain. In the mesh that has been constructed, if a node is moved, all elements sharing it are aAected, and thus all surrounding elements are aAected. This is analogous to incompressible 9uid 9ow or heat transfer problems at steady state.

1 2

(a)

2!

2

1

(b)

!

(c)

!

(d)

!

(e)

Fig. 5. Treatment of misshaped subdomains case by case: (a) 3-1-1-1 subdomain resolution, (b) 3-2-1 or 2-3-1 resolution, (c) 2-2-2 resolution, (d) 2-2-1-1 resolution, (e) 2-1-2-1 resolution, (f) 2-1-1-1-1 resolution and (g) 1-1-1-1-1-1 resolution.

914

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

(f)

(g)

Fig. 5. Continued.

As in 9uid mechanics, a simple way to solve this problem is to smooth the node coordinates using the successive overrelaxation method or SOR. In this method, each node is individually adjusted by setting its new position to the average position of those nodes linked to it [14]. This is performed repeatedly until no more change in the nodal positions is evident. Generally, 10 iterations are su cient. To perform SOR, it is necessary that each node knows its neighbors. Additionally, when a node position changes, that change must be re9ected in all elements that share the node. These requirements are easily fullled using the data structures inherent in the C++ implementation of this algorithm as described below. 3. C++ implementation The mesh generation algorithm described above can be best implemented using the C++ objectiveoriented programming language. This is because object-oriented programming languages take the advantage of the data structure inherent in mesh generation processes. The implementation involves the use of the basic object of line segments, a loop of which forms a domain. The repeated use of these line and domain objects accomplishes the recursive decomposition of the geometry into quadrilateral elements. While there have been many mesh generators available, FORTRAN was generally used for the convenience of downstream data processing. In our approach, C++ was used to implement the mesh generation algorithm, as the embedded data structure greatly increases the computational e ciency of mesh generation. Since the full details are documented in a recent thesis [13], only the essential features of the C++ implementation used in our program are summarized below. 3.1. Line segment object Implementation using C++ entails creating various objects. One crucial object for meshing is a line segment, the data structure of which has the following properties: • The line segment can point to the next and previous segments in the loop. • The segment stores its endpoints.

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

915

• The segment can construct its own nodes, based on knowing the desired node density and whether the node count in the potential loop to be constructed is even or odd. • The line can break itself into two segments, keeping one portion and returning the remainder, so that it can be used in a new loop. Note that this function also divides up the nodes properly, so that the right nodes go to the right line. • The line also stores whether or not the segment is an edge segment or an interior segment. 3.2. Element object The element list consists of element objects with the following signicant member variables and functions: • • • • • • •

Global integer address of each element. Global integer addresses of each node in the element. Floating point storage for node coordinates. Local integer addresses of each node in element. Memory to save addresses of neighbor elements. Integer 9ag saving which sides of element, if any, are on a boundary. Function that, given the address of another element, determines if the other element is a neighbor and what side of the element the neighbor borders on. This function is reciprocal—when an element nds a neighbor and updates its own information, the neighbor’s information is updated, too. • A pointer to another element. This is used to construct linked lists of elements, as is used in the program. 3.3. Domain construction The mesh generation starts with the construction of a domain, which is done by drawing the domain counterclockwise until the end and beginning segments are reconnected. Holes or other geometric objects can be fed into the interior of the domain. This process is shown in Fig. 6. It is not important where the start point of the loop is, so long as a counterclockwise path is always traced. As each segment is drawn, it is added to a linked list of segments for storage and manipulation. Upon creation, the segment also creates its nodes. All boundary segments have an even number of nodes to ensure that the domain starts with an even number of nodes. For any one domain, such as shown in Fig. 7a, a corresponding doubly linked list is created (see Fig. 7b). This allows for easy traversal in either direction along the list. With the list constructed, the line splitting subroutine, GetSplitLine, is called (see Fig. 7d). A pointer to the rst element in the linked list of segments, pElement is passed to the subroutine. Since the subroutine is recursive, depth control is required. In the program, the recursion depth is controlled by checking that the variable “steps” is greater than zero. Also, the program needs to make sure that an empty list was not passed to the routine. This will generate a list joined front to back to form a loop out of the initial list in Fig. 7b. The segments of this loop are counted. Memory is allocated to store the loop vertices and the vertices are stored into arrays x and y.

916

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

Fig. 6. Clockwise construction of a domain with a hole inside.

The program then loops through each element of the list, creating and evaluating split lines. If the loop shown in Fig. 7b were the loop containing the domain data, the algorithm would compare all possible lines formed between segment 1’s nodes and the nodes of all other segments. Then, all possible dividing lines formed by the nodes of segment 2 and the segments after it would be compared, and so on, until all combinations of segments and nodes have been evaluated. This process is depicted graphically in Fig. 7c, where pTop and pNext are pointers to list segment. Once there are no pNexts left, pTop advances itself until all possible segment combinations have been exhausted. 3.4. Determination of the best split line With line segments and domains constructed above, best split lines are determined. The basic procedure is shown in Fig. 7e. Each of the criteria listed in Eqs. (1) – (8) is calculated and checked individually and the line that minimizes the base function (Eq. (9)) is chosen to be the split line. In determining the best split line, additional checks on meshing criteria must also be made for each line, which are described below. 3.4.1. Rejected cases The validity of the split line is ensured by performing the following tests: • The splitting line may not intersect any segment of the domain, as shown in Fig. 8. • The splitting line may not pass outside of the domain to be split, as shown in Fig. 8. • The splitting line must have at least three nodes to form a domain greater than zero. If any combination of segments or points fail these tests, the line is immediately rejected and the next potential splitting line is evaluated.

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

917

3.4.2. Intersection test To ensure that a potential segment does not pass through another segment, a test is performed using parametric line segments, which are dened below, p = p0 + u(p1 − p0 )

and

l = l0 + v(l1 − l0 ):

(10)

Here, p and l are points on line segments running from p0 to p1 and l0 to l1 , respectively. The parameters u and v range from 0 to 1. The point of intersection where p = l is determined by

5

4 3

6 CCW

2 1 (a)

! NULL

1

2

3

4

5

NULL

6

(b)

pNext pNext pNext pNext

pNext

pNext pNext pNext

6

pNext

Iteration 5

5

pTop

pNext pNext

Iteration 4

pTop

pNext

pNext

Iteration 3

pTop

Iteration 2

4

pNext

3

pTop

pTop

Iteration 1

2

pNext

1

(c)

Fig. 7. Representation of a typical subdomain: (a) clockwise traversal of the domain, (b) doubly linked list representation, (c) order of segment loop traversal, (d) initializations for line splitting algorithm and (e) line splitting algorithm.

918

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930 Subroutine GetSplitLine(pElement, Steps)

pTop = pElement

Steps < 0

True

End

False pTop = NULL

True

End

False Link Segments into Loop Count Loop Segments

Store Segment Vertices in x and y arrays

Go to step 2 (d)

Fig. 7. Continued.

calculating u and v using the following equations:       (px − px ) −(lx − lx )   u   lx − px   0  1 0 1 0   0  ;   =   y  (p1 − p0y ) −(ly1 − ly0 )   v   ly0 − p0y 

(11)

where superscripts denote the x and y components of the vectors. If either u or v is between 0 and 1, inclusive, the splitting line intersects another segment and is rejected. 3.4.3. Exterior segment and zero area rejection To detect lines that fall completely outside the domain to be split, a test similar to a computer graphics scan line conversion algorithm for polygon lling is used [15]. A test point is created by using the parametric equations above with u = 0:5. A horizontal test line is created that passes through the test point. The test line is then checked for intersections with each of the line segments in the polygon. If the number of intersections to either the left or the right of the test point is odd, the point lies in the polygon. This procedure is illustrated in Fig. 9. If the split line passes the above two tests, it is considered valid and is kept for further comparison.

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

2

PTop->next = pElement

False

Line1_pts = pTop # nodes. Fill x1 and y1 arrays with pTop node coordinates.

True pNext = pTop->next pTop = pTop->next True

pNext->next = pElement False

Line2_pts = pNext # nodes. Fill x2 and y2 arrays with pNext node coordinates.

False

pNext = pNext->next

i < Line1_pts

i=0 i++

True False

j < Line2_pts

j=0 j++

Valid split line found?

True

False

p0[0] = x1[i]; p0[1] = y1[i]; p1[0] = x2[j]; p1[1] = x2[j];

True End

Segment is best split line

Call Splitdomain

False

True Store segment data and score (e)

Fig. 7. Continued.

919

920

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

5

4 6

3 2 1

Fig. 8. Invalid splitting lines. These lines intersect other lines in the domain or fall outside the domain.

2 1

4 3

Fig. 9. Test of interior nodes using scan line conversion concept.

3.5. Domain subdivision Once the dividing line has been determined, the domain must be split into two subdomains. This is done by subroutine SplitDomain, which sends the addresses of the segments containing the nodes that will form the new dividing line and the points that form the endpoints of the dividing line. This subdomain then breaks the loop into two new lists of segments, which form new subdomains. This process continues until one of the 6-node cases discussed in Section 2.2 is encountered and then the special treatment is used to further break these 6-node elements into quadrilateral elements. Once all the elements are generated, the elements pass through the mesh smoothing subroutine, which implements the Laplacian algorithm with successive overrelaxation. 4. Results and discussion The mesh generation algorithm as described above has been applied to generate quadrilateral elements over a wide range of geometric congurations. A selection of these results is given below.

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

921

The mesh generation computations consist of several steps: (1) entry of data as line segments, (2) construction of domains based on line segments, (3) recursive subdivision of domains based on previous domains, treating domains with six nodes as special cases, (4) storing elements as they are constructed, (5) relaxing the mesh for better elements, and (6) storing the element and node data for further processing. Computations for mesh generation are e cient. For 534 number of elements with 599 nodes, for example, it took about 1:87 CPU s on a Pentium II PC (450 MHz) and 1.04 on a Pentium III PC (800 MHz).

4.1. Meshed cases 4.1.1. Rectangular domain For the rst case, the simple rectangular domain of Fig. 10a is meshed. We use this case also to illustrate the general procedure taken in the described algorithm for mesh generation. Note that the initial domain shows all initial nodes as small squares on the domain boundary. Running the algorithm gives the result shown in Fig. 10b, which is the same as any structured mesh generator such as transnite mapping would produce as well. As should be expected for this simple square domain, a very regular mesh is obtained as intended here, consisting only of horizontal and vertical lines. This is a consequence of satisfying the angle criterion for line splitting. Mesh smoothing produces the same result as that shown above, since the mesh already consists of perfect quadrilaterals. Mesh renement with desired mesh distribution is obtained by simply breaking a domain’s set of line segments into multiple segments of diAerent node densities. An example of mesh with node density control is shown in Fig. 11a, where a pre-set node density is placed at the top right corner of the square. The resulting mesh without smoothing is shown in Fig. 11b. Clearly, it suAers from very distorted elements, though they are all quadrilateral. These distorted elements are improved signicantly by mesh smoothing, as is evident in Fig. 11c. 4.1.2. Polygonal bracket domain For the second example, a domain consisting of a bracket-type polygon is meshed. In a bracket of this type, the stress concentrations would be located at the interior corners. It is desirable to have more elements in this location. This is accomplished by increasing the node density in the line segments containing those corners. Additionally, the element structure and size can be further controlled by creating leaders extending them into the domain. The result of this meshing scheme after smoothing is given in Fig. 12. Note that the node density is very high in the regions where high density is desired. 4.1.3. Rectangular domain with a rectangular hole In this case, a rectangular domain with a rectangular hole is created. Note the bi-directional leader used to construct the hole. To improve the above mesh, an area of high element density is desired around the inner hole. This is achieved by specifying the node density on the interior segments as well as adding leaders into the domain. The resulting mesh after smoothing, is shown in Fig. 13.

922

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

(a)

(b)

Fig. 10. Mesh generation in a square domain: (a) distribution of boundary nodes and (b) meshed square domain.

4.1.4. Other complex cases The similar procedures are applied to generate meshes in more complex geometries of either simply connected and multiply connected domains, as illustrated in Figs. 14 –18. Fig. 14 depicts a rectangular domain with a diamond-shaped hole located at the center. The leaders are added to further control node spacing. The nal mesh shown is improved by smoothing. Fig. 15 shows the mesh generated in a domain with a large hole inside. The algorithm should also produce an axially symmetric mesh, since the domain has this symmetry. As expected, the

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

923

(b)

(a)

(c)

Fig. 11. Quadrilateral mesh generation in a square domain with varying node distribution: (a) square domain with varied boundary node density, (b) mesh before smoothing and (c) mesh after smoothing.

mesh exhibits the same symmetry as the initial domain. The results without smoothing produce very distorted elements, which are improved signicantly by the smoothing algorithm. Meshes shown in Fig. 16 are created for a circle with an inner diamond-shaped hole. The hole has a greater node density along its segments as well as leaders to promote good element generation. Evidently, the resulting mesh with smoothing produces excellent quadrilaterals.

924

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

Fig. 12. Meshed polygonal domain with varying element density after smoothing.

Fig. 13. Meshed rectangular domain with a rectangular hole after smoothing.

Fig. 17 considers the case where a small hole has been added oA center from a large circle. The desired result is a heavy concentration of nodes near the hole. To facilitate this, high node density has been specied near the hole and leaders are inserted which gradually decrease the node density as their proximity to the hole increases. The nal result with smoothing is a mesh with spacing as desired. As the last example, a domain with two holes of varying shapes is considered. The domain is multiply connected. The meshes with a given distribution are generated and the results with smoothing are shown in Fig. 18. For this case, as for other complex cases described above, meshing without improved algorithm would require human inspection and interference to correct many misshaped and sometimes non-quadrilateral elements. While the improved algorithm ensures that the non-quadrilateral elements do not occur, smoothing is critical and must be applied to straighten out

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

Fig. 14. Smoothed mesh distribution in a square domain with a diamond-shaped hole at the center.

Fig. 15. Distribution of quadrilateral elements in a square domain with a circular hole.

925

926

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

Fig. 16. Smoothed mesh distribution in a circular domain with a diamond-shaped hole.

Fig. 17. Smoothed mesh in a circular domain with an eccentric round hole.

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

927

Fig. 18. Smoothed mesh distribution in a domain with multiple holes. Table 1 Comparison of algorithm execution speeds Shape

Element count

Node count

Mesh time (s)

Smooth time (s)

Square Rened square Polygon Rened polygon Rectangle hole Rened rectangle hole Diamond hole Round hole Circle Circle with diamond hole Circle hole in circle Two holes

168 621 384 1124 420 641 2310 560 640 684 3674 8079

195 679 443 1204 486 716 2392 650 673 732 3818 8387

0.38 1.98 1.04 4.67 1.38 2.86 12.41 2.52 1.81 3.18 33.78 96.39

0.38 4.23 1.65 13.51 1.92 4.34 60.36 3.24 4.61 5.16 183.12 1333.04

the misshaped elements. For this case, in particular, the smoothing algorithm signicantly enhances the quality of the mesh. 4.2. Algorithm speed To evaluate the speed of this algorithm, many diAerent geometric domains were meshed and rened (or smoothed) on Pentium II 450 MHz machine. Table 1 compares the computational e ciency and CPU times for meshing and smoothing for various cases on the Pentium PC machine. A few

928

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

samples produce very regular meshes, similar to that structured mesh generators produce, and as such no smoothing is needed. Smoothing has been performed in all cases, however, so that a quantitative comparison can be made. Also, the smoothing value shown is obtained after ten iterations of smoothing. In general, only three or four iterations are su cient. It is noted that smoothing time increases drastically with an increase in node numbers. For the last case in the table, the smoothing time was about 14 times that required for meshing. In our program implementation, the basic successive overrelaxation scheme was used. Apparently, an optimized successive overrelaxation scheme should help to speed up the computation for mesh generation. 4.3. Discussion Here, selected examples illustrating the capabilities of the algorithm have been presented, along with computational e ciency. The program is very 9exible, successfully meshing a variety of shapes, including those with holes. While for most tested cases, the algorithm took just over 3 s to smooth a domain, the CPU time required for smoothing increases drastically with the number of nodes and elements (see Table 1). The algorithm presented, including the smoothing stage, produces domains that are generally well structured with well-shaped quadrilateral elements. As with any other meshing algorithms, the degree of structure to which the elements obtain is heavily dependent upon the domain being meshed as well as node distributions. Domains with sharp corners will most certainly have misshapen elements. However, this problem can be alleviated with smaller elements, followed by mesh smoothing, as shown in the above gures. In this regard, the improvements prove to be extremely successful and completely eliminate misshaped elements with any human interference. This algorithm has been implemented using object-oriented methods. The objects used are line segments and elements. The use of these data structures allows the programmer great 9exibility in algorithm implementation. Many actions, such as line splitting and element storage, would certainly be more di cult using a non-object-oriented method. For example, in this implementation, a domain is represented by a doubly linked list forming a loop. This is an ideal way to think of a polygon, especially when the polygon is to be traversed. Additionally, the ability for a segment to split itself and preserve both segments is very useful in domain subdivision. Element creation also benets greatly from the object-oriented paradigm. The elements have the ability to compare themselves to one another and preserve relationships when they are determined. Additionally, the element’s ability to save all of its pertinent data is of great use for writing les for use in nite element programs. Not only is programming simplied when objects are used, but speed is increased. When comparing two nodes for best splitting lines, line segments do not compare nodes on the same segment. This immediately eliminates all nodes on the segment the starting point of the potential splitting line is on. By automatically doing this, a great saving in computation time is achieved. While the algorithm is designed and has been implemented for 2-D geometries, its extension to 3-D mesh generation is possible. One of the simplest extensions would be 3-D mesh generation by extrusion. There are many applications in which 3-D computational domains are simply the geometric congurations the parallel planes of which along the third dimension are either the same or similar in structure. For example, a simple cylinder with and=or without a through hole of varying size along the direction of the vertical axis can be readily meshed by rst generating a 2-D plane structure using

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

929

the algorithm presented and extruded along the vertical axis. For more general cases, however, the basic mesh criteria, along with appropriate treatment schemes for poorly meshed elements, must be developed. This development is envisioned to involve the denition of geometric relations of splitting planes in addition to those for lines. The C++ implementation of these extensions, however, should be rather straightforward and requires the addition of additional components in the data structure for 3-D treatment. 5. Concluding remarks This paper has presented an e cient algorithm for creating unstructured quadrilateral meshes for nite element or unstructured numerical computations. The algorithm was developed by improving upon an existing recursive mesh generation reported in literature and implemented using the C++ object-oriented programming language. The improvements were based on the treatment of exhausted cases of misshaped elements for which a general remedy is not possible or does not exist. These improvements have led to an increased computational e ciency and to a complete automation of meshing operation without a need for a nal mesh resolution by human interference. Tests show that for all the cases studied, except the regular, structured mesh, these improvements are required to make the mesh generation automatic. The implementation of the algorithm using object-oriented methods proves to be very benecial both in terms of programming eAort and computational e ciency. The objects used are line segments and elements. The use of these data structures allows the programmer great 9exibility in algorithm implementation. Numerical experience with meshing over a variety of geometries shows that not only is programming simplied when objects are used, but also computational speed is increased. When comparing two nodes for best splitting lines, line segments do not compare nodes on the same segment. This immediately eliminates all nodes on the segment the starting point of the potential splitting line is on. By automatically doing this, a great saving in computation time is achieved. Acknowledgements The authors gratefully acknowledge the support of this work by NASA (Grant Nos.: NAG-1477 and NAG8-1693). References [1] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, 4th Edition, McGraw-Hill, New York, 1988. [2] F. Bassi, S. Rebay, A high-order accurate discontinuous nite element method for the numerical solution of the compressible Navier–Stokes equations, J. Comp. Phy. 131 (1997) 267–279. [3] S.P. Song, B.Q. Li, Coupled boundary=nite element analysis of magnetic levitation processes: free surface deformation and thermal phenomena, J. Heat Transfer 120 (1998) 492–504. [4] S.P. Song, B.Q. Li, Coupled boundary=nite element solution of magnetothermal problems, Int. J. Num. Methods Eng. 44 (1999) 1055–1077. [5] S.J. Owen, A survey of unstructured mesh generation technology, Proceedings, Seventh International Meshing Roundtable, Sandia National Laboratories, Albuquerque, NM, 1998, pp. 239 –267.

930

M. Bastian, B.Q. Li / Finite Elements in Analysis and Design 39 (2003) 905 – 930

[6] T.D. Blacker, M.B. Stephenson, Paving: a new approach to automated quadrilateral mesh generation Int. J. Num. Methods Eng. 32 (1991) 811–847. [7] D.R. White, P. Kinney, Redesign of the paving algorithm: robustness enhancements through element by element meshing, Proceedings, Sixth International Meshing Roundtable, Sandia National Laboratories, Albuquerque, NM, 1997, pp. 323–335. [8] J. Heinrich, D. Pepper, Intermediate Finite Element Method: Fluid Flow and Heat Transfer Applications, Taylor & Francis, Castleton, 1999. [9] FIDAP Manual: Preprocessing, 1996. [10] G. Farin, Curves and Surfaces for Computer-Aided Geometric Design: a Practical Guide, 4th Edition, Academic Press, San Diego, CA, 1997. [11] J. Sarrate, A. Huerta, E cient unstructured quadrilateral mesh generation, Int. J. Num. Methods Eng. 49 (2000) 1327–1350. [12] I. Horton, Beginning Visual C++ 6.0, Wrox Press, Canada, 1998. [13] M. Bastian, MS Thesis, Washington State University, Pullman, WA, 2001. [14] S.C. Chapra, R.P. Canale, Numerical Methods for Engineers, 3rd Edition, McGraw-Hill, New York, NY, 1998. [15] J.D. Foley, Introduction to Computer Graphics, Addison-Wesley Publishing Company, Reading, MA, 1997.