An Automatic Method for Complete Triangular Mesh Conversion into Quadrilateral Mesh for Multiple Domain Geometry

This research developed an automatic two-dimensional finite element meshing system to resolve practical engineering problems in the fields of geology, hydrology, and water resources. This system first used the Delaunay triangulation method to create reasonable-density triangular mesh and then converted it into quadrilateral mesh by combining proper pairs of adjacent triangles. A series of combination patterns aiming at three cases were established. The effect of the number of boundary edges on the subsequent meshing procedures were studied and summarized. For the geometry with multiple domains an adjustment method is proposed to completely eliminate the residual triangles during quadrilateral meshing through adjusting the number of boundary edges in each loop to be even. A special boundary loop identification method is proposed for priority treatment. Corresponding treatment methods aimed at three different situations are established for common boundary loops. For a certain boundary loop with an odd number of boundary edges, the appropriate edge for new point insertion is determined by the position properties and relative density errors. Practical applications confirm that the method proposed in this paper could successfully implement the full conversion from the triangular mesh to the quadrilateral mesh.


INTRODUCTION
To study complex problems in practical applications of water flow, thermal transport, salinity, sediment, and pollutant transport, numerical simulations have been employed for several decades. The finite element method was introduced into engineering applications and gained popularity and effectiveness in dealing with various types of troublesome issues. The finite element method can approximate irregular geometry more accurately because of its flexibility in choosing computational grid sizes.
The whole finite element simulation and numerical analysis process includes three major stages: pre-process, computation, and post-process. Finite element mesh generation is the essential aspect of preprocessing, which is still sometimes a bottleneck in the overall simulations.
Because of the explicit definition requirements for each nodal point and element connectivity, the finite element meshing output preparation becomes very tedious as the computational geometry of a physical problem becomes very large and complicated. In practical applications, such as computational hydraulics, the numerical results are converged in order to insure the computational grids should be refined globally or locally. An automatic mesh generator is required to ease the burden in the finite element computations. This generator for finite element simulations should generate node coordinates and the connectivity of the elements with given specifications. The interaction and mapping relationships between one-dimensional river/ stream, junctions with or without storage and two-dimensional overland flows (or surface runoff) are required as output data according to the users' requirements. In view of these demands, key finite element meshing techniques are researched in detail in this paper to develop an effective generator to solve practical engineering applications in geology, hydraulics, hydrology, water resources, fluid mechanics, and mechanics, etc. The finite element mesh is divided into two categories, two-dimensional mesh, and three-dimensional mesh, according to the solution space. The two-dimensional mesh further contains triangular mesh, quadrilateral mesh and hybrid mesh (Ito 2013).
Among the numerous kinds of triangular mesh generation methods, the Delaunay triangulation method is the most widely used method. It has two crucial criteria: maximum-minimum angle criterion and in-circle criterion. Sibson (1978) proved that the Delaunay triangulation result for an arbitrary set of given points in a common plane possessed the optimal advantage. The sum of minimal internal angles of all Delaunay triangles obtains the maximal value and the maximal internal angles will produce the minimal value. The maximum-minimum angle criterion is able to avoid producing long and thin elements with small internal angles automatically. The in-circle criterion signifies that the circum-circle of each Delaunay triangle will not contain any other points in the given set. Bowyer-Watson algorithm (Bowyer 1981;Watson 1981) just takes advantage of the in-circle criterion. In recent years, Delaunay triangulation algorithm has attracted many research scholars because of its excellent flexibility and extensibility (Žalik 2005;Ebeida et al. 2011). Therefore, this paper adopts Delaunay triangulation method to generate triangular meshes, which are used as the background grid to create quadrilateral meshes.
The most extensively applied quadrilateral meshing methods are grouped into two main categories: direct method and indirect method (Ho-Le 1988). The direct method generates quadrilateral meshes directly based on the geometric features of the analyzed geometry. Direct methods include the medial-axis method (Tam and Armstrong 1991), quad tree method (Yerry and Shephard 1983), and advancing front method (Blacker and Stephenson 1991), etc. The indirect method is also known as the triangle-to-quadrilateral conversion method. That is, triangular mesh is generated first and then converted into quadrilateral mesh (Shimada and Itoh 1995;Borouchaki et al. 1996).
The indirect method takes advantage of the benefits of triangle meshing, thus has more advantages than the direct method in controlling mesh density and element size, meshing efficiency and generation versatility, etc (Itoh and Shimada 2002). There are two implementation styles for the indirect method: the first is the division method; the second is the combination method. The division method subdivides a single triangle into three quadrilaterals by adding three new points at the midpoints of the three edges and one point in the triangle centroid. This method is simple, easy to implement and very versatile. Because a large number of irregular nodes are not connected to four quadrilaterals, the element quality cannot be guaranteed. The addition of the new points and elements has a serious effect on the mesh density dis-tribution, possibly leading to a generated mesh that cannot satisfy the meshing density requirements specified by the users. In addition, the insertion of the new points on curved boundaries might cause curvature errors and area loss.
The combination method generates quadrilateral meshes through combining each appropriate pair of adjacent triangles into a single quadrilateral. This method is relatively difficult in algorithms because the choice of triangle pairs has a great impact on the conversion effect and mesh quality. However, the combination method is able to produce reasonable quality meshes that can accurately capture the geometric features of the input geometry. Since there is no new point to be inserted, the density control of users for mesh elements and nodes can be fully satisfied. As a result this paper uses the combination method to produce quadrilateral meshes.
The implementation procedures for quadrilateral mesh generation using the combination method are often troublesome and time-consuming for complex geometries, especially multiple-domain geometries. The complete conversion from triangles into quadrilaterals for each sub-domain has become a challenging research issue. The conformity of mesh elements and nodes between two adjacent sub-domains must be perfectly ensured. Until now, most literatures focused mainly on research on conversion strategies for triangles and resulting mesh quality optimization. Few papers have been published that pay attention to the number control of mesh elements and nodes for complicated geometries with multiple sub-domains.
This research conducted an intense examination of the key techniques and treatment strategies in the triangle and quadrilateral meshing procedures, containing Voronoi diagram construction, boundary point generation, the insertion of new points, the conversion from triangles into quadrilaterals and the treatment of residual triangles. The odd and even parity properties and changing rules between the numbers of boundary edges and mesh elements in each subdomain were studied and summarized through numerous practical examples. An adjustment method is proposed for controlling the number of mesh elements and triangle conversion pairs by controlling the number of boundary edges in each loop. Corresponding treatment methods for adjusting the number of boundary edges for special and common loops are established. The relative density error concept is proposed to determine the appropriate boundary edge where the new point is inserted. The reliability and effectiveness of the methods presented in this paper are demonstrated using several practical examples.

KEY TECHNIQUES IN TWO-DIMENSIONAL MESH GENERATION
The two-dimensional finite element meshing system developed in this paper adopts the Delaunay triangulation method to generate density-controlled triangular meshes.
Based on the Delaunay meshes an improved conversion method is applied to further generate quadrilateral meshes by combining appropriate pairs of adjacent triangles.

Theoretical Foundation
The Delaunay triangle meshing program was developed based on the Voronoi diagram. The curved input geometry boundaries are constructed using B-splines according to the point spacing specified by the users.

The Voronoi Diagram
Suppose that a set of points S in a two-dimensional plane are given, with n distinct points in all. For each point P i , its Voronoi territory, labeled as V(P i ), is defined as the locus of points in the common plane whose distance is closer to P i than to any other points of S. V(P i ) is expressed as (Bowyer 1981): The Voronoi cell of each point in set S can be expressed as a convex polygon using Eq. (1). As a result, the whole plane where S is located is divided into n Voronoi polygons, each of which is associated with a unique point. The n convex polygons form the Voronoi diagram of the point set S. The straight line duality of the Voronoi diagram creates Delaunay triangles. Each Voronoi vertex corresponds to a unique Delaunay triangle and vice versa.

B-Spline
B-spline is short for basis spline, first proposed by Schoenberg (1946). It possesses strong adaptability and flexibility, and is widely used in the fields of computer-aided design and manufacturing.
Assuming that P(u) represents the position vector of an arbitrary point on a curved boundary, the k-th order B-spline along variable u is defined by: where B i indicates the position vectors of the control points. N i, k (u) denotes the normalized B-spline blending functions. While constructing a k-th order B-spline, the Cox-de Boor formula (Cox 1972;De Boor 1972) of the i-th blending function is defined as: where u i represents the knots of the B-splines, with u i ≤ u i + 1 . Note that the entire value of any term in Eq. (4) is assigned zero when the denominator is zero.

Input Data
The meshing system developed in this paper requires two input data aspects, the formed points, and the boundary edges of the geometry.

The Formed Points
Assume there are np points in total in the set S to describe the boundary features of the input geometry. First of all the position coordinates in the x-and y-axes of the Cartesian coordinate system for each point are required to be given, as P i (x i , y i ) ! S, i = 1, 2, 3, …, np.
In addition, users usually put forward different requirements for the density distribution of mesh elements. Some geometric features that need particular analysis and simulation require a large density to ensure computation accuracy, such as the large-curvature boundaries, slight-size features, and the domains of particular interest to users. The geometric features that do not need focus require a relatively small density to save computation time and storage space, such as the boundaries with a smooth curvature or the domains of no interest to users. The density of mesh elements and points should therefore be specified previously by users. To facilitate computation the mesh density is controlled by the point spacing. The point spacing in this paper is defined as the lengths of the edges connected to the point.

Boundary Edges
The data information stored for each boundary edge is composed of the serial numbers of its two endpoints as well as its curved or straight status. Suppose that there are nb boundary edges, then each edge is denoted as b i (e 0 , e 1 ), i = 1, 2, 3, …, nb, where e 0 and e 1 are the two endpoints.
It should be noted that in order to ensure successful Delaunay triangulation performance the input boundary edges for each meshed region must follow such a rule: the outer edges are arranged in counterclockwise order and the inner edges are arranged in clockwise order. Therefore, the boundary edges on the intersecting lines between two connected geometric sub-domains are overlapped. If a boundary edge in a sub-domain is b i (e 0 , e 1 ), then its corresponding overlapped edge in the neighboring sub-domain is b j (e 1 , e 0 ), j ≠ i. Figure 1 shows a geometric model that requires meshing, containing two domains with different types of materials, the outer domain A and the inner domain B. These two domains should be meshed separately and compatibly. P 1 , P 2 , …, P 12 represent the data points of the geometry, with np = 12. These 12 data points create 16 boundary edges, nb = 16. Since the input geometry includes two connected sub-domains, four pairs of overlapped edges are present on the intersecting boundaries, as P 9 -P 10 , P 10 -P 11 , P 11 -P 12 , P 12 -P 9 in sub-domain A and P 10 -P 9 , P 11 -P 10 , P 12 -P 11 , P 9 -P 12 in sub-domain B. The decimal in the parentheses after the label for each point indicates its required point spacing. If the specified spacing of a point is less than the actual length of the boundary edge connected to it, new points should be generated on this edge.

The Identification of the Geometric Features
In order to achieve effective finite element mesh generation, identifying the boundary loops, characteristic points and boundary segments of the geometric model must be accomplished first.
The boundary loop is defined as the closed ring made up with end-to-end boundary edges. Each geometric model subdomain is bounded by one or more boundary loops. For example, the geometric model as shown in Fig. 1 has two boundary loops, the outer loop L 1 and the inner loop L 2 . The outer subdomain A is bounded by both of these two loops and the inner sub-domain B is formed only by the inner loop L 2 .
The characteristic point of the geometry is defined in this paper as the point that is shared by three or more boundary edges. In general, the characteristic point is connected to more than two sub-domains.
Based on the characteristic point concept this paper establishes a segmentation rule for the geometric boundaries as follows: (1) The straight edge is considered a single segment itself; (2) The end-to-end curved edges are integrated to form a whole segment; (3) If a characteristic point is encountered, its connected edges are regarded as belonging to different segments. For the geometric model in Fig. 1, all of the edges on loop L 1 are curved and form a whole segment, while the edges on loop L 2 are straight and each is treated as an individual segment. Figure 2 shows the flow chart for describing the whole Delaunay triangulation process, including a total of nine steps. The techniques for generating boundary points and inserting new nodes are two key research aspects.

The Generation of Boundary Points
The boundary points are generated according to the required point spacing specified by users. For an arbitrary edge b i (i= 1, 2, …, nb), suppose that its two endpoints are P a and P b with the required point spacing as ps(P a ) and ps(P b ).
First, find the midpoint of b i , denoted as P m . If b i is a straight edge, P m is just the midpoint of P P a b . If b i belongs to a curved segment, P m is obtained by interpolation at u m = (u a + u b )/2 using Eqs. (2), (3), and (4).
Second, compute the practical point spacing at P m using: Third, determine the lower value between ps(P a ) and ps(P b ), denoted as ps(min).
Next, check whether psl (P m ) is larger than ps(min). If yes, accept P m . If no, reject it. If P m is accepted, a required  point spacing, ps(P m ), needs to be assigned to it using the following interpolation equation: The method stated above is used to generate boundary points for the geometry shown in Fig. 1 (Fig. 3a). The number of the boundary points is increased to 67 and the number of boundary edges is increased to 78. Figure 3b provides the current boundary contour of the geometry, where the black points represent the boundary points on loop L 2 .

The Insertion of New Points
Since each Voronoi vertex corresponds to a Delaunay triangle, a new point insertion can be realized through the destruction and reconstruction of the Voronoi structures. The whole process is divided into the following four steps (Bowyer 1981): (1) Identify the Delaunay triangle containing the new point using the right-hand rule. For an arbitrary Delaunay triangle the introduced point forms three new triangles with its three edges. Next, check whether these three new triangles have the same orientations. If yes, the introduced point is inside the selected triangle. If no, the introduced point is outside the selected triangle and located on the side with the different orientation.
(2) Determine the Voronoi vertices to be deleted. According to the in-circle criterion the Voronoi vertices whose circum-circles contain the newly inserted point should be deleted. (3) Construct a local Voronoi diagram and Delaunay triangles. When the Voronoi vertices identified in step (2) are deleted, an empty convex polygon appears in the deleted region. The edges forming the empty polygon are the residual edges of deleted triangles in practice. Every residual edge gives rise to a new Voronoi vertex with the introduced point. (4) Merge the local Voronoi diagram into the global diagram. Figure 3c shows the resulting mesh after inserting all of the boundary points into the Voronoi diagram. Figure 3d shows the final mesh after generating and smoothing the interior points, which is used as the background grid to produce quadrilateral mesh in the subsequent procedures.

Key Techniques in Quadrilateral Mesh Generation
The quadrilateral mesh is generated indirectly by combining adjacent pairs of triangles. Figure 4 shows the whole quadrilateral mesh generation process, consisting mainly of four aspects: the initial combination of triangles, the treat-ment of residual triangles, topological optimization, and node smooth.

The Initial Combination of Triangles
In order to ensure that there are as many triangles combined into quadrilaterals, this paper employed the searching advancing-front edges method to create quadrilaterals layer by layer (Owen et al. 1998(Owen et al. , 1999.
The inherent boundaries of the input geometry are regarded as the first layer of advancing-front edges for quadrilateral creation. Each appropriate pair of adjacent triangles is combined into a quadrilateral by taking the advancing-front edge as the base. The possible combination patterns based on each advancing-front edge provide two options: between the two connected triangles sharing the other two edges of the base triangle, select one to form a quadrilateral.
There are many factors affecting the selection of combination patterns, such as the priority levels of the elements, the number and positions of advancing-front edges, etc. This paper made an intensive research on the various types of connection relationships among connected triangles. Five combination patterns aimed at three cases were established (Fig. 5) (Sun et al. 2015). The first case (Fig. 5a) is suitable for the triangle with two advancing-front edges. Since the base triangle fit for this case has only one combination candidate, it possesses the highest priority level. This case should be identified and converted first so as to avoid omission during combination. The second case (Fig. 5b) aims at the situation where two adjacent triangles contain two connected advancing-front edges. The third case (Fig. 5c) is suitable for the common base triangle with only one advancing-front edge.

The Movement of Residual Triangles
In the resulting mesh after initial combination, some residual triangles scattered in the quadrilaterals always remain. These residual triangles are treated and eliminated through movement and recombination. For an arbitrary residual triangle, the movement process includes the following three steps (Sun et al. 2015): (1) Find the moving target of the residual triangle, that is, the residual triangle whose position is closest to it. Compute the center of the target triangle, labeled as P m .
(2) Determine the movement path of the residual triangle.
This paper introduced the determinant value of the normalized Jacobian matrix (the scaled Jacobian) (Knupp 2003;Garimella et al. 2004) to determine the movement path of the residual triangle. Three new triangles with rotation directions are created first by the three edges of the residual triangle and the center of the target triangle, P m .
Then, compute the scaled Jacobians of these three new triangles at point P m by: where (x p , y p ) represent the coordinates of P m , (x 0 , y 0 ) and (x 1 , y 1 ) represent the starting and ending points of the three edges of the residual triangle.
Compare the three values of the scaled Jacobians computed using Eqs. (7) and (8). The edge that forms the triangle with the smallest negative value is specified as the movement path.
(3) Move and recombine the residual triangles. This operation is carried out through the destruction and reconstruction of adjacent quadrilaterals (Heighway 1983). First, find the connected triangle that shares the movement path determined in step (2).
Second, identify the quadrilateral in which the connected triangle belongs. In this case, there are two ways to combine the residual triangle: one is its connected triangle sharing the movement path, the other is the new triangle derived from the common-edge exchanging of the identified quadrilateral. Suppose that these two ways are both applied and compute the distances from the new residual triangle to P m respectively. Accept the way that yields the smaller value.
The movement process is repeated until two residual triangles are connected with a common edge and thus combined into a quadrilateral.

THE ODD AND EVEN PARITY CONTROL ON THE NUMBER OF RESIDUAL TRIANGLES
In theory, every two residual triangles in the same meshed domain are able to be recombined into a quadrilateral through moving to each other along an appropriate path. However, in practical applications, the specific problems are more complicated. The number of residual triangles and the boundary features of the geometric domains both play an important role in the feasibility and integrity of quadrilateral conversion.

The Treatment Strategies of Residual Triangles
Due to the diversity of the boundary features for the geometry, the complexity in treatment strategies of residual triangles has different degrees. In general, the number and distribution of the sub-domains have a considerable affect on the ways to handle residual triangles.

The Situation for the Geometric Model with One Sub-Domain
If the geometric model to be analyzed contains only one single domain, there will be no overlapped edge on the input boundaries. In this situation, a full conversion from triangles into quadrilaterals is able to be achieved regardless whether an odd or even number of triangles is generated in the Delaunay mesh.
When the Delaunay triangular mesh consists of an even number of triangles, the resulting mesh after adjacent triangle combination still contains an even number of residual triangles. After movement, all of the residual triangles are recombined in pairs and there is no residual triangle left in the final quadrilateral mesh When the Delaunay triangular mesh contains an odd number of triangles, the resulting mesh after initial combination still contains an odd number of residual triangles. With the movement and recombination of residual triangles, one and only one residual triangle will remain in the finally generated mesh. In this case, the only residual triangle can be treated using three stages: (1) Move it onto the nearest boundary edge of the geometry.
(2) Insert a new point on the nearest boundary edge. If this edge belongs to a curved segment, B-splines should be constructed first.
(3) Convert the only residual triangle into quadrilaterals according to the transformation template (see Fig. 6). In this way the complete conversion from triangles into quadrilaterals is realized.

The Situation for the Geometric Model with Multiple Sub-Domains
If the meshed geometry contains two or more subdomains, then the edges lying on the intersecting boundaries between neighboring sub-domains are overlapped. As proper pairs of adjacent triangles are combined into quadrilaterals, the integrity of the combination is determined to the maximum extent by the number of the triangles in each sub-domain.
After the geometry is triangulated, some sub-domains may contain even numbers of triangles and some may contain odd numbers. Since the triangles are combined in pairs, the number of the triangles that have been converted is certainly even. If a sub-domain contains an even number of Delaunay triangles, it will leave an even number of residual triangles in the resulting mesh after preliminary combination. All of the residual triangles in this sub-domain are able to be eliminated using pairwise movement and recombination. However, if a sub-domain contains an odd number of Delaunay triangles, there will be an odd number of residual triangles left in this meshed region after preliminary combination. It should be noted that the movement of each residual triangle is only restricted to the sub-domain that it belongs to. Thus, a single residual triangle will remain in this sub-domain after treatment. If this sub-domain is located on the boundaries of the geometry, the single residual triangle can be moved to the nearest outer edge and converted into quadrilaterals using the transformation template in Fig. 6. If this sub-domain is inside the geometric model, all of the boundary edges surrounding it are overlapped. In this case, while the single residual triangle is moved to the nearest overlapped edge and topologically handled, it will lead to an incompatible meshing result. This is because of new point insertion on the overlapped edge during applying the conversion template in Fig. 6. The newly inserted point cannot be simultaneously compatible with the elements in both connected sub-domains and will thereby produce unqualified meshes. Figure 7a shows the resulting mesh after preliminary combination for the triangular mesh in Fig. 3d. The figure shows that there are 17 residual triangles in the outer domain A and 3 residual triangles in the inner domain B. The num-bers of residual triangles in these two sub-domains are both odd. These residual triangles are then moved and recombined within their own sub-domains. The resulting mesh is shown in Fig. 7b. It is obvious that there is only one residual triangle left in each sub-domain. The treatment method for the remaining residual triangle in domain A is as follows: (1) Move it onto the outer boundaries of the geometric model (Fig. 7c); (2) Covert it into quadrilaterals according to the conversion template (Figs. 7d, e). However, this treatment method is not suitable for the only residual triangle in domain B. This is because the conversion template cannot ensure conformity for both sub-domains.
If the single residual triangle is located in the interior domain of the geometry, it cannot be converted into quadrilaterals using any topological optimization template under the premise of ensuring mesh conformity. For this reason, this problem is resolved by controlling the number of the residual triangles in the previous steps.

The Relationships Between the Number of Boundary Edges and the Number of Residual Triangles in Each Sub-Domain
Intensive research was conducted in this work on the relationships between the number of boundary edges and the number of mesh elements for a certain sub-domain using numerous triangle and quadrilateral meshing examples. Several crucial rules are summarized as follows: (1) If a sub-domain consists of an odd number of boundary edges it will generate an odd number of triangles during Delaunay triangulation. Similarly, if a sub-domain consists of an even number of boundary edges, it will generate an even number of Delaunay triangles.
(2) Once a new point is inserted to the current Voronoi diagram using the method described in section 2.3.2, the number of the Delaunay triangles in this sub-domain will be changed by 2n (n is an integer). In other words, the odd and even parity of the number of the triangles in each domain will not be changed with the addition of new points, regardless whether boundary points or interior points are inserted.
(3) If a sub-domain in the Delaunay triangular mesh contains an odd number of triangles, the resulting mesh after initial combination also contains an odd number of residual triangles. After movement and recombination of residual triangles, a residual triangle still remains in this sub-domain. If a sub-domain contains an even number of triangles in the Delaunay triangular mesh, the number of residual triangles in the initially combined mesh will also be even. There will be no residual triangle left in the resulting mesh after moving and recombining residual triangles.
As a result the number of boundary edges for the given geometry has a straightforward affect on a series of subsequent meshing procedures. In summary, for a certain subdomain, an odd number of boundary edges will lead to the following results directly or indirectly: (1) This sub-domain contains an odd number of triangles in the resulting mesh after inserting boundary points; (2) This sub-domain contains an odd number of triangles after inserting interior points; (3) The number of residual triangles in the mesh after initial conversion is odd; (4) There will be one residual triangle left in this sub-domain after moving and recombining the residual triangles.
Similarly, for a certain sub-domain, an even number of the boundary edges will lead to the following results directly or indirectly: (1) This sub-domain contains an even number of triangles in the resulting mesh after inserting boundary points; (2) This sub-domain contains an even number of triangles after inserting interior points; (3) The number of residual triangles in the resulting mesh after initial conversion is even; (4) There is no residual triangle left in this subdomain after treatment and all of the triangles are converted into quadrilaterals.

The Adjustment for the Number of Boundary Edges
Based on the rules summarized in section 3.2, this paper proposed an adjustment method to control the number of the residual triangles appearing in quadrilateral conversion through adjusting the number of boundary edges in the previous steps.
Since each sub-domain of the geometry may contain several boundary loops and each loop is probably shared by other sub-domains, the control on the number of boundary edges for different sub-domains may influence each other.
In order to avoid this cross interaction, this paper proposes adjusting the number of the boundary edges on each loop to even. In this way there will always be an even number of boundary edges in each sub-domain and thus the number of generated triangles will be even as well. If the number of boundary edges in a loop is odd an additional point is required to be inserted and an appropriate edge for the insertion of the new point should be determined at the same time.

The Factors Affecting the Selection of Boundary Edge
For the geometric model with multiple domains, the intersecting boundary between two adjacent sub-domains is composed of a collection of overlapped edges. The selection of the boundary edge for new point insertion is an important issue. If the boundary edge is not selected correctly, the conformity of the mesh between adjacent sub-domains will possibly be destroyed. The selection of the boundary edge is affected by two factors: (1) The position properties of boundary loops: The boundary loops for a complex geometric model surround a sub-domain formed by chains of edges lying in different positions. Some edges may be located on the outer boundaries of the geometry and only belong to this sub-domain, while other edges may lie on the intersection boundaries between two connected sub-domains and are overlapped. As the new point is inserted on the overlapped boundary edge, it must be compatible with both connected sub-domains. Thus, the new point insertion will impact the odd and even parity of the number of boundary edges in both of these sub-domains. (2) The topological connections of boundary loops: If the geometric model contains several sub-domains, different boundary loops usually contain different numbers of overlapped boundary edges and the complexity of their topological connections with adjacent sub-domains varies. When the points belonging to a boundary loop have relatively complicated topological connections the number of control methods becomes more difficult. For example, if a sub-domain is surrounded by boundary loops with one or more characteristic points, especially when they are inside the geometric model, the new point insertion process will be hindered by the mutual influence among the sub-domains connected to the characteristic points. In order to resolve this problem the boundary loop types should be identified and adjusted preferentially.

The Identification of Special Boundary Loops and Common Boundary Loops
Certain boundary loops with special geometric topologies need to be identified and adjusted preferentially. Only when these special loops are adjusted to have an even number of boundary edges should other common loops be checked.
(1) The priority adjustment for special boundary loops: The special boundary loop required for priority adjustment is identified as the loop that conforms to the following two conditions: (a) all of the boundary edges lying on it are overlapped; and (b) it includes at least one characteristic point. If a special boundary loop consists of an odd number of boundary edges, a new point must be preferentially inserted.
(2) The adjustment of the common boundary loops: When the special boundary loops required for priority treatment have been identified and adjusted to contain an even number of boundary edges, the next step is to check the odd and even parity of the number of boundary edges on the common loops. This paper proposes corresponding adjustment rules aimed at three different situations for the treatment of common boundary loops. The first rule is suitable for the case where the boundary loop has no overlapped edges and lies on the outer boundaries of the geometry. It should be noticed that, if a boundary loop contains no overlapped edge it will generally contain no characteristic point. If the boundary loop fit for this case consists of an odd number of boundary edges the new point can be inserted on any one of the edges belonging to it.
The second situation aims at the boundary loops whose boundary edges are all overlapped and have no characteristic point. Such boundary loops are located in the interior of the geometry independently. The geometric meaning is expressed as the entire intersection lines between two adjacent sub-domains. The odd and even parity control on the number of boundary edges for this type of boundary loop has an identical affect on its two connected sub-domains. If this type of boundary loop consists of an odd number of boundary edges, the new point can be inserted on any one of the edges belonging to it.
The third rule is fit for the situation where the boundary loop contains overlapped edges and also non-overlapped edges. In this case, these two types of boundary edges have different priority for new point insertion. In order to avoid bringing adverse impact on the odd and even parity control of boundary edges in the adjacent sub-domain, the nonoverlapped boundary edges have a higher priority level for new point insertion.

The Further Determination of Boundary Edge for Inserting New Points
When a boundary loop is required for adjustment, regardless whether it is a special loop or common loop, there is more than one option for the edge that is chosen for new point insertion. It is necessary to find the most appropriate edge from the candidate options.
In order to reduce the error between the practical point spacing and the required point spacing to the minimal level, this paper proposes a judgment method to determine the most appropriate edge for new point insertion based on the relative point spacing error. This method is realized using the following procedures: First, find all of the candidate boundary edges that belong to the certain loop and are able to be selected for new point insertion. Those are all of the boundary edges on special loops or the first two cases of common loops, and the non-overlapped edges in the third case of common loops.
When all of the candidate boundary edges are identified in the first step, compute their relative density errors. The relative density error of each edge is defined by the point spacing error of its two endpoints. Assuming the starting and ending points of an edge are labeled as e 0 and e 1 , its relative density error is computed using the following equation: where ps(e 0 ) and ps(e 1 ) represent the required point spacing of points e 0 and e 1 , and psl (e 0 ) and psl (e 1 ) represent their practical point spacing. |e 0 e 1 | indicates the length of the boundary edge. Among all candidate edges sort out the one that obtains the smallest relative density error value. This edge is selected as the edge where the new point is generated. Figure 8 shows the flow chart for describing the whole boundary edge number control process for each loop, where i represents the current loop to be adjusted and nl represents the number of boundary loops. It can be seen that boundary point generation contains two aspects: the first is based on the required point spacing, and the second is based on the number control.
For the geologic model shown in Fig. 1, the outer loop L 1 contains 8 boundary edges and the inner loop L 2 contains 4 boundary edges. After the boundary points are generated based on the required point spacing, the number of edges on loop L 1 is increased to an even number, 56, and that on loop L 2 is increased to an odd number, 11 (see Figs. 3a,b). That is to say, there are a total of 67 boundary edges in the outer sub-domain A (bounded by loops L 1 and L 2 ) and 11 boundary edges in the inner sub-domain B (bounded by loop L 2 ). Each of these two sub-domains consists of an odd number of boundary edges. In this case, a residual triangle will remain in each sub-domain of the resulting mesh with initial conversion (Fig. 7b). The adjustment method proposed in this paper is applied to resolve this problem. Since loop L 2 contains an odd number of edges, a new point P is inserted on it to change the odd and even property (Fig. 9a). The number of the edges in loop L 2 is increased to 12, an even number. A total of 68 boundary edges in sub-domain A and 12 boundary edges in sub-domain B are produced. Both of these two sub-domains contain even numbers of boundary edges. The Delaunay triangular mesh with number adjustment is shown in Fig. 9b. The resulting mesh after quadrilateral conversion is shown in Fig. 9c, with 18 residual triangles in sub-domain A and 2 residual triangles in sub-domain B. The numbers of residual triangles in these two sub-domains are both even and can be recombined. The residual triangles are handled and the mesh quality is improved (Figs. 9d, e). Neither of the sub-domains contains a residual triangle and all of the triangles are converted into quadrilaterals.

APPLICATIONS
The meshing generator developed in this work is used in this section to create quadrilateral meshes for complex geometries using several practical engineering application sub-domains. The adjustment method proposed is applied to implement full quadrilateral conversion by controlling the odd and even number parity of the boundary edges on each loop.
The interaction between one-dimensional river and two-dimensional overland flows involves two cases: the first occurs between overland and river nodes, and the second between overland and junction nodes. Each river node is associated with two overland nodes. Each junction node is associated with several overland nodes. Two treatment strategies are used to handle this problem realized by establishing the mapping relationships between one-dimensional river and two-dimensional overland. The first treatment strategy (Fig. 10) locates the two mapping nodes belonging to two separate regions on the same position as the original river node. The second treatment (Fig. 11) locates the two mapping nodes for each river point on their practical positions on the connected overland. Figure 10 shows an actual example to resolve realistic geology and hydrology problems for the Kaoping River drainage basin in Taiwan. Figure 10a shows the boundary map of the analyzed district captured from satellite images. The meshed region consists of four domains separated by river reaches intersecting at two junctions, J 1 and J 2 (two characteristic points of the geometry). Each domain should be meshed separately and the points on the intersecting boundaries between adjacent domains must be conformal. As boundary points are generated, each loop is adjusted to have an even number of edges and the overlapped edges have the lower priority level for adding new points. Figure 10b provides the triangular mesh generated using the Delaunay triangulation method, containing 10186 nodes and 20071 triangles. Appropriate pairs of adjacent triangles are combined into quadrilaterals based on the conversion patterns in Fig. 5. Figure 10c shows the resulting mesh after initial combination. Each sub-domain contains an even number of residual triangles. In this case all of the residual triangles can be moved and recombined into quadrilaterals. Figure 10d provides the final mesh after treating residual triangles and quality improvement, consisting of 10222 nodes and 9823 quadrilaterals in all. This mesh can accurately describe the boundary features of the input geometry and achieve the full conversion of quadrilaterals. It should be noted that each of the nodes on the river reaches is associated with two overland nodes. Figure 11 provides a practical quadrilateral meshing problem for handling realistic geology, hydrology, and water resources for the Choshui River of Taiwan. The boundary map of the analyzed district is captured from satellite images (Fig. 11a), containing seven meshing regions with three river reaches intersecting at a junction. While generating boundary points each loop is adjusted to have an even number of edges. Figure 11b shows the triangular mesh produced by the Delaunay triangulation method based on the specified point spacing, containing a total of 14109 nodes and 27480 triangles. Note that each mesh domain consists of an even number of triangles. The quadrilateral mesh is created based on this triangular mesh by combining appropriate pairs of connected triangles in each domain (Fig. 11c). In each domain an even number of residual triangles remain. Figure 11d shows the final quadrilateral mesh after residual triangle treatment and quality improvement, consisting of 14133 nodes and 13430 quadrilaterals in all. All of the elements in the finally obtained mesh are quadrilaterals without any residual triangles left. If additional quadrilateral meshes are requested to be generated in the river and junction area, the resulting mesh as shown in Fig. 11e is achieved, where the amplified view shows the local mesh in the vicinity of the junction where three river reaches meet. This mesh includes 14283 nodes and 14070 elements, with all of the blank areas completely covered by quadrilaterals.

CONCLUSIONS
This paper developed an automatic mesh generation system for numerical simulations in the fields of geology, hydrology, and water resources using C++ and Fortran programming. The developed mesh generator serves as a preprocessor for finite element models in solving two-di-mensional subsurface flow and water transport problems. It is designed to generate three-point triangular or four-point quadrilateral elements for two-dimensional domains.
The key techniques in the meshing procedures were studied in detail, such as the generation of boundary points, the insertion of new points, combination of triangles and the treatment of residual triangles, etc.
(1) Five combination patterns aimed at three cases were proposed to convert a triangular mesh into a quadrilateral mesh. The application of these patterns realized that there were as many triangles combined into quadrilaterals during the initial conversion. The quality of the generated quadrilaterals was ensured well.
(2) The treatment strategies for residual triangles for arbitrary types of geometric models were researched and presented. Intensive research was conducted on the effect for boundary edges surrounding each sub-domain on the subsequent meshing procedures.
(3) For the geometry with several sub-domains, an adjustment method was proposed to control the number of boundary edges on each loop to ensure that each sub-domain contained an even number of Delaunay triangles. This method was able to implement the full conversion from triangles into quadrilaterals. (4) The number of boundary edges on each loop was adjusted by inserting new points. The adjustment rules for special loops and three cases of common loops were presented. The most appropriate edge for new point insertion was determined by the relative point spacing error. This treatment approach could ensure reasonable density distribution for mesh elements and nodes. (5) The two-dimensional triangular and quadrilateral mesh can be stretched along the z-coordinate axis into a threedimensional triangular prism and hexahedral mesh. The proposed program was also designed to accommodate the need for aquifers with heterogeneous systems by identifying the type of material in each layer of elements.