- Email: [email protected]
Automatic triangulation over three-dimensional parametric surfaces based on advancing front method
Finite Elements in Analysis and Design 41 (2005) 892 – 910 www.elsevier.com/locate/finel
Automatic triangulation over three-dimensional parametric surfaces based on advancing front method Baohai Wu∗ , Shangjin Wang School of Energy and Power Engineering, Xi’an Jiaotong University, 710049 Xi’an, China Received 20 May 2004; accepted 4 November 2004 Available online 22 January 2005
Abstract A new mesh generation method is proposed for creating adaptive triangle meshes over three-dimensional parametric surfaces based on the advancing front method. The method is advanced from two opposite boundaries instead of the completed edges. This method solves the problem that will lead to triangulation procedure failure or quality decrease of the surface triangulation results because of poor surface corner conditions. A robust and refined triangular element formation procedure is established, which proceeds directly in three dimensions and then the element is mapped into the parametric space in order to carry out the succeeding functions in two-dimensional space, such as validity checking and convergence checking. A simple but effective convergence checking scheme is also presented in this paper and misjudgement of the convergence points can be avoided. Numerical experiments demonstrate that the triangulation results with high quality can be obtained easily by using the method. 䉷 2004 Elsevier B.V. All rights reserved. Keywords: Mesh generation; Triangulation; Three-dimensional parametric surface; Advancing front method
1. Introduction Mesh generation is a bridge connecting between geometric definition and numerical analysis of the geometric solids. Discretization and triangulation of surfaces is one of the most important research aspects in this topic. It is widely used in finite element methods (FEM), stereolitography (SLA), CAD/CAM and parametric surface rendering, etc. [1]. For an application in which a complicated domain is involved, it is quite often that the time spent on the mesh creation is longer than that in the numerical ∗ Corresponding author. Tel.: +86 29 82663777x8311.
E-mail address: [email protected] (B. Wu). 0168-874X/$ - see front matter 䉷 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2004.11.003
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
893
analysis [2]. With the rapid advancement in the areas of computer science and numerical analysis, the role of a fully automatic mesh generator becomes increasingly important. Furthermore, the result of surface triangulation is the input data of the tetrahedral grids and the quality of tetrahedral grids directly depends on the quality of the surface discretization [3]. In recent years, significant efforts have been made towards developing efficient and robust algorithms for the automatic triangulation over three-dimensional parametric surface. It appears that the following two methods are most effective and reliable to generate high quality meshes among many methods available to date [2]: (1) Delaunay based algorithms [4–7]. In this approach, a coarse mesh is usually generated to cover the problem domain at first. This coarse mesh is then systemically refined by repetitive insertions of nodes and by the modifications of the element connectivity until the required gradation effect is achieved. (2) Advancing front methods [3,8–11]. In this method, the boundary of the problem domain will be discretized first to form current front, and then new nodes and elements are generated one by one based on the current front. The front is refreshed when it is completed and the procedures mentioned above are reduplicated until the front is null, it means that the triangulation of the surface is finished. One of the advantages of the advanced front method (AFM) is that new triangle element formation is coinstantaneous with new node generation, and this advantage makes it possible to control the shape and size of the element through adjusting the location of the node [12]. However a lot of intersection checking between the new generated triangle element and existing elements must be computed in order to ensure that the triangle elements are valid. It is very difficult to finish this work in three-dimensional space because the triangle elements are always not coplanar and the new generated element perhaps penetrates through the existing elements. So a method is needed to place anisotropic triangles in two dimensions that will map back to isotropic triangles in three dimensions. For three-dimensional parametric surface, a metric map based on the first fundamental form of the surface was developed to finish this transformation [9]. However, the interior node derived in the parametric space cannot always obtain a triangle element with high quality in three-dimensional space because that some useful information, such as advancing direction and distance will not be maintained in the parametric space after metric mapping. For these reasons, this paper introduces a new approach to generate interior nodes which proceed directly in the three-dimensional space and then the new triangle is mapped into the parametric space in order to carry out validity checking and convergence checking of the new element in two-dimensional space. For regular parametric surface, its normalized parametric area is a unit square with four boundaries. The triangle elements are advanced from the four edges and two triangles are generated at each corner in the conventional AFM. It is suitable when the corner angle is about 120◦ , intersection or crack will be generated if the corner angle is far less or much more than 120◦ . This situation always leads to failure of the surface triangulation or quality decrease of the triangulation results at least. At the same time, the anisotropic triangle elements in the parametric space perhaps intersect with each other in an irregular way, which causes the difficulties of validity checking and convergence checking of the surface triangulation. Hence an improved method is presented in this paper, which is advanced from two opposite sides instead of the completed edges. Number of triangle elements at surface corners is determined by each corner angle. The method provides more flexibility to control the element shape and reduces the possibility of intersection between the triangle elements in the parametric space. In this study, a new approach for convergence checking of the surface triangulation is proposed as well, which is based on examining whether or not a point is inside a quadrangle.
894
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
This paper focuses on the triangulation of three-dimensional parametric surfaces based on the advancing front method. The rest of the paper is organized as follows: Section 2 formalizes the defects existing in the conventional advancing front method and the corresponding improvements are given. In AFM scheme, a good triangulation result can only be obtained from an optimal boundary curve discretization and this part is discussed in detail in Section 3. Section 4 describes the procedure of new node generation and new triangle element formation. Validity checking of triangle element is presented in Section 5 and convergence checking of the triangulation procedure is presented in Section 6. Section 7 gives some examples together with analysis. Finally, conclusions are discussed in Section 8.
2. Overall algorithm In any advancing front method, the algorithm begins with a set of boundary loop segments, which is defined as the initial front. Triangles are constructed from the front segments and grow towards the interior of the problem domain. Fig. 1 illustrates the conventional case using this method in the parametric space. The triangles are advanced from the four parametric edges and two triangles are generated in each corner. It is suitable when the corner angle formed by the two adjacent edges in three-dimensional space is about 120◦ . Intersection or crack will be produced if the corner angle is far less or much more than 120◦ (see Fig. 2, corners 1 and 2). Compromises will be introduced between the two corresponding triangles in order to continue the triangulation procedure. As shown in Fig. 2, the two new generated nodes which lead to intersection or crack are merged to one, as the dashed lines illustrate. Quality of the two corresponding triangles decreases consequently due to the compromise. Such phenomenon is determined by the geometric property of the surface since corner angles of the surface are not always about 120◦ . In three-dimensional space, the next front is similar to the offset of the current front, so the corner angles of the new front will not vary quickly. Accumulation of errors caused by the compromise in the imperfect corner will result in failure of the triangulation procedure or quality decrease of the triangulation results at least. In the conventional AFM, the elements formed by four edges will
3
4
v
1
2 u
Fig. 1. Conventional AFM in the parametric space.
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
4
895
3
v
1
2 u Fig. 2. Error producing at the surface corners.
ii 4
ii
i
i
3
v
1
2 u
Fig. 3. The developed AFM in the parametric space.
intersect with each other in an irregular way, which makes it very difficult to finish the validity checking and convergence checking of the triangle element, especially near the convergence area. Misjudgement of convergence nodes may be encountered and crack areas may be remained in the final triangulation results. In order to overcome these drawbacks, the conventional AFM is developed to triangulate threedimensional parametric surface in the present study. The initial front is not the four boundaries of the normalized parametric area but the two opposite edges, as shown in Fig. 3. The current front consists of two parts, as i and ii in Fig. 3. The improved method can avoid the defects mentioned above while reserving the advantages of the conventional AFM in the ideal case. In addition, the convergence checking
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
896
I i II ii
III v
iii
u (b)
(a)
Fig. 4. One triangle element is formed at the surface corner: (a) three-dimensional space and (b) parametric space.
I
II III i
ii
iii
v
u (a)
(b)
Fig. 5. Two triangle elements are formed at the surface corner: (a) three-dimensional space and (b) parametric space.
of triangulation procedure is simplified. Choosing two longer opposite boundaries from the four edges of the surface and computing the corner angles, if the corner angle 80◦ , one triangle is formed at this corner, as illustrated in Fig. 4, the front is advanced from I to II and to III in three dimensions and fronts i, ii and iii are formed in the parametric space. If the corner angle 80◦ < 140◦ , two triangles divide this corner (see Fig. 5) and if 140◦ < 180◦ , three triangles will be formed at the corner. This improvement enhances the flexibility of the triangulation at the surface corners and reduces the possibility of intersection between triangle elements in the parametric space, and therefore it is more robust than the conventional advancing front method.
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
897
3. Boundary curve discretization Cuilliere [13] presented a method to discretize the three-dimensional parametric curves. Given an arbitrary nodal spacing function Ed (x, y, z) across the Euclidean space, it represents a desired curvilinear arc length between two nodes. The total number of segments which need to be created along the curve with respect to the nodal spacing function is L ds N1 = . (1) 0 Ed (s) However the number of segments N1 obtained from Eq. (1) is not an integer in most cases, which means that the length of the last segment of the curve discretization cannot satisfy the nodal spacing function Ed (s), as shown in Fig. 6, the curve length of the last segment tn−1 tn is far less than the given nodal spacing function. Then the real total segment number of the curve N is taken to be equal to the nearest integer to N1 , and Ed (s) is adjusted accordingly so that the error is shared by all discretization segments. Thus, the optimal solution can be obtained with an increase of a fraction of N1 equal to As =
N . N1
(2)
And in each discretized segment [si , si+1 ] si+1 ds = As , Ed (s) si
(3)
As = 1 in ideal case. Assume the surface is r(u, v) = (x(u, v), y(u, v), z(u, v)) and c(t) = r(u(t), v(t)), (u0 u(t) u1 , u0 u(t) u1 ) is one of the curves on the given surface r(u, v). Then
u(t) = 21 [(u1 − u0 )t + (u1 + u0 )], v(t) = 21 [(v1 − v0 )t + (v1 + v0 )], t ∈ [−1, 1], 2 2 ds dv dv du du = E +G + 2F , dt dt dt dt dt
t2
(4)
tn-2 tn-1 tn
t1
t0 Fig. 6. Curve discretization.
898
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
E, F and G are the first fundamental quantities of the surface [14], 2 2 2 jr jr jx jy jz = + + , E= ju ju ju ju ju jr jr jx jx jy jy jz jz F= = + + , ju jv ju jv ju jv ju jv 2 2 2 jr jr jx jy jz G= = + + . jv jv jv jv jv And Eq. (3) can be rewritten as ti+1 ds 1 dt = As . Ed (t) dt ti
(5)
(6)
Romberg integration method is utilized for Eq. (6) in this study and the result is satisfying. For boundary curves, there is only one variable parameter and the other is constant, so Eq. (4) can be simplified. For the boundary v= constant, 2 ds du = E . (7) dt dt And for the boundary u= constant, 2 ds dv = G . dt dt
(8)
The following work is to quickly determine the parametric increment t = ti+1 − ti with respect to As . Cuilliere [13] obtained the increment through the Gauss integration scheme by giving a little parametric increment t. This scheme not only increases the computing efforts but also accurate solution cannot be obtained easily. Newton method is adopted in this paper, iterative and these drawbacks are avoided. t ds (t) = 1 Let f (t) = ti Ed1(t) ds dt − A , then f s dt Ed (t) dt . The implementation details of Newton iterative method are listed as follows: 1. Giving an initial value t0 and f (t0 ), f (t0 ) are calculated. 2. Computing tk+1 = tk − f (tk )/f (tk ). 3. The iteration is finished if |tk+1 − tk | < 1 or |f (tk+1 )| < 2 and t = tk+1 ; else k = k + 1, turn to 2. The curve will be discretized according to the given node spacing function Ed (x, y, z) and all discretized nodes on the surface and in the parametric space are obtained together. 4. New triangle element generation The key challenge of generating the new triangle element is to determine the new node according to the base segment of the current front. There are two methods to determine the new node, one proceeds in the parametric space [2,9] and the other proceeds in the three-dimensional space [3,10,11]. Iterative
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910 A rv m
a
899
r(u,v0)
yL M
b zL ru
xL
C c
B
r(u0 ,v)
Fig. 7. Finding the new node in the three-dimensional space.
process or time-consuming computing is needed for the latter method, and the final node is located on the tangent plane because finding an accurate node on the surface to satisfy the given condition is very difficult. Based on these facts, this paper presents a simple but effective approach to find the new node, which is calculated directly in tangent plane and then the new triangle is mapped into the parametric space. Giving an arbitrary line segment of the current front ab in the parametric space, the coordinates of the two endpoints are (ua , va ) and (ub , vb ) respectively and the corresponding points on the surface are A and B, as shown in Fig. 7. M is the middle point of the curve segment AB and m is the corresponding location in the parametric space. Let xL = AB/ AB , yL is the unit surface normal at M, and zL = xL × yL . dAB is the length of the line segment AB in the three-dimensional space. In order to obtain a high√ quality triangle element, the distance between the new node C and AB, i.e. MC , should be dMC = 23 dAB . For a regular surface, ru × rv = 0, the parametric increments of u and v from the middle point M can be determined by solving the following two linear equations (see Fig. 7): u(ru · xL ) + v(rv · xL ) = 0, u(ru · zL ) + v(rv · zL ) = dMC ,
(9) (10)
where ru and rv are the tangential vectors along u-parametric direction and v-parametric direction respectively, and ru = jr/ju, rv = jr/jv . Eq. (9) represents that the distance dMC has zero projection on the xL direction and Eq. (10) indicates the distance dMC marching along zL with dMC . By solving Eqs. (9) and (10), we can find the parameter increments u and v in the following way: u =
−dMC (rv · xL ) , (ru · xL )(rv · zL ) − (ru · zL )(rv · xL )
(11)
v =
dMC (ru · xL ) . (ru · xL )(rv · zL ) − (ru · zL )(rv · xL )
(12)
Therefore the parametric coordinates of the new node C are uc = um + u, vc = vm + v.
900
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
Connecting points A, B and C, the new triangle element in three dimensions is constructed and the corresponding triangle element in the parametric space is obtained through connecting points a, b and c, as denoted by the dashed lines in Fig. 7.
5. Validity checking of the new triangle element In AFM, validity checking of the new triangle element is essential to surface triangulation. It is possible for the new triangle element to intersect with the existing elements and redundant or low quality triangle elements will be generalized if two adjacent nodes in the new front are too close or too far. Validity checking involves two aspects, i.e. intersection checking between triangle elements and nodes merging or inserting. Both are discussed in detail in the following sections. 5.1. Intersection checking between triangle elements When a new triangle element is formed, the intersection between the new element and existing elements must be checked and it is finished in the parametric space in present study. The essence of triangles intersection in the two-dimensional space is the intersection between line segments. Li [15] provided a fast approach to estimate if the two line segments intersect with each other. Giving two planar line segments P1 P2 and Q1 Q2 (see Fig. 8), P1 Z and Q1 Z are perpendicular to the plane determined by P1 P2 and Q1 Q2 . Then the sufficient and necessary conditions for the two line segments to intersect at one point other than endpoints are ((P1 P2 × P1 Q1 ) · P1 Z)((P1 P2 × P1 Q2 ) · P1 Z) < 0
and
((Q1 Q2 × Q1 P1 ) · Q1 Z )((Q1 Q2 × Q1 P2 ) · Q1 Z ) < 0.
(13)
Z'
Q1 Z
P2
P1
Q2
Fig. 8. Intersection checking between two line segments.
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
(a)
901
(b)
Fig. 9. Intersection occurrence and modification among triangle elements: (a) intersection occurring and (b) modification for the intersection case.
When a new triangle element is formed, the two edges connected with the new node should be checked whether they intersect with all existing triangle elements in the parametric space. If the intersections occur (see Fig. 9(a)), the two nodes should be merged to one new node, and the new node can be taken as the mid-point of the two nodes, as illustrated in Fig. 9(b). 5.2. Nodes merging or inserting The criterion of nodes merging or inserting is to check the distance between two adjacent nodes on the new generated front. If the distance is larger than a given upper limit, a new node is inserted; if the distance is less than a given lower limit, these two nodes should be merged to one. As the triangulation procedure proceeds in the parametric space and the distance between two adjacent nodes is calculated in the three dimensions, the metric map is introduced based on the first fundamental quantities of the surface. It gives an approximate distance map from the three-dimensional space to the parametric space. For a point p in the parametric space, P is the corresponding point on the surface, the tangent plane Tp of the surface at P can be considered as a Euclidean space with the metric map M [2,6,9], E F M= , (14) F G where E, F and G are the first fundamental quantities of the surface as defined in Eq. (5). In the neighborhood of point P, the difference between the distance with this metric map on Tp and the corresponding curve distance on the surface is only a high order small quantity and can be neglected. Since EG − F 2 > 0 [14], the map matrix is positive definite symmetric and it can be rewritten as 1 0 [e1 e2 ]T , M = [e1 e2 ] (15) 0 2 where e1 , e2 and 1 , 2 are the eigenvectors and eigenvalues of the matrix M, respectively. If h1 = h2 , the metric is isotropic and the directions are immaterial in such a case. Otherwise, the metric is anisotropic
902
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
so that e√ 1 and e2 will define the stretching directions of the elements in the parametric space. Given hi = 1/ i (i = 1, 2), hi and ei (i = 1, 2) are then computed to form an ellipse E with the center at C. It will be a unit circle in the tangent plane Tp if it is metric mapped into the three-dimensional space, as √ illustrated in Fig. 10. For a circle with radius r in Tp , the corresponding principal radii will be hi = r/ i (i = 1, 2) in the parametric space. Let r = 0.5 AB and searching ellipse E1 is formed, if the adjacent nodes are inside E1 , these nodes should be merged into one new node and the new node can be taken as the centroid of these nodes. Let r = 1.4 AB and searching ellipse E2 is formed, if the adjacent node is outside this searching ellipse, a new node should be inserted between these two nodes, and two corresponding triangle elements are formed. As shown in Fig. 11, polygonal line 1234 is the current front and 5, 6 and 7 are the new generated nodes, E1 and E2 are the searching ellipses with the center at node 6. Fig. 11 shows that node 5 is outside the ellipse E2 and node 7 is inside the ellipse E1 . This indicates that the three-dimensional distance between node 5 and 6 on the surface is larger than the upper limit and the distance between node 7 and 6 on the surface is less than the lower limit. Therefore, a e2 h2 h1
e1
C
A
B
Fig. 10. Searching ellipse in the parametric space.
5
E2 E1
8 1 6
2
9 7
3
4
Fig. 11. Nodes merging or inserting.
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
903
new node 8 is inserted between node 5 and 6 and two corresponding triangle elements 258 and 268 are formed. The new node 8 can be taken as the middle point of the node 5 and 6. Node 6 and 7 should be merged into one node, i.e. node 9 in Fig. 11. Node 9 can also be taken as the middle point of node 6 and 7. All of the formed triangle elements in the parametric space are illustrated by the dashed lines in Fig. 11.
6. Convergence checking Convergence checking is crucial for surface triangulation. In this study, a new approach is proposed by the authors for checking triangulation convergence. The approach is based on examining whether or not a point is inside a quadrangle. Given a simple quadrangle p1 p2 p3 p4 and a point p, as illustrated in Fig. 12, a horizontal half line is made from the point p in order to identify whether the point p is inside the quadrangle p1 p2 p3 p4 . The number of intersection points between the half line and the quadrangle p1 p2 p3 p4 is calculated. If the number is odd, p is inside the quadrangle; otherwise p is outside the quadrangle [16]. The convergence checking approach of the triangulation procedure is discussed in detail as follows. As illustrated in Fig. 13, polygonal line 12345 and 1 2 3 4 are parts of the current front, and 6, 7, 8, 9, 5 , 6 , 7 are new generated nodes. The two adjacent triangle elements of the current front on one side are joined to form a quadrangle. As shown in Fig. 13, the left triangle elements can form 3 quadrangles, i.e., 1276, 2387 and 3498. Then checking whether or not the new generated nodes of the front on the other side, i.e. nodes 5 , 6 and 7 in Fig. 13, are inside these quadrangles. If the node is inside the quadrangle, the vertex of the quadrangle which is closest (in the three-dimensional space) to the node can be determined. There may be two different situations. One is that the vertex is located on the current front, as illustrated in Fig. 13, node 7 is inside the quadrangle 3498 and closest to vertex 4, then new node 7 is obtained, which is taken as the mid-point of node 7 and 4. In addition, node 3 and 8, node 4 and 9 should be treated in the same way and nodes 3 , 4 are obtained. All triangle elements associated with these nodes, i.e. 4, 7 , 3 , 8, 4 and 9 should be adjusted accordingly, as shown in Fig. 14. The other case is that the closest vertex is one of the new generated nodes, as node 6 in Fig. 13, it is inside the quadrangle 2387 and is closest to vertex 7. This situation is relatively simpler than the former case. A new node 6 is determined by the mean p1
p p4 p2 p3 Fig. 12. Examining whether or not a point is inside a quadrangle.
904
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
1 6 5'
2
1'
7 3
4
2'
6' 8
7'
5
3'
4'
9 Fig. 13. Convergence checking approach.
1 6 2
5'
1'
7 6'
3
7 '' 5
2' 3 ''
4 '' Fig. 14. The closest vertex is on current front.
value of node 6 and 7 to replace these two nodes and the corresponding triangle elements are adjusted, and the final situation is depicted in Fig. 15. After all of the nodes have been checked, the one-side front which has formed quadrangles is fixed and the other-side front is refreshed while the nodes which are inside the quadrangles should be removed from the current front, such as node 6 and 7 . The procedure discussed above is repeated until the front is null, which indicates that the triangulation procedure is finished. Two problems may arise during the above procedure, and they must be dealt with properly. One is that triangle element with zero area may be generated during the convergence checking process. This kind of element should be discarded. The other is that some nodes on the one-side front which is fixed may be omitted, as node 7 in Fig. 16(a). The adjacent nodes of 7 are 6 and 8 and the mid-point of 6 and 8 is calculated. Then node 7 is replaced by the new node 7 , which is determined by the mean value of node 7 and the above mid-point. Two new triangle elements are formed accordingly, i.e. 67 2 and 7 82 in Fig. 16(b).
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
905
1 6 5'
2
1'
6 ''
3
2' 3 ''
7 ''
4 ''
5
Fig. 15. The closest vertex is one of the new generated nodes.
1 6 2
1'
7 3
2'
8
4
3'
5
4'
(a)
2
6 7'
3
8
2'
(b)
Fig. 16. Node is omitted when the triangulation is finished: (a) node 7 is ommitted and (b) two new triangle elements are formed.
7. Examples As an example, triangulation is implemented for a centrifugal compressor impeller. The impeller consists of 12 blades and one hub, as shown in Fig. 17. The blade is a three-dimensional free-form surface and the hub is a revolution surface. To assess the mesh quality, quality evaluation indicators are introduced. For a triangle ABC, quality is used to evaluate the individual element, 1 and 2 quality are used to assess the total triangle elements [12]. √
=2 3 1 =
CA 2
N T
i=1 i
NT
,
CA × CB , + AB 2 + BC 2 NT
2 = N T
i=1 1/i
.
(16) (17)
In the ideal case, = 1 and = 0 if A, B and C are co-linear. The hub is first triangulated. Chosen two circular boundaries as the initial fronts, the triangulation is advanced to the interior from these two boundaries. The results are shown in Fig. 18. Total 2350 triangles are created, 1 = 0.9885 and 2 = 0.9802. For this revolution surface, its normalized parametric area
906
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
Fig. 17. Centrifugal compressor impeller.
Fig. 18. Triangulation results of the hub.
is a unit square with four boundaries, u ∈ [0, 1], v ∈ [0, 1]. If triangulation was advanced from these four boundaries according to the conventional AFM, a redundant boundary line would be generalized which is actually nonexistent for this revolution surface. Triangulation will be much harder because of the boundary and the corner angles brought by this nonexistent boundary. However this drawback can be overcome by advancing triangulation from the real two circular boundaries.
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
907
The blade is a free-form surface and is modeled by B-spline, as shown in Fig. 19. Boundary 1 and 3 are u-parametric curves corresponding to u ∈ [0, 1], v = 0 and u ∈ [0, 1], v = 1, respectively, 2 and 4 are v-parametric curves corresponding to u = 1, v ∈ [0, 1] and u = 0, v ∈ [0, 1], respectively. The curve lengths of boundaries 1, 2, 3, and 4 are 47.1051, 329.9315, 119.3394 and 227.6364 mm, respectively by integrating Eqs. (7) and (8). The nodal spacing function Ed (x, y, z) is set as a constant of 5 mm and the actual nodal spacing distribution is 5.2339, 4.9990, 4.9725 and 4.9486 mm, respectively. The boundaries discretization results are shown in Fig. 20. Fig. 21 shows the triangulation results of the blade surface. The total number of triangle elements in Fig. 21 is 1702, 1 = 0.9942 and 2 = 0.9869. Majority of triangle elements with poor quality exist near the boundary 1, which is determined by the surface geometry property. For this blade surface, the corner angles formed by boundaries 1 and 2, 1 and 4 are both some 90◦ ; the angle formed by boundaries 3 and 4 is some 120◦ ; and the last angle formed by 2 and 3 is about 60◦ . It is obvious that the quality of the triangle elements near the boundary 1 is poor. The computed values of 1 and 2 are 0.8217 and 0.8125 and the minimum value of is 0.7324 in this area. Quality of the triangle elements near boundaries with poor corner conditions can be improved through subdivisions of the triangle elements. 2
1 4
Fig. 19. The centrifugal compressor blade.
Fig. 20. The boundaries discretization results of the blade.
Fig. 21. Triangulation results of the blade.
3
908
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
1
2
3
Fig. 22. The triangle status in the parametric space for trimmed surface.
Fig. 23. Triangulation results of the trimmed blade.
Triangulation of trimmed surface is also important in many branches of CAD/CAM and Graphics. A trimmed parametric surface consists of a parametric surface and a set of properly ordered trimming curves lying within the parametric rectangle of the surface. There are three statuses of triangles in the parametric space under the procedure of surface triangulation, i.e. IN, OUT and ON [17]. As shown in Fig. 22, the circle is the trimming curve in the parametric rectangle. Status IN means the triangle is in the valid parametric region, as triangle 1 in Fig. 22; status OUT means the triangle is in the trimmed area completely, as triangle 2; and status ON (triangle 3) means that the triangle intersects with the trimming curve. The triangles with status IN are reserved and the triangles with status OUT are discarded. The triangles with status ON will be trimmed so that the parts inside the valid parametric region are reserved and the others are discarded. The reserved parts of these triangles will be merged or divided according to their states respectively [18]. In order to simplify the computation, the trimming curve may be approximated by the polygonal line. Fig. 23 shows the triangulation results of a trimmed blade with the trimming curves consisting of a circle and a freeform curve in the parametric region. Total 1450 triangles are generalized and 1 = 0.9432, 2 = 0.9023. For the trimmed surface, the triangulation results of the whole initial surface can be regarded as the initial mesh and trimming of the initial mesh can be looked upon the post treatment, therefore the whole procedure is fast and robust.
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
909
8. Conclusions and discussion A developed advancing front method is presented to triangulate the three-dimensional parametric surface. The details of implementation, validity and convergence checking are given in the paper. The method is advanced from two longer opposite boundaries and the number of triangle elements at surface corner is determined by the surface corner angle. This developed method solves the problem caused by poor surfaces corner condition in the conventional AFM and increases the flexibility on the corner and boundary treatment. The new node is generated directly in three dimensions and the validity checking of the new triangle element is performed in the parametric space, therefore the advantages of AFM in two dimensions are retained sufficiently. A simple but effective convergence checking approach is presented as well, which is based on examining whether a point is inside a quadrangle or not. Numerical experiments show that the developed method is robust and easy to carry out and high quality triangulation results have been achieved.
Acknowledgements The work reported in this paper was supported by the Research Fund for the Doctoral Program of Higher Education of China under Grant No. 20020698042, and the National Science Foundation of China under Grant No. 50146014 (2002). References [1] M. Vigo, N.P. Garciia, P.B. Crosa, Directional adaptive surface triangulation, Comput. Aided Geom. Design 16 (1999) 107–126. [2] C.K. Lee, Automatic adaptive mesh generation using metric advancing front approach, Eng. Comput. 16 (1999) 230–263. [3] J.C. Cuilliere, An adaptive method for the automatic triangulation of 3D parametric surfaces, Comput. Aided Design 30 (2) (1998) 139–149. [4] P.L. Baehmann, S.L. Wittchen, M.S. Shephard, Robust, geometrically based, automatic two-dimensional mesh generation, Internat. J. Numer. Methods Engng. 24 (1987) 1043–1078. [5] T.P. Fang, L. Piegl, Delaunay triangulation using a uniform grid, IEEE Comput. Graph. Appl. 13 (1993) 36–47. [6] H. Borouchaki, P.L. George, F. Hecht, Delaunay mesh generation governed by metric specifications, Part I, Algorithms, Finite Elements Anal. Des. 25 (1997) 61–83. [7] M. Vigo, N. Pla, Computing directional constrained Delaunay triangulations, Comput. Graph. 24 (2000) 181–190. [8] P.L. George, E. Seveno, The advancing front mesh generation method revisited, Internat. J. Numer. Methods Engng. 37 (1994) 3605–3619. [9] J.R. Tristano, S.J. Owen, S.A. Canann, Advancing front surface mesh generation in parametric space using a Riemannian surface definition, in: Proceedings of the Seventh International Meshing Roundtable, Dearborn MI, USA, 1998. [10] T.S. Lau, S.H. Lo, Finite element mesh generation over analytical curved surfaces, Comput. & Structures 59 (2) (1996) 301–309. [11] V. Francois, J.C. Cuilliere,Automatic mesh pre-optimization based on the geometric discretization error,Adv. Eng. Software 31 (2000) 763–774. [12] Z.Q. Guan, C. Song, Y.X. Gu, Recent advances of research on finite element mesh generation methods, J. Comput. Aided Des. Comput. Graph. 15 (1) (2003) 1–14. [13] J.C. Cuilliere, A direct method for the automatic discretization of 3D parametric curves, Comput. Aided Design 29 (1997) 639–647.
910
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
[14] D.J. Meng, K. Liang, Differential Geometric, Science Press, Beijing, China, 2002. [15] Y.H. Li, Y.X. Hua, A method to quickly accessing the intersection of line segments, Bull. Surve. Mapping Eng. (2003) 30–31. [16] P.D. Zhou, Computational Geometric—Algorithms Analysis and Design, Tsinghua University Press, Beijing, 2000. [17] L.A. Piegl, W. Tiller, Geometry-based triangulation of trimmed NURBS surfaces, Comput. Aided Design 30 (1) (1998) 11–18. [18] M.W. Ji, C.G.Yang, J.G. Sun, A fast algorithm for triangulating the parametric space of trimmed NURBS surfaces, Chinese J. Comput. 19 (6) (1996) 450–456.
Automatic triangulation over three-dimensional parametric surfaces based on advancing front method Baohai Wu∗ , Shangjin Wang School of Energy and Power Engineering, Xi’an Jiaotong University, 710049 Xi’an, China Received 20 May 2004; accepted 4 November 2004 Available online 22 January 2005
Abstract A new mesh generation method is proposed for creating adaptive triangle meshes over three-dimensional parametric surfaces based on the advancing front method. The method is advanced from two opposite boundaries instead of the completed edges. This method solves the problem that will lead to triangulation procedure failure or quality decrease of the surface triangulation results because of poor surface corner conditions. A robust and refined triangular element formation procedure is established, which proceeds directly in three dimensions and then the element is mapped into the parametric space in order to carry out the succeeding functions in two-dimensional space, such as validity checking and convergence checking. A simple but effective convergence checking scheme is also presented in this paper and misjudgement of the convergence points can be avoided. Numerical experiments demonstrate that the triangulation results with high quality can be obtained easily by using the method. 䉷 2004 Elsevier B.V. All rights reserved. Keywords: Mesh generation; Triangulation; Three-dimensional parametric surface; Advancing front method
1. Introduction Mesh generation is a bridge connecting between geometric definition and numerical analysis of the geometric solids. Discretization and triangulation of surfaces is one of the most important research aspects in this topic. It is widely used in finite element methods (FEM), stereolitography (SLA), CAD/CAM and parametric surface rendering, etc. [1]. For an application in which a complicated domain is involved, it is quite often that the time spent on the mesh creation is longer than that in the numerical ∗ Corresponding author. Tel.: +86 29 82663777x8311.
E-mail address: [email protected] (B. Wu). 0168-874X/$ - see front matter 䉷 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2004.11.003
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
893
analysis [2]. With the rapid advancement in the areas of computer science and numerical analysis, the role of a fully automatic mesh generator becomes increasingly important. Furthermore, the result of surface triangulation is the input data of the tetrahedral grids and the quality of tetrahedral grids directly depends on the quality of the surface discretization [3]. In recent years, significant efforts have been made towards developing efficient and robust algorithms for the automatic triangulation over three-dimensional parametric surface. It appears that the following two methods are most effective and reliable to generate high quality meshes among many methods available to date [2]: (1) Delaunay based algorithms [4–7]. In this approach, a coarse mesh is usually generated to cover the problem domain at first. This coarse mesh is then systemically refined by repetitive insertions of nodes and by the modifications of the element connectivity until the required gradation effect is achieved. (2) Advancing front methods [3,8–11]. In this method, the boundary of the problem domain will be discretized first to form current front, and then new nodes and elements are generated one by one based on the current front. The front is refreshed when it is completed and the procedures mentioned above are reduplicated until the front is null, it means that the triangulation of the surface is finished. One of the advantages of the advanced front method (AFM) is that new triangle element formation is coinstantaneous with new node generation, and this advantage makes it possible to control the shape and size of the element through adjusting the location of the node [12]. However a lot of intersection checking between the new generated triangle element and existing elements must be computed in order to ensure that the triangle elements are valid. It is very difficult to finish this work in three-dimensional space because the triangle elements are always not coplanar and the new generated element perhaps penetrates through the existing elements. So a method is needed to place anisotropic triangles in two dimensions that will map back to isotropic triangles in three dimensions. For three-dimensional parametric surface, a metric map based on the first fundamental form of the surface was developed to finish this transformation [9]. However, the interior node derived in the parametric space cannot always obtain a triangle element with high quality in three-dimensional space because that some useful information, such as advancing direction and distance will not be maintained in the parametric space after metric mapping. For these reasons, this paper introduces a new approach to generate interior nodes which proceed directly in the three-dimensional space and then the new triangle is mapped into the parametric space in order to carry out validity checking and convergence checking of the new element in two-dimensional space. For regular parametric surface, its normalized parametric area is a unit square with four boundaries. The triangle elements are advanced from the four edges and two triangles are generated at each corner in the conventional AFM. It is suitable when the corner angle is about 120◦ , intersection or crack will be generated if the corner angle is far less or much more than 120◦ . This situation always leads to failure of the surface triangulation or quality decrease of the triangulation results at least. At the same time, the anisotropic triangle elements in the parametric space perhaps intersect with each other in an irregular way, which causes the difficulties of validity checking and convergence checking of the surface triangulation. Hence an improved method is presented in this paper, which is advanced from two opposite sides instead of the completed edges. Number of triangle elements at surface corners is determined by each corner angle. The method provides more flexibility to control the element shape and reduces the possibility of intersection between the triangle elements in the parametric space. In this study, a new approach for convergence checking of the surface triangulation is proposed as well, which is based on examining whether or not a point is inside a quadrangle.
894
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
This paper focuses on the triangulation of three-dimensional parametric surfaces based on the advancing front method. The rest of the paper is organized as follows: Section 2 formalizes the defects existing in the conventional advancing front method and the corresponding improvements are given. In AFM scheme, a good triangulation result can only be obtained from an optimal boundary curve discretization and this part is discussed in detail in Section 3. Section 4 describes the procedure of new node generation and new triangle element formation. Validity checking of triangle element is presented in Section 5 and convergence checking of the triangulation procedure is presented in Section 6. Section 7 gives some examples together with analysis. Finally, conclusions are discussed in Section 8.
2. Overall algorithm In any advancing front method, the algorithm begins with a set of boundary loop segments, which is defined as the initial front. Triangles are constructed from the front segments and grow towards the interior of the problem domain. Fig. 1 illustrates the conventional case using this method in the parametric space. The triangles are advanced from the four parametric edges and two triangles are generated in each corner. It is suitable when the corner angle formed by the two adjacent edges in three-dimensional space is about 120◦ . Intersection or crack will be produced if the corner angle is far less or much more than 120◦ (see Fig. 2, corners 1 and 2). Compromises will be introduced between the two corresponding triangles in order to continue the triangulation procedure. As shown in Fig. 2, the two new generated nodes which lead to intersection or crack are merged to one, as the dashed lines illustrate. Quality of the two corresponding triangles decreases consequently due to the compromise. Such phenomenon is determined by the geometric property of the surface since corner angles of the surface are not always about 120◦ . In three-dimensional space, the next front is similar to the offset of the current front, so the corner angles of the new front will not vary quickly. Accumulation of errors caused by the compromise in the imperfect corner will result in failure of the triangulation procedure or quality decrease of the triangulation results at least. In the conventional AFM, the elements formed by four edges will
3
4
v
1
2 u
Fig. 1. Conventional AFM in the parametric space.
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
4
895
3
v
1
2 u Fig. 2. Error producing at the surface corners.
ii 4
ii
i
i
3
v
1
2 u
Fig. 3. The developed AFM in the parametric space.
intersect with each other in an irregular way, which makes it very difficult to finish the validity checking and convergence checking of the triangle element, especially near the convergence area. Misjudgement of convergence nodes may be encountered and crack areas may be remained in the final triangulation results. In order to overcome these drawbacks, the conventional AFM is developed to triangulate threedimensional parametric surface in the present study. The initial front is not the four boundaries of the normalized parametric area but the two opposite edges, as shown in Fig. 3. The current front consists of two parts, as i and ii in Fig. 3. The improved method can avoid the defects mentioned above while reserving the advantages of the conventional AFM in the ideal case. In addition, the convergence checking
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
896
I i II ii
III v
iii
u (b)
(a)
Fig. 4. One triangle element is formed at the surface corner: (a) three-dimensional space and (b) parametric space.
I
II III i
ii
iii
v
u (a)
(b)
Fig. 5. Two triangle elements are formed at the surface corner: (a) three-dimensional space and (b) parametric space.
of triangulation procedure is simplified. Choosing two longer opposite boundaries from the four edges of the surface and computing the corner angles, if the corner angle 80◦ , one triangle is formed at this corner, as illustrated in Fig. 4, the front is advanced from I to II and to III in three dimensions and fronts i, ii and iii are formed in the parametric space. If the corner angle 80◦ < 140◦ , two triangles divide this corner (see Fig. 5) and if 140◦ < 180◦ , three triangles will be formed at the corner. This improvement enhances the flexibility of the triangulation at the surface corners and reduces the possibility of intersection between triangle elements in the parametric space, and therefore it is more robust than the conventional advancing front method.
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
897
3. Boundary curve discretization Cuilliere [13] presented a method to discretize the three-dimensional parametric curves. Given an arbitrary nodal spacing function Ed (x, y, z) across the Euclidean space, it represents a desired curvilinear arc length between two nodes. The total number of segments which need to be created along the curve with respect to the nodal spacing function is L ds N1 = . (1) 0 Ed (s) However the number of segments N1 obtained from Eq. (1) is not an integer in most cases, which means that the length of the last segment of the curve discretization cannot satisfy the nodal spacing function Ed (s), as shown in Fig. 6, the curve length of the last segment tn−1 tn is far less than the given nodal spacing function. Then the real total segment number of the curve N is taken to be equal to the nearest integer to N1 , and Ed (s) is adjusted accordingly so that the error is shared by all discretization segments. Thus, the optimal solution can be obtained with an increase of a fraction of N1 equal to As =
N . N1
(2)
And in each discretized segment [si , si+1 ] si+1 ds = As , Ed (s) si
(3)
As = 1 in ideal case. Assume the surface is r(u, v) = (x(u, v), y(u, v), z(u, v)) and c(t) = r(u(t), v(t)), (u0 u(t) u1 , u0 u(t) u1 ) is one of the curves on the given surface r(u, v). Then
u(t) = 21 [(u1 − u0 )t + (u1 + u0 )], v(t) = 21 [(v1 − v0 )t + (v1 + v0 )], t ∈ [−1, 1], 2 2 ds dv dv du du = E +G + 2F , dt dt dt dt dt
t2
(4)
tn-2 tn-1 tn
t1
t0 Fig. 6. Curve discretization.
898
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
E, F and G are the first fundamental quantities of the surface [14], 2 2 2 jr jr jx jy jz = + + , E= ju ju ju ju ju jr jr jx jx jy jy jz jz F= = + + , ju jv ju jv ju jv ju jv 2 2 2 jr jr jx jy jz G= = + + . jv jv jv jv jv And Eq. (3) can be rewritten as ti+1 ds 1 dt = As . Ed (t) dt ti
(5)
(6)
Romberg integration method is utilized for Eq. (6) in this study and the result is satisfying. For boundary curves, there is only one variable parameter and the other is constant, so Eq. (4) can be simplified. For the boundary v= constant, 2 ds du = E . (7) dt dt And for the boundary u= constant, 2 ds dv = G . dt dt
(8)
The following work is to quickly determine the parametric increment t = ti+1 − ti with respect to As . Cuilliere [13] obtained the increment through the Gauss integration scheme by giving a little parametric increment t. This scheme not only increases the computing efforts but also accurate solution cannot be obtained easily. Newton method is adopted in this paper, iterative and these drawbacks are avoided. t ds (t) = 1 Let f (t) = ti Ed1(t) ds dt − A , then f s dt Ed (t) dt . The implementation details of Newton iterative method are listed as follows: 1. Giving an initial value t0 and f (t0 ), f (t0 ) are calculated. 2. Computing tk+1 = tk − f (tk )/f (tk ). 3. The iteration is finished if |tk+1 − tk | < 1 or |f (tk+1 )| < 2 and t = tk+1 ; else k = k + 1, turn to 2. The curve will be discretized according to the given node spacing function Ed (x, y, z) and all discretized nodes on the surface and in the parametric space are obtained together. 4. New triangle element generation The key challenge of generating the new triangle element is to determine the new node according to the base segment of the current front. There are two methods to determine the new node, one proceeds in the parametric space [2,9] and the other proceeds in the three-dimensional space [3,10,11]. Iterative
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910 A rv m
a
899
r(u,v0)
yL M
b zL ru
xL
C c
B
r(u0 ,v)
Fig. 7. Finding the new node in the three-dimensional space.
process or time-consuming computing is needed for the latter method, and the final node is located on the tangent plane because finding an accurate node on the surface to satisfy the given condition is very difficult. Based on these facts, this paper presents a simple but effective approach to find the new node, which is calculated directly in tangent plane and then the new triangle is mapped into the parametric space. Giving an arbitrary line segment of the current front ab in the parametric space, the coordinates of the two endpoints are (ua , va ) and (ub , vb ) respectively and the corresponding points on the surface are A and B, as shown in Fig. 7. M is the middle point of the curve segment AB and m is the corresponding location in the parametric space. Let xL = AB/ AB , yL is the unit surface normal at M, and zL = xL × yL . dAB is the length of the line segment AB in the three-dimensional space. In order to obtain a high√ quality triangle element, the distance between the new node C and AB, i.e. MC , should be dMC = 23 dAB . For a regular surface, ru × rv = 0, the parametric increments of u and v from the middle point M can be determined by solving the following two linear equations (see Fig. 7): u(ru · xL ) + v(rv · xL ) = 0, u(ru · zL ) + v(rv · zL ) = dMC ,
(9) (10)
where ru and rv are the tangential vectors along u-parametric direction and v-parametric direction respectively, and ru = jr/ju, rv = jr/jv . Eq. (9) represents that the distance dMC has zero projection on the xL direction and Eq. (10) indicates the distance dMC marching along zL with dMC . By solving Eqs. (9) and (10), we can find the parameter increments u and v in the following way: u =
−dMC (rv · xL ) , (ru · xL )(rv · zL ) − (ru · zL )(rv · xL )
(11)
v =
dMC (ru · xL ) . (ru · xL )(rv · zL ) − (ru · zL )(rv · xL )
(12)
Therefore the parametric coordinates of the new node C are uc = um + u, vc = vm + v.
900
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
Connecting points A, B and C, the new triangle element in three dimensions is constructed and the corresponding triangle element in the parametric space is obtained through connecting points a, b and c, as denoted by the dashed lines in Fig. 7.
5. Validity checking of the new triangle element In AFM, validity checking of the new triangle element is essential to surface triangulation. It is possible for the new triangle element to intersect with the existing elements and redundant or low quality triangle elements will be generalized if two adjacent nodes in the new front are too close or too far. Validity checking involves two aspects, i.e. intersection checking between triangle elements and nodes merging or inserting. Both are discussed in detail in the following sections. 5.1. Intersection checking between triangle elements When a new triangle element is formed, the intersection between the new element and existing elements must be checked and it is finished in the parametric space in present study. The essence of triangles intersection in the two-dimensional space is the intersection between line segments. Li [15] provided a fast approach to estimate if the two line segments intersect with each other. Giving two planar line segments P1 P2 and Q1 Q2 (see Fig. 8), P1 Z and Q1 Z are perpendicular to the plane determined by P1 P2 and Q1 Q2 . Then the sufficient and necessary conditions for the two line segments to intersect at one point other than endpoints are ((P1 P2 × P1 Q1 ) · P1 Z)((P1 P2 × P1 Q2 ) · P1 Z) < 0
and
((Q1 Q2 × Q1 P1 ) · Q1 Z )((Q1 Q2 × Q1 P2 ) · Q1 Z ) < 0.
(13)
Z'
Q1 Z
P2
P1
Q2
Fig. 8. Intersection checking between two line segments.
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
(a)
901
(b)
Fig. 9. Intersection occurrence and modification among triangle elements: (a) intersection occurring and (b) modification for the intersection case.
When a new triangle element is formed, the two edges connected with the new node should be checked whether they intersect with all existing triangle elements in the parametric space. If the intersections occur (see Fig. 9(a)), the two nodes should be merged to one new node, and the new node can be taken as the mid-point of the two nodes, as illustrated in Fig. 9(b). 5.2. Nodes merging or inserting The criterion of nodes merging or inserting is to check the distance between two adjacent nodes on the new generated front. If the distance is larger than a given upper limit, a new node is inserted; if the distance is less than a given lower limit, these two nodes should be merged to one. As the triangulation procedure proceeds in the parametric space and the distance between two adjacent nodes is calculated in the three dimensions, the metric map is introduced based on the first fundamental quantities of the surface. It gives an approximate distance map from the three-dimensional space to the parametric space. For a point p in the parametric space, P is the corresponding point on the surface, the tangent plane Tp of the surface at P can be considered as a Euclidean space with the metric map M [2,6,9], E F M= , (14) F G where E, F and G are the first fundamental quantities of the surface as defined in Eq. (5). In the neighborhood of point P, the difference between the distance with this metric map on Tp and the corresponding curve distance on the surface is only a high order small quantity and can be neglected. Since EG − F 2 > 0 [14], the map matrix is positive definite symmetric and it can be rewritten as 1 0 [e1 e2 ]T , M = [e1 e2 ] (15) 0 2 where e1 , e2 and 1 , 2 are the eigenvectors and eigenvalues of the matrix M, respectively. If h1 = h2 , the metric is isotropic and the directions are immaterial in such a case. Otherwise, the metric is anisotropic
902
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
so that e√ 1 and e2 will define the stretching directions of the elements in the parametric space. Given hi = 1/ i (i = 1, 2), hi and ei (i = 1, 2) are then computed to form an ellipse E with the center at C. It will be a unit circle in the tangent plane Tp if it is metric mapped into the three-dimensional space, as √ illustrated in Fig. 10. For a circle with radius r in Tp , the corresponding principal radii will be hi = r/ i (i = 1, 2) in the parametric space. Let r = 0.5 AB and searching ellipse E1 is formed, if the adjacent nodes are inside E1 , these nodes should be merged into one new node and the new node can be taken as the centroid of these nodes. Let r = 1.4 AB and searching ellipse E2 is formed, if the adjacent node is outside this searching ellipse, a new node should be inserted between these two nodes, and two corresponding triangle elements are formed. As shown in Fig. 11, polygonal line 1234 is the current front and 5, 6 and 7 are the new generated nodes, E1 and E2 are the searching ellipses with the center at node 6. Fig. 11 shows that node 5 is outside the ellipse E2 and node 7 is inside the ellipse E1 . This indicates that the three-dimensional distance between node 5 and 6 on the surface is larger than the upper limit and the distance between node 7 and 6 on the surface is less than the lower limit. Therefore, a e2 h2 h1
e1
C
A
B
Fig. 10. Searching ellipse in the parametric space.
5
E2 E1
8 1 6
2
9 7
3
4
Fig. 11. Nodes merging or inserting.
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
903
new node 8 is inserted between node 5 and 6 and two corresponding triangle elements 258 and 268 are formed. The new node 8 can be taken as the middle point of the node 5 and 6. Node 6 and 7 should be merged into one node, i.e. node 9 in Fig. 11. Node 9 can also be taken as the middle point of node 6 and 7. All of the formed triangle elements in the parametric space are illustrated by the dashed lines in Fig. 11.
6. Convergence checking Convergence checking is crucial for surface triangulation. In this study, a new approach is proposed by the authors for checking triangulation convergence. The approach is based on examining whether or not a point is inside a quadrangle. Given a simple quadrangle p1 p2 p3 p4 and a point p, as illustrated in Fig. 12, a horizontal half line is made from the point p in order to identify whether the point p is inside the quadrangle p1 p2 p3 p4 . The number of intersection points between the half line and the quadrangle p1 p2 p3 p4 is calculated. If the number is odd, p is inside the quadrangle; otherwise p is outside the quadrangle [16]. The convergence checking approach of the triangulation procedure is discussed in detail as follows. As illustrated in Fig. 13, polygonal line 12345 and 1 2 3 4 are parts of the current front, and 6, 7, 8, 9, 5 , 6 , 7 are new generated nodes. The two adjacent triangle elements of the current front on one side are joined to form a quadrangle. As shown in Fig. 13, the left triangle elements can form 3 quadrangles, i.e., 1276, 2387 and 3498. Then checking whether or not the new generated nodes of the front on the other side, i.e. nodes 5 , 6 and 7 in Fig. 13, are inside these quadrangles. If the node is inside the quadrangle, the vertex of the quadrangle which is closest (in the three-dimensional space) to the node can be determined. There may be two different situations. One is that the vertex is located on the current front, as illustrated in Fig. 13, node 7 is inside the quadrangle 3498 and closest to vertex 4, then new node 7 is obtained, which is taken as the mid-point of node 7 and 4. In addition, node 3 and 8, node 4 and 9 should be treated in the same way and nodes 3 , 4 are obtained. All triangle elements associated with these nodes, i.e. 4, 7 , 3 , 8, 4 and 9 should be adjusted accordingly, as shown in Fig. 14. The other case is that the closest vertex is one of the new generated nodes, as node 6 in Fig. 13, it is inside the quadrangle 2387 and is closest to vertex 7. This situation is relatively simpler than the former case. A new node 6 is determined by the mean p1
p p4 p2 p3 Fig. 12. Examining whether or not a point is inside a quadrangle.
904
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
1 6 5'
2
1'
7 3
4
2'
6' 8
7'
5
3'
4'
9 Fig. 13. Convergence checking approach.
1 6 2
5'
1'
7 6'
3
7 '' 5
2' 3 ''
4 '' Fig. 14. The closest vertex is on current front.
value of node 6 and 7 to replace these two nodes and the corresponding triangle elements are adjusted, and the final situation is depicted in Fig. 15. After all of the nodes have been checked, the one-side front which has formed quadrangles is fixed and the other-side front is refreshed while the nodes which are inside the quadrangles should be removed from the current front, such as node 6 and 7 . The procedure discussed above is repeated until the front is null, which indicates that the triangulation procedure is finished. Two problems may arise during the above procedure, and they must be dealt with properly. One is that triangle element with zero area may be generated during the convergence checking process. This kind of element should be discarded. The other is that some nodes on the one-side front which is fixed may be omitted, as node 7 in Fig. 16(a). The adjacent nodes of 7 are 6 and 8 and the mid-point of 6 and 8 is calculated. Then node 7 is replaced by the new node 7 , which is determined by the mean value of node 7 and the above mid-point. Two new triangle elements are formed accordingly, i.e. 67 2 and 7 82 in Fig. 16(b).
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
905
1 6 5'
2
1'
6 ''
3
2' 3 ''
7 ''
4 ''
5
Fig. 15. The closest vertex is one of the new generated nodes.
1 6 2
1'
7 3
2'
8
4
3'
5
4'
(a)
2
6 7'
3
8
2'
(b)
Fig. 16. Node is omitted when the triangulation is finished: (a) node 7 is ommitted and (b) two new triangle elements are formed.
7. Examples As an example, triangulation is implemented for a centrifugal compressor impeller. The impeller consists of 12 blades and one hub, as shown in Fig. 17. The blade is a three-dimensional free-form surface and the hub is a revolution surface. To assess the mesh quality, quality evaluation indicators are introduced. For a triangle ABC, quality is used to evaluate the individual element, 1 and 2 quality are used to assess the total triangle elements [12]. √
=2 3 1 =
CA 2
N T
i=1 i
NT
,
CA × CB , + AB 2 + BC 2 NT
2 = N T
i=1 1/i
.
(16) (17)
In the ideal case, = 1 and = 0 if A, B and C are co-linear. The hub is first triangulated. Chosen two circular boundaries as the initial fronts, the triangulation is advanced to the interior from these two boundaries. The results are shown in Fig. 18. Total 2350 triangles are created, 1 = 0.9885 and 2 = 0.9802. For this revolution surface, its normalized parametric area
906
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
Fig. 17. Centrifugal compressor impeller.
Fig. 18. Triangulation results of the hub.
is a unit square with four boundaries, u ∈ [0, 1], v ∈ [0, 1]. If triangulation was advanced from these four boundaries according to the conventional AFM, a redundant boundary line would be generalized which is actually nonexistent for this revolution surface. Triangulation will be much harder because of the boundary and the corner angles brought by this nonexistent boundary. However this drawback can be overcome by advancing triangulation from the real two circular boundaries.
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
907
The blade is a free-form surface and is modeled by B-spline, as shown in Fig. 19. Boundary 1 and 3 are u-parametric curves corresponding to u ∈ [0, 1], v = 0 and u ∈ [0, 1], v = 1, respectively, 2 and 4 are v-parametric curves corresponding to u = 1, v ∈ [0, 1] and u = 0, v ∈ [0, 1], respectively. The curve lengths of boundaries 1, 2, 3, and 4 are 47.1051, 329.9315, 119.3394 and 227.6364 mm, respectively by integrating Eqs. (7) and (8). The nodal spacing function Ed (x, y, z) is set as a constant of 5 mm and the actual nodal spacing distribution is 5.2339, 4.9990, 4.9725 and 4.9486 mm, respectively. The boundaries discretization results are shown in Fig. 20. Fig. 21 shows the triangulation results of the blade surface. The total number of triangle elements in Fig. 21 is 1702, 1 = 0.9942 and 2 = 0.9869. Majority of triangle elements with poor quality exist near the boundary 1, which is determined by the surface geometry property. For this blade surface, the corner angles formed by boundaries 1 and 2, 1 and 4 are both some 90◦ ; the angle formed by boundaries 3 and 4 is some 120◦ ; and the last angle formed by 2 and 3 is about 60◦ . It is obvious that the quality of the triangle elements near the boundary 1 is poor. The computed values of 1 and 2 are 0.8217 and 0.8125 and the minimum value of is 0.7324 in this area. Quality of the triangle elements near boundaries with poor corner conditions can be improved through subdivisions of the triangle elements. 2
1 4
Fig. 19. The centrifugal compressor blade.
Fig. 20. The boundaries discretization results of the blade.
Fig. 21. Triangulation results of the blade.
3
908
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
1
2
3
Fig. 22. The triangle status in the parametric space for trimmed surface.
Fig. 23. Triangulation results of the trimmed blade.
Triangulation of trimmed surface is also important in many branches of CAD/CAM and Graphics. A trimmed parametric surface consists of a parametric surface and a set of properly ordered trimming curves lying within the parametric rectangle of the surface. There are three statuses of triangles in the parametric space under the procedure of surface triangulation, i.e. IN, OUT and ON [17]. As shown in Fig. 22, the circle is the trimming curve in the parametric rectangle. Status IN means the triangle is in the valid parametric region, as triangle 1 in Fig. 22; status OUT means the triangle is in the trimmed area completely, as triangle 2; and status ON (triangle 3) means that the triangle intersects with the trimming curve. The triangles with status IN are reserved and the triangles with status OUT are discarded. The triangles with status ON will be trimmed so that the parts inside the valid parametric region are reserved and the others are discarded. The reserved parts of these triangles will be merged or divided according to their states respectively [18]. In order to simplify the computation, the trimming curve may be approximated by the polygonal line. Fig. 23 shows the triangulation results of a trimmed blade with the trimming curves consisting of a circle and a freeform curve in the parametric region. Total 1450 triangles are generalized and 1 = 0.9432, 2 = 0.9023. For the trimmed surface, the triangulation results of the whole initial surface can be regarded as the initial mesh and trimming of the initial mesh can be looked upon the post treatment, therefore the whole procedure is fast and robust.
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
909
8. Conclusions and discussion A developed advancing front method is presented to triangulate the three-dimensional parametric surface. The details of implementation, validity and convergence checking are given in the paper. The method is advanced from two longer opposite boundaries and the number of triangle elements at surface corner is determined by the surface corner angle. This developed method solves the problem caused by poor surfaces corner condition in the conventional AFM and increases the flexibility on the corner and boundary treatment. The new node is generated directly in three dimensions and the validity checking of the new triangle element is performed in the parametric space, therefore the advantages of AFM in two dimensions are retained sufficiently. A simple but effective convergence checking approach is presented as well, which is based on examining whether a point is inside a quadrangle or not. Numerical experiments show that the developed method is robust and easy to carry out and high quality triangulation results have been achieved.
Acknowledgements The work reported in this paper was supported by the Research Fund for the Doctoral Program of Higher Education of China under Grant No. 20020698042, and the National Science Foundation of China under Grant No. 50146014 (2002). References [1] M. Vigo, N.P. Garciia, P.B. Crosa, Directional adaptive surface triangulation, Comput. Aided Geom. Design 16 (1999) 107–126. [2] C.K. Lee, Automatic adaptive mesh generation using metric advancing front approach, Eng. Comput. 16 (1999) 230–263. [3] J.C. Cuilliere, An adaptive method for the automatic triangulation of 3D parametric surfaces, Comput. Aided Design 30 (2) (1998) 139–149. [4] P.L. Baehmann, S.L. Wittchen, M.S. Shephard, Robust, geometrically based, automatic two-dimensional mesh generation, Internat. J. Numer. Methods Engng. 24 (1987) 1043–1078. [5] T.P. Fang, L. Piegl, Delaunay triangulation using a uniform grid, IEEE Comput. Graph. Appl. 13 (1993) 36–47. [6] H. Borouchaki, P.L. George, F. Hecht, Delaunay mesh generation governed by metric specifications, Part I, Algorithms, Finite Elements Anal. Des. 25 (1997) 61–83. [7] M. Vigo, N. Pla, Computing directional constrained Delaunay triangulations, Comput. Graph. 24 (2000) 181–190. [8] P.L. George, E. Seveno, The advancing front mesh generation method revisited, Internat. J. Numer. Methods Engng. 37 (1994) 3605–3619. [9] J.R. Tristano, S.J. Owen, S.A. Canann, Advancing front surface mesh generation in parametric space using a Riemannian surface definition, in: Proceedings of the Seventh International Meshing Roundtable, Dearborn MI, USA, 1998. [10] T.S. Lau, S.H. Lo, Finite element mesh generation over analytical curved surfaces, Comput. & Structures 59 (2) (1996) 301–309. [11] V. Francois, J.C. Cuilliere,Automatic mesh pre-optimization based on the geometric discretization error,Adv. Eng. Software 31 (2000) 763–774. [12] Z.Q. Guan, C. Song, Y.X. Gu, Recent advances of research on finite element mesh generation methods, J. Comput. Aided Des. Comput. Graph. 15 (1) (2003) 1–14. [13] J.C. Cuilliere, A direct method for the automatic discretization of 3D parametric curves, Comput. Aided Design 29 (1997) 639–647.
910
B. Wu, S. Wang / Finite Elements in Analysis and Design 41 (2005) 892 – 910
[14] D.J. Meng, K. Liang, Differential Geometric, Science Press, Beijing, China, 2002. [15] Y.H. Li, Y.X. Hua, A method to quickly accessing the intersection of line segments, Bull. Surve. Mapping Eng. (2003) 30–31. [16] P.D. Zhou, Computational Geometric—Algorithms Analysis and Design, Tsinghua University Press, Beijing, 2000. [17] L.A. Piegl, W. Tiller, Geometry-based triangulation of trimmed NURBS surfaces, Comput. Aided Design 30 (1) (1998) 11–18. [18] M.W. Ji, C.G.Yang, J.G. Sun, A fast algorithm for triangulating the parametric space of trimmed NURBS surfaces, Chinese J. Comput. 19 (6) (1996) 450–456.