- Email: [email protected]
Improvement of quadrilateral meshes for discretization of tunnels
Computers and Geotechnics 31 (2004) 47–56 www.elsevier.com/locate/compgeo
Improvement of quadrilateral meshes for discretization of tunnels A. Menendez Dıaz
a,*
, C. Gonz alez Nicieza b, A.E. Alvarez Vigil
c
a
b
Department of Construction Engineering and Manufacturing, Mining Engineering School, University of Oviedo, Independencia 13, 33004 Asturias, Spain Department of Mining Engineering, Mining Engineering School, University of Oviedo, Independencia 13, 33004 Asturias, Spain c Department of Mathematics, Mining Engineering School, University of Oviedo, Independencia 13, 33004 Asturias, Spain Received 7 February 2003; received in revised form 12 October 2003; accepted 16 October 2003
Abstract Automatic mesh generators of 2D polygonal domains, especially those that use structured meshes, generate elements (triangles or quadrilaterals), which, in relatively complex models, tend to lose their equilateral form. In order to avoid this, the posterior smoothing of elements is necessary to obtain meshes with minimum geometrical distortions. In this paper, one of these smoothing techniques is developed for use in underground modelling. A system of equations is generated in a similar mode to structural analysis, where the degrees of freedom are the smoothed values of each node coordinate, and the local stiffness matrix is a measure of the geometrical connectivity of each element in the mesh. This facilitates the smoothing of non-convex domains with interior islands, which allows the imposition of directionally constrained nodes on the boundary of the domain. At present, this method is being utilized to compliment our mesh generation technique for the discretization of tunnels in 2D. Ó 2003 Elsevier Ltd. All rights reserved.
1. Introduction The generation of quadrilateral meshes in finite difference or finite element programs, constitutes a commonly used tool in the modelling of physical, chemical or purely mathematical phenomena. In our program of analysis, Almec2D, there are certain characteristics to be controlled in the mesh: the local density of points, the smoothness of the point, and the shape of the elements. High densities give more accuracy, but the computation takes longer. This leads to adaptive meshing methods. Large variations in grid density or shape can cause numerical problems [3]. As mesh generation becomes more automated, and as mesh sizes increase, the need to create meshes completely free of badly distorted elements also increases. Irregular elements of this kind produces inaccurate results or instability in the posterior analysis process. To mitigate this problem, it is common to improve the quality of the initial mesh before carrying out the *
Corresponding author. Tel.: +33-98-5104266; fax: +34-98-5104265. E-mail address: [email protected] (A. Menendez Dıaz). 0266-352X/$ - see front matter Ó 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2003.10.001
analysis. The smoothing methods analyse the relocation, addition or elimination of nodes in the mesh in order to minimise the geometrical distortions. In Almec2D, a structured meshing technique is used to obtain the initial quadrilateral mesh (see Fig. 1). First, the geometric model is divided into convex regions delimited by four boundary polylines (see Fig. 1(a)), later, each region is subdivided into quadrilaterals (quads), using Coons patches [5] to interpolate inner nodes (see Fig. 1(b)); and finally the boundary nodes are reconnected in readiness for the readjustment of the topology of the global mesh. The generated quads are very regular in the inner regions, but the shape of those near the boundaries is determined by the geometric constraints of the boundary. Thus, smoothing techniques are necessary to readjust the mesh to reasonable levels of geometric distortion near these boundaries. After making a brief description of the studies made in this field, a new optimization smoothing technique for quadrilateral meshes is developed in this paper, showing its field of application, the form in which it is implemented numerically, and illustrating the results in the smoothing of a series of particular meshes. The quality
48
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
(a)
(b)
(c)
Fig. 1. The parametric mesh of a vertical tunnel cross section: (a) boundary polylines defining the regions; (b) meshed regions; (c) reconnected region of the initial mesh.
of the resulting meshes is evaluated by defining a metric estimator of the degree of uniformity, and fundamental conclusions are drawn. All of the above is being applied in order to improve the structured meshes obtained with the analysis program Almec [12], and although this paper considers two dimensions, the concepts can be generalised to three dimensions.
2. Previous work and bibliographical review The smoothing techniques can be grouped on the basis of how they affect the mesh. On the one hand, there are those that require the insertion of new nodes [15], whilst others are based on the interchange of edges [4] or in the identification and elimination of distorted elements [2]. Lastly, are the methods that modify the location of part of the nodes without changing the connectivity of the mesh, but improving substantially the degree of regularity [6,9]. For these, smoothing Laplacian or optimization-based techniques can be used. The techniques based on Laplacian smoothing relocate the nodes recursively, based on the relative position of each node with respect to its adjacent nodes. This technique works quite well for meshes in convex domains, but it does not work so well in domains with multiple concavities, or that present interior islands. Iterative processes are relatively slow because the distortions of elements propagate to subsequent ones, requiring the establishment of a minimum number of iterations to obtain a reasonable result. In order to make the process faster, it is possible to introduce weights to each neighbouring node, on the basis of geometric factors, or to establish restrictions in the movement of the nodes [9]. Whichever option is taken, the method of moving nodes on the basis of their proximity to their neighbours, continues to be a heuristic procedure. In contrast, the methods based on optimization-based smoothing [7,13] establish a function that measures the global distortion of the mesh, and they determine the new coordinates of the nodes by minimising this function. Optimizationbased smoothing techniques are more expensive than the Laplacian ones, but they give better results in concave zones. Mixed methods are used for very big meshes,
synchronizing Laplacian and Optimization-based smoothing techniques, using the latter locally [3]. Finally, smoothing techniques which try to take advantage of physical phenomena have also been used. Lohner [11] developed a method that considers the edges of the mesh composed of an assembly of springs that join the nodes between them, and over which forces that tend to obtain a more regular form in each element are applied. Shimada [14] viewed the nodes as being the centre of an ideal bubble, which is relocated on the basis of its contact with other bubbles until equilibrium is achieved.
3. Method for quadrilateral mesh smoothing We have developed an optimization smoothing technique based on a physical analogy, where each element of the mesh is considered a deformable material, such that the smoothed coordinates are simply the degrees of freedom. Therefore, an elementary geometrical stiffness matrix for each quadrilateral has been defined, and by assembling the elements we obtain a matrix of the global stiffness of the system, by which the location of the new nodes of the mesh are calculated. The optimization-based smoothing technique is implemented in the following stages: 1. Estimation of the degree of uniformity of the initial mesh. 2. Calculation of the element stiffness matrix for each quad. 3. Assembly of the matrix of global stiffness of the mesh. 4. Imposition of the boundary conditions. 5. Computation of MckeeÕs method for the solution of the equationÕs system. 6. Estimation of the degree of uniformity of the smoothed mesh. In addition, it is possible to modify the global stiffness matrix in order to consider different boundary conditions to restrict the movements of the nodes. Below we develop the system of equations for the re-solution of the problem. 3.1. Element stiffness matrix of a quad The main idea of the smoothing technique is to consider the mesh as a deformable solid whose stiffness
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
depends on the degree of geometric distortion of each element. For its determination, the methodology developed by Balendran is followed [1], but we introduce a shape factor. For a quadrilateral mesh, let a given quad have vertices Pi Pj Pk Pl (see Fig. 2(a)). The ideal local relationship between the coordinates of the vertices is when the two diagonals are perpendicular to each other and, in addition, intersect at the midpoint. Since this will rarely happen (see Fig. 2(b)), we will have to take these nodes to positions Pi0 Pj0 Pk0 Pl0 (smoothed points) so that the ideal conditions are fulfilled (see Fig. 2(c)). If we consider one of the diagonals as fixed, for example Pj Pl (see Fig. 3(a)), points Pi and Pk are moved until they lie on the median line of Pj Pl . The distance of the new smoothed points Pi0 and Pk0 to the midpoint of the diagonal can vary, depending on whether we want the improved quads to tend towards rhombic or square forms (this is controlled by a shape factor l, the ratio between the magnitudes of the segments Pi0 Pk0 and Pj Pl ). l¼
49
(b)
(a)
Fig. 3. Smoothed element by deforming the initial quad.
x0k yk0
¼
1 xj þ xl 1 ðyl yj Þ l xl xj 2 yj þ yl 2
ð8Þ
and then xj þ lyj þ 2x0k xl lyl ¼ 0;
ð9Þ
lxj yj þ 2yk0 þ lxl yl ¼ 0:
ð10Þ
Once the mesh is smoothed, we will have
Pi0 Pk0 k~ P0 ~ P 0k k ¼ i : Pj Pl k~ Pj ~ P lk
ð1Þ
From the metric properties of Fig. 3(a) in node Pi0 it can be concluded that: ~ P 0i ¼ ~ P m þ 12 k~ P 0i ~ P 0k k ~ u;
ð2Þ
where Pm is the mid point of diagonal Pj Pl , and ~ u is the unit vector perpendicular to this diagonal. 1 ðyl yj Þ ~ u¼ ; ð3Þ xl xj P lk k~ Pj ~ x and y being the values of the coordinates of the points indicated in each subscript. If we express each component of Eq. (2) we have: 0 1 xj þ xl 1 ðyl yj Þ xi ¼ þ l ð4Þ xl xj yi0 2 yj þ yl 2 and then 2x0i xj lyj xl þ lyl ¼ 0;
ð5Þ
2yi0
ð6Þ
þ lxj yj lxl yl ¼ 0:
For node Pk0 we have ~ P 0k ¼ ~ P m 12 k~ P 0i ~ P 0k k ~ u;
ð7Þ
(a)
x0t xt ¼ 0 yt0 yt ¼ 0
t ¼ i; j; k; l
ð11Þ
enabling us to deduce from the previous relations the first, second, fifth and sixth rows of the matrix equation: 3 2 03 2 xi 2 0 1 l 0 0 1 l 6 07 7 6 0 2 l 1 0 0 l 1 7 6 yi 7 6 7 6 07 6 6 1 k xj 7 2 0 1 k 0 0 7 6 7 7 6 6 6 07 7 6 k 1 0 6 y 2 k 1 0 0 7 6 j7 6 7 7 6 6 0 x0k 7 0 1 l 2 0 1 l 7 6 7 7 6 6 6 07 6 0 yk 7 0 l 1 0 2 l 1 7 7 6 6 7 6 07 6 7 4 1 k 0 0 1 k 2 0 5 6 4 xl 5 1 2 3 0 607 6 7 6 7 607 6 7 607 6 7 ¼ 6 7: 607 6 7 607 6 7 6 7 405
k
0
(b) Fig. 2. Smoothing of a quad.
(c)
0
0
k
1
0
2
yl0
ð12Þ
50
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
The rest of equations can be determined using the other diagonal Pi Pk as fixed (see Fig. 3(b)) where k¼
Pj0 Pl0
ð13Þ
:
Pi Pk
Considering the symmetry of the problem k ffi l. When l ¼ 1 we will have dependency relationships that will transform a quad into a square, while for the rest of the cases the smoothed quad will tend to have rhombic forms. In this study l is approximated by a geometric constant, and it depend only of the initial geometry of the quad, as will be justified in Section 5. The elementary stiffness matrix k e ½8 8 of a generic quad is determined by the previous matrix relationship. 3 2 2 0 1 l 0 0 1 l 6 0 2 l 1 0 0 l 1 7 7 6 6 1 l 2 0 1 l 0 0 7 7 6 6 l 1 0 2 l 1 0 0 7 e 7: 6 k ¼6 0 1 l 2 0 1 l 7 7 6 0 6 0 0 l 1 0 2 l 1 7 7 6 4 1 l 0 0 1 l 2 0 5 l 1 0 0 l 1 0 2 ð14Þ This symmetric matrix can be expressed in terms of nodal sub-matrices as 2 3 kii kij kik kil 6 kji kjj kjk kjl 7 7 ke ¼ 6 ð15Þ 4 kki kkj kkk kkl 5: kli kjl klk kll
kie being the matrix of elementary stiffness of the quad i properly referred to its degrees of freedom considering the global identification of the nodes. The degrees of freedom are the smoothing coordinates, and they have been put in the order: ½x01 y10 x02 y20 x0m ym0 ;
ð20Þ
where m is the number of nodes. The assembled matrix is presented in Fig. 4. We have assumed that the numbering sequence of the nodes is i < k < l < j, in order to indicate the global assembling of the matrix elements. Actually, when assembling the matrix of global stiffness, we consider how each element affects the surrounding elements, thereby relating the degrees of freedom of each quad with its adjacent ones. Intuitively, it can be observed that the greater the number of quads sharing a particular node, the greater the stiffness of that node and, therefore, the greater the weight in the determination of the smoothed coordinates of the adjacent nodes. Once the global stiffness matrix of the system is assembled, we must solve the sparse system of linear algebraic equations defined by K~ c0 ¼ ~ b;
ð21Þ T ½x01 y10 x02 y20 x0m ym0
where ~ c0 ¼ the smoothed coordinates, right hand side, and K the dimension 2m 2m.
is the solution vector with T ~ b ¼ ½b1 b2 b3 b4 b2m the global stiffness matrix of
3.3. Assigning boundary conditions
Letting 2 D¼ 0
0 2
O¼
0 0 0 0
L¼
1 l
l ; 1
ð16Þ
The system of equation (Eq. (21)) would be indeterminate unless some of the nodes of the mesh were fixed.
we have kii ¼ kjj ¼ kkk ¼ kll ¼ D
kij ¼ kjk ¼ kkl ¼ kli ¼ L;
kik ¼ kki ¼ kjl ¼ klj ¼ O
kji ¼ kkj ¼ klk ¼ kil ¼ LT ;
x'i y'i
x'k y'k
x'l y'l x' j y' j
ð17Þ T
where L is the transpose matrix of L. This symmetric matrix k e can be expressed in terms of the nodal submatrices in the form: 2 3 D L 0 LT 6 LT D L 0 7 7: ð18Þ ke ¼ 6 40 LT D L 5 T L 0 L D
x'i y'i
ki i
ki k
ki l ki j
x'k y'k
kk i
kkk
kk l kk j
3.2. Assembly of the global stiffness matrix
x'l y' l x'j y'j
kl i
kl k
kl l kl j
kj i
k jk
kj l kj j
The matrix of global stiffness K for n quads will be assembled as X K¼ kie ; ð19Þ i¼1;n
Fig. 4. Assembling of the matrix of global stiffness K.
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
The fixed conditions are applied to the nodes of the boundary of the zone to be smoothed, or to any fixed inner nodes. If a node is free (non-fixed), the corresponding degrees of freedom in vector ~ b will be zero. Specific values will only appear if the node is fixed. The value of the latter can be physically interpreted as a force that prevents the point from moving due to the initial geometric distortion of the mesh. In order to handle the fixed value c0fix in the solution of the system of equations we proceed as in the methods of structural analysis; the stiffness in the main diagonal K½n; n of degree of freedom n is replaced by a very high fictitious stiffness (k1 1020 ), and simultaneously this fictitious stiffness multiplied by the constant value is introduced in the RHS of the system of equations (see Eq. (22)). b½n ¼ k1 c0fix :
K½n; n ¼ k1
ð22Þ
A fictitious stiffness of this value can be considered for numerical effects as infinity. When the system of equations is solved, the outcome of the value of the smoothed coordinate for that point is the same as the initial fixed value c0 ½n ¼ c0fix . In this way, the x or y coordinates of a node can be fixed since each degree of freedom is independent. In the case of the outer boundaries and inner boundaries of islands or isolated quads and segregated zones, the normal procedure is to fix both degrees of freedom of the node. A great advantage of this smoothing technique, is that boundary conditions are introduced only by altering the global stiffness of the nodes involved. This method is not limited to boundaries parallel to the x and y axes (see Fig. 7). If there is a fixed lineal dependency between the degrees of freedom n and g, given by the ratio rfix and the constant value sfix , see Eq. (23), we input this relationship using the fictitious stiffness given by Eq. (24) for K and Eq. (25) for ~ b. c0 ½n þ rfix c0 ½g ¼ sfix ; K½n; n ¼ k1 K½g; g ¼
2 rfix
ð23Þ
K½n; g ¼ K½g; n ¼ rfix k1 k1 ;
b½n ¼ sfix k1
b½g ¼ rfix sfix k1 :
ð24Þ ð25Þ
3.4. Resolution of the system of equations In Eq. (21) most of the terms of the matrix of global stiffness K are zero, since those degrees of freedom whose nodes do not share the same quad, are not bound by any stiffness, and therefore the input in the matrix of stiffness corresponding to those degrees of freedom is zero. A judicious reordering of the rows and columns for the coefficients matrix can lead to enormous savings in computer execution time and storage. Thus, we have
51
proposed a permutation that allows the rearrangement of the degrees of freedom, such that the same matrix can be stored in the most economic form. The system considered to implement the permutation has been the Mckee algorithm, while for the storage of the final matrix, the band and envelope method has been used [8]. Once the matrix K is stored, the system of equations is solved through the conventional Cholesky method, restoring the original equation order later, with the purpose of reassigning the values of smoothed coordinates to each node. If the matrix K is not symmetric and positive definite, some form of pivoting (row and/or column interchanges) is necessary to ensure the numerical stability of CholeskyÕs method. In order to avoid such difficulties, we have assumed k ¼ l and 0 < l 6 1; in this conditions K is symmetric and positive definite. 3.5. Two element mesh example In order to illustrate the method used in this work, we have considered the mesh of two quads of the Fig. 5. Assembly of the stiffness matrices for the two quads and using Eq. (22), leads to the following 12 equations: 2
4
6 6 0 6 6 2 6 6 6 0 6 6 0 6 6 6 0 6 6 1 6 6 6 l 6 6 1 6 6 6 l 6 6 0 4 0
0 4 0 2 0 0 l 1 l 1 0
2
0
0
0
1
0
2
0
0
l 1
l
1 l
0
1
0
0
3 2
x01
3
7 6 7 0 7 6 y10 7 7 6 7 6 07 4 0 1 l 0 0 0 0 1 l 7 7 6 x2 7 7 6 07 0 4 l 1 0 0 0 0 l 1 7 6 y2 7 7 6 7 6 07 1 l k1 0 1 l 0 0 0 0 7 7 6 x3 7 7 6 07 l 1 0 k1 l 1 0 0 0 0 7 6 y3 7 76 7 6 07 0 0 1 l k1 0 0 0 0 0 7 7 6 x4 7 7 6 07 6 7 0 0 l 1 0 k1 0 0 0 0 7 7 6 y4 7 6 07 0 0 0 0 0 0 k1 0 1 l 7 7 6 x5 7 7 6 07 6 7 0 0 0 0 0 0 0 k1 l 1 7 7 6 y5 7 7 0 7 1 l 0 0 0 0 1 l k1 0 5 6 4 x6 5 y60 l 1 0 0 0 0 l 1 0 k1 3
0 2 0 7 6 60 7 7 6 60 7 7 6 7 6 60 7 7 6 6 0k 7 6 17 7 6 6 1k1 7 7 ¼6 6 0k 7: 6 17 7 6 6 0k1 7 7 6 6 2k 7 6 17 7 6 6 0k1 7 7 6 6 2k 7 4 15 1k1
l
ð26Þ
Setting l ¼ 1 and k1 ¼ 1020 and solving the equations gives T ~ c0 ¼ ½0; 1; 1; 1; 0; 1; 0; 0; 2; 0; 2; 1 :
ð27Þ
52
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
Fig. 5. Initial mesh and smoothed mesh of two quads.
3.6. Degree of uniformity of the mesh Before and after the mesh has been smoothed, its degree of uniformity is calculated. To do this, it is important to select a suitable metric. When the verification of the value of the degree of uniformity of an element is above a certain threshold, the corresponding nodes are not included, and therefore, the size of the matrix of global stiffness is reduced. In this paper, we have taken a as the degree of uniformity, for a triangle with nodes Pi Pj Pk , as proposed by Lee and Lo [10] (see Fig. 6(a)), that is ! ! pffiffiffi kPk Pi Pk Pj k a ¼ 2 3 ! ð28Þ ! 2 ! 2 : 2 kPk Pi k þ kPi Pj k þ kPj Pk k To determine the degree of uniformity b of a quad (see Fig. 6(b) and (c)), it is divided into the two families of triangles bordered by their diagonals. The degree of triangular uniformity for each one of the four triangles is sorted in descending order of magnitude (a1 P a2 P a3 P a4 ). Finally, the value of the uniformity of the quad is a3 a4 b¼ : ð29Þ a1 a2
4. Mesh quality improvement results for the smoothing technique The scheme of smoothing developed here has been implemented in C++ within the Almec2D program, and it has been used to treat different meshes generated by this program. In Figs. 7–10 several tunnel models are illustrated indicating on the left the initial mesh, and on the right of the final smoothed mesh. As can be observed, the improvement in the mesh is very significant, and it tries to introduce the morphology of the tunnel cross-section very gradually. The minimum
(a)
(b) Fig. 7. Initial mesh (a) and smoothed mesh (b) of model A.
This degree of uniformity is 6 1 and approaches unity as the quad approaches a square.
(a) Fig. 6. Degree of uniformity a of the meshÕs elements.
(b) Fig. 8. Initial mesh (a) and smoothed mesh (b) of model B.
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
(a)
53
(b) Fig. 9. Initial mesh (a) and smoothed mesh (b) of model C.
Fig. 11. Degree of uniformity before smoothing (model D).
(a)
(b) Fig. 10. Initial mesh (a) and smoothed mesh (b) of model D.
and average values of b for each model are collected in Table 1, before and after making the smoothing. In the last column, we point out the increase of b in percentage. An increase of the average degree of uniformity is seen in all the models. In the best case, it is of the order of 4.5% (model A), although improvement values from 1.5% to 3.0% seem more representative (models C and D), as these are of the same order as improvements obtained with other smoothing algorithms [3,10]. The average value of the degree of uniformity does not clearly illustrate the improvement that takes place in the mesh. In order to overcome this difficulty we have to present the degree of local uniformity of each quadrilateral by mean of a map of isovalues. The zones are represented by colours ranging from red for near zero b values to blue for values close to unity. This legend and the degree of uniformity before smoothing are represented in Fig. 11, and after smoothing in Fig. 12. The improvement in the mesh can be observed clearly especially in the quadrilaterals adjacent to the boundaries of the sub-meshed into which the initial model was divided (see Fig. 1(a)).
Fig. 12. Degree of uniformity after smoothing (model D).
The results obtained for the rest of models are illustrated in Figs. 13–15.
5. Analysis of the shape factor In all the previous examples we have used a shape factor l defined by Eq. (32). Nevertheless, to complete the work we have to justify this choice. To do this, different l definitions have been studied in order to see
Table 1 Degree of uniformity of the meshes Model
A B C D
Number of nodes
440 1089 1790 1114
Number of quads
400 1000 1688 1038
Db (%)
Degree of uniformity b Before smoothing
After smoothing
Min
Average
Min
Average
0.334 0.042 0.254 0.391
0.699 0.728 0.884 0.858
0.363 0.050 0.278 0.356
0.744 0.737 0.905 0.878
4.5 0.9 2.1 2.8
54
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
Fig. 13. Degree of uniformity of initial mesh (a) and smoothed mesh (b) of model A.
Fig. 14. Degree of uniformity of initial mesh (a) and smoothed mesh (b) of model B.
Fig. 15. Degree of uniformity of initial mesh (a) and smoothed mesh (b) of model C.
which one gives the greatest degree of uniformity. In particular, four cases have been studied, as shown in Table 2. In case 1, l is constant for all quadrilaterals, in case 2 the relationship between the minimum and maxi-
mum diagonals is considered, in case 3 the relationship between the minimum and maximum side, and finally case 4, l is taken equal to the value of the degree of uniformity of the quadrilateral b given by Eq. (27).
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56 Table 2 Degree of uniformity and shape factor for each case Case
Maximum b
Shape factor l
1
0.772
l¼1
ð30Þ
2
0.795
minðPj Pl ; Pi Pk Þ l¼ maxðPj Pl ; Pi Pk Þ
ð31Þ
3
0.799
4
0.793
l¼
minðPi Pj ; Pj Pk ; Pk Pl ; Pl Pi Þ maxðPi Pj ; Pj Pk ; Pk Pl ; Pl Pi Þ
l¼b
ð32Þ ð33Þ
55
minimizes the user time required to produce an acceptable mesh in a simple way. It is a direct method that does not need iterative techniques of calculation and whenever the boundary conditions are well defined, the degree of uniformity improves remarkably. The shape factor helps to increase the degree of uniformity, especially when Eq. (32) is taken for its definition. The disadvantage is the high cost of memory when working with a matrix of global stiffness, but even so, this cost is less than that needed to make the final finite element analysis. The use of MckeeÕs method to solve the system of equations has been specially useful to mitigate this difficulty. This method is used within Almec2D to obtain more regular meshes, and it is one of the techniques for improving the numerical results of the calculation module. The possibility of extending this method to hexahedral meshes in 3D for the Almec3D program is being analyzed.
Acknowledgements Fig. 16. Degree of uniformity of cases 1–2–3–4.
In order to study several meshes, the simplest mesh has been taken as a reference, (the one corresponding to model A), and the density of nodes has been increased in the radial direction to obtain a series of eight meshes, identified by their number of nodes: 120, 200, 280, 360, 440, 520, 600 and 680 nodes. In Fig. 16 the values obtained for the degree of uniformity for eight consecutive meshes are represented. In the second column of Table 2, the maximum values reached by the degree of uniformity b are indicated. Initially b increases proportionally with the increase in the number of nodes, converging on a value of 0.8. However, it can be observed that from 500 nodes onwards, no further improvement is obtained in the degree of uniformity, regardless of the value l. Case 1 behaves better for meshes with a smaller number of nodes, but cases 2, 3 and 4 give the best results for meshes of more than 450 nodes. The greatest increase of average uniformity degree is obtained for case 3 (see Table 2), and for this reason we have used the value of l given by Eq. (32).
6. Conclusions An efficient 2D smoothing method has been developed. The method works with quadrilateral meshes by considering a criterion of global optimization that
The authors gratefully acknowledge the support of Robin Walker for the preparation of this paper in English.
References [1] Balendran B. A direct smoothing method for surface meshes. In: Proceedings of the Eighth International Meshing Roundtable, South Lake Tahoe, CA, USA 1999. p. 189–93. [2] Blacker TD, Stephenson MB. Paving: a new approach to automated quadrilateral mesh generation. Int J Numer Meth Eng 1991;32:811–47. [3] Canann Scott A, Tristano JR, Staten ML. An approach to combined Laplacian and optimization-based smoothing for triangular, quadrilateral, and quad-dominant meshes. In: Seventh International Meshing Roundtable. Sandia National Labs; 1998. p. 479–494. [4] Edelsbrunner H, Shah N. Incremental topological flipping works for regular triangulations. In: Proceedings of the Eighth ACM Symposium on Computational Geometry, 1992. p. 43–52. [5] Faux ID, Pratt MJ. Computational geometry for design and manufacture. Chichester, West Sussex, England: John Wiley & Sons; 1987. [6] Field D. Laplacian smoothing and Delaunay triangulations. Commun Numer Meth Eng 1988;4:709–12. [7] Freitag L. On combining Laplacian and optimization-based mesh smoothing techniques AMD trends in unstructured mesh generation. In: ASME, vol. 220. 1997. p. 37–43. [8] George A, Liu J. Computer solution of large sparse positive definite systems. New Jersey: Prentice-Hall; 1981. [9] George PL, Borouchaki H. Delaunay triangulation and meshing. Application to finite elements. Paris, France: Editions Hermes; 1998. p. 230–4.
56
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
[10] Lee CK, Lo SH. A new scheme for the generation of a graded quadrilateral mesh. Comput Struct 1994;52(5):847–57. [11] Lohner R, Morgan K, Zienkiewicz OC. Adaptive grid refinement for the compressible Euler equations. In: Accuracy estimates and adaptive refinements in finite element computations. New York: Wiley; 1986. p. 281–97. [12] Menendez Dıaz A, Gonzalez Nicieza C, Alvarez Vigil A, Morıs Menendez G. Malladores Geometricos para el Disesimno Estructural por Diferencias Finitas, Disegno di Macchine e Progettazione Industriale, Junio 1996. p. 19–32.
[13] Parthasarathy V, Kodiyalam S. A constrained optimization approach to finite element mesh smoothing. Finite Elements Anal Des 1991;9:309–20. [14] Shimada K. Physically-based mesh generation: Automated triangulation of surfaces and volumes via bubble packing, Tesis Doctoral, ME Department, Massachusetts Institute of Technology, Cambridge, MA, 1993. [15] Staten ML, Canann SA. Post refinement element shape improvement for quadrilateral meshes, trends in unstructured mesh generation. ASME 1997;220:9–16.
Improvement of quadrilateral meshes for discretization of tunnels A. Menendez Dıaz
a,*
, C. Gonz alez Nicieza b, A.E. Alvarez Vigil
c
a
b
Department of Construction Engineering and Manufacturing, Mining Engineering School, University of Oviedo, Independencia 13, 33004 Asturias, Spain Department of Mining Engineering, Mining Engineering School, University of Oviedo, Independencia 13, 33004 Asturias, Spain c Department of Mathematics, Mining Engineering School, University of Oviedo, Independencia 13, 33004 Asturias, Spain Received 7 February 2003; received in revised form 12 October 2003; accepted 16 October 2003
Abstract Automatic mesh generators of 2D polygonal domains, especially those that use structured meshes, generate elements (triangles or quadrilaterals), which, in relatively complex models, tend to lose their equilateral form. In order to avoid this, the posterior smoothing of elements is necessary to obtain meshes with minimum geometrical distortions. In this paper, one of these smoothing techniques is developed for use in underground modelling. A system of equations is generated in a similar mode to structural analysis, where the degrees of freedom are the smoothed values of each node coordinate, and the local stiffness matrix is a measure of the geometrical connectivity of each element in the mesh. This facilitates the smoothing of non-convex domains with interior islands, which allows the imposition of directionally constrained nodes on the boundary of the domain. At present, this method is being utilized to compliment our mesh generation technique for the discretization of tunnels in 2D. Ó 2003 Elsevier Ltd. All rights reserved.
1. Introduction The generation of quadrilateral meshes in finite difference or finite element programs, constitutes a commonly used tool in the modelling of physical, chemical or purely mathematical phenomena. In our program of analysis, Almec2D, there are certain characteristics to be controlled in the mesh: the local density of points, the smoothness of the point, and the shape of the elements. High densities give more accuracy, but the computation takes longer. This leads to adaptive meshing methods. Large variations in grid density or shape can cause numerical problems [3]. As mesh generation becomes more automated, and as mesh sizes increase, the need to create meshes completely free of badly distorted elements also increases. Irregular elements of this kind produces inaccurate results or instability in the posterior analysis process. To mitigate this problem, it is common to improve the quality of the initial mesh before carrying out the *
Corresponding author. Tel.: +33-98-5104266; fax: +34-98-5104265. E-mail address: [email protected] (A. Menendez Dıaz). 0266-352X/$ - see front matter Ó 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2003.10.001
analysis. The smoothing methods analyse the relocation, addition or elimination of nodes in the mesh in order to minimise the geometrical distortions. In Almec2D, a structured meshing technique is used to obtain the initial quadrilateral mesh (see Fig. 1). First, the geometric model is divided into convex regions delimited by four boundary polylines (see Fig. 1(a)), later, each region is subdivided into quadrilaterals (quads), using Coons patches [5] to interpolate inner nodes (see Fig. 1(b)); and finally the boundary nodes are reconnected in readiness for the readjustment of the topology of the global mesh. The generated quads are very regular in the inner regions, but the shape of those near the boundaries is determined by the geometric constraints of the boundary. Thus, smoothing techniques are necessary to readjust the mesh to reasonable levels of geometric distortion near these boundaries. After making a brief description of the studies made in this field, a new optimization smoothing technique for quadrilateral meshes is developed in this paper, showing its field of application, the form in which it is implemented numerically, and illustrating the results in the smoothing of a series of particular meshes. The quality
48
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
(a)
(b)
(c)
Fig. 1. The parametric mesh of a vertical tunnel cross section: (a) boundary polylines defining the regions; (b) meshed regions; (c) reconnected region of the initial mesh.
of the resulting meshes is evaluated by defining a metric estimator of the degree of uniformity, and fundamental conclusions are drawn. All of the above is being applied in order to improve the structured meshes obtained with the analysis program Almec [12], and although this paper considers two dimensions, the concepts can be generalised to three dimensions.
2. Previous work and bibliographical review The smoothing techniques can be grouped on the basis of how they affect the mesh. On the one hand, there are those that require the insertion of new nodes [15], whilst others are based on the interchange of edges [4] or in the identification and elimination of distorted elements [2]. Lastly, are the methods that modify the location of part of the nodes without changing the connectivity of the mesh, but improving substantially the degree of regularity [6,9]. For these, smoothing Laplacian or optimization-based techniques can be used. The techniques based on Laplacian smoothing relocate the nodes recursively, based on the relative position of each node with respect to its adjacent nodes. This technique works quite well for meshes in convex domains, but it does not work so well in domains with multiple concavities, or that present interior islands. Iterative processes are relatively slow because the distortions of elements propagate to subsequent ones, requiring the establishment of a minimum number of iterations to obtain a reasonable result. In order to make the process faster, it is possible to introduce weights to each neighbouring node, on the basis of geometric factors, or to establish restrictions in the movement of the nodes [9]. Whichever option is taken, the method of moving nodes on the basis of their proximity to their neighbours, continues to be a heuristic procedure. In contrast, the methods based on optimization-based smoothing [7,13] establish a function that measures the global distortion of the mesh, and they determine the new coordinates of the nodes by minimising this function. Optimizationbased smoothing techniques are more expensive than the Laplacian ones, but they give better results in concave zones. Mixed methods are used for very big meshes,
synchronizing Laplacian and Optimization-based smoothing techniques, using the latter locally [3]. Finally, smoothing techniques which try to take advantage of physical phenomena have also been used. Lohner [11] developed a method that considers the edges of the mesh composed of an assembly of springs that join the nodes between them, and over which forces that tend to obtain a more regular form in each element are applied. Shimada [14] viewed the nodes as being the centre of an ideal bubble, which is relocated on the basis of its contact with other bubbles until equilibrium is achieved.
3. Method for quadrilateral mesh smoothing We have developed an optimization smoothing technique based on a physical analogy, where each element of the mesh is considered a deformable material, such that the smoothed coordinates are simply the degrees of freedom. Therefore, an elementary geometrical stiffness matrix for each quadrilateral has been defined, and by assembling the elements we obtain a matrix of the global stiffness of the system, by which the location of the new nodes of the mesh are calculated. The optimization-based smoothing technique is implemented in the following stages: 1. Estimation of the degree of uniformity of the initial mesh. 2. Calculation of the element stiffness matrix for each quad. 3. Assembly of the matrix of global stiffness of the mesh. 4. Imposition of the boundary conditions. 5. Computation of MckeeÕs method for the solution of the equationÕs system. 6. Estimation of the degree of uniformity of the smoothed mesh. In addition, it is possible to modify the global stiffness matrix in order to consider different boundary conditions to restrict the movements of the nodes. Below we develop the system of equations for the re-solution of the problem. 3.1. Element stiffness matrix of a quad The main idea of the smoothing technique is to consider the mesh as a deformable solid whose stiffness
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
depends on the degree of geometric distortion of each element. For its determination, the methodology developed by Balendran is followed [1], but we introduce a shape factor. For a quadrilateral mesh, let a given quad have vertices Pi Pj Pk Pl (see Fig. 2(a)). The ideal local relationship between the coordinates of the vertices is when the two diagonals are perpendicular to each other and, in addition, intersect at the midpoint. Since this will rarely happen (see Fig. 2(b)), we will have to take these nodes to positions Pi0 Pj0 Pk0 Pl0 (smoothed points) so that the ideal conditions are fulfilled (see Fig. 2(c)). If we consider one of the diagonals as fixed, for example Pj Pl (see Fig. 3(a)), points Pi and Pk are moved until they lie on the median line of Pj Pl . The distance of the new smoothed points Pi0 and Pk0 to the midpoint of the diagonal can vary, depending on whether we want the improved quads to tend towards rhombic or square forms (this is controlled by a shape factor l, the ratio between the magnitudes of the segments Pi0 Pk0 and Pj Pl ). l¼
49
(b)
(a)
Fig. 3. Smoothed element by deforming the initial quad.
x0k yk0
¼
1 xj þ xl 1 ðyl yj Þ l xl xj 2 yj þ yl 2
ð8Þ
and then xj þ lyj þ 2x0k xl lyl ¼ 0;
ð9Þ
lxj yj þ 2yk0 þ lxl yl ¼ 0:
ð10Þ
Once the mesh is smoothed, we will have
Pi0 Pk0 k~ P0 ~ P 0k k ¼ i : Pj Pl k~ Pj ~ P lk
ð1Þ
From the metric properties of Fig. 3(a) in node Pi0 it can be concluded that: ~ P 0i ¼ ~ P m þ 12 k~ P 0i ~ P 0k k ~ u;
ð2Þ
where Pm is the mid point of diagonal Pj Pl , and ~ u is the unit vector perpendicular to this diagonal. 1 ðyl yj Þ ~ u¼ ; ð3Þ xl xj P lk k~ Pj ~ x and y being the values of the coordinates of the points indicated in each subscript. If we express each component of Eq. (2) we have: 0 1 xj þ xl 1 ðyl yj Þ xi ¼ þ l ð4Þ xl xj yi0 2 yj þ yl 2 and then 2x0i xj lyj xl þ lyl ¼ 0;
ð5Þ
2yi0
ð6Þ
þ lxj yj lxl yl ¼ 0:
For node Pk0 we have ~ P 0k ¼ ~ P m 12 k~ P 0i ~ P 0k k ~ u;
ð7Þ
(a)
x0t xt ¼ 0 yt0 yt ¼ 0
t ¼ i; j; k; l
ð11Þ
enabling us to deduce from the previous relations the first, second, fifth and sixth rows of the matrix equation: 3 2 03 2 xi 2 0 1 l 0 0 1 l 6 07 7 6 0 2 l 1 0 0 l 1 7 6 yi 7 6 7 6 07 6 6 1 k xj 7 2 0 1 k 0 0 7 6 7 7 6 6 6 07 7 6 k 1 0 6 y 2 k 1 0 0 7 6 j7 6 7 7 6 6 0 x0k 7 0 1 l 2 0 1 l 7 6 7 7 6 6 6 07 6 0 yk 7 0 l 1 0 2 l 1 7 7 6 6 7 6 07 6 7 4 1 k 0 0 1 k 2 0 5 6 4 xl 5 1 2 3 0 607 6 7 6 7 607 6 7 607 6 7 ¼ 6 7: 607 6 7 607 6 7 6 7 405
k
0
(b) Fig. 2. Smoothing of a quad.
(c)
0
0
k
1
0
2
yl0
ð12Þ
50
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
The rest of equations can be determined using the other diagonal Pi Pk as fixed (see Fig. 3(b)) where k¼
Pj0 Pl0
ð13Þ
:
Pi Pk
Considering the symmetry of the problem k ffi l. When l ¼ 1 we will have dependency relationships that will transform a quad into a square, while for the rest of the cases the smoothed quad will tend to have rhombic forms. In this study l is approximated by a geometric constant, and it depend only of the initial geometry of the quad, as will be justified in Section 5. The elementary stiffness matrix k e ½8 8 of a generic quad is determined by the previous matrix relationship. 3 2 2 0 1 l 0 0 1 l 6 0 2 l 1 0 0 l 1 7 7 6 6 1 l 2 0 1 l 0 0 7 7 6 6 l 1 0 2 l 1 0 0 7 e 7: 6 k ¼6 0 1 l 2 0 1 l 7 7 6 0 6 0 0 l 1 0 2 l 1 7 7 6 4 1 l 0 0 1 l 2 0 5 l 1 0 0 l 1 0 2 ð14Þ This symmetric matrix can be expressed in terms of nodal sub-matrices as 2 3 kii kij kik kil 6 kji kjj kjk kjl 7 7 ke ¼ 6 ð15Þ 4 kki kkj kkk kkl 5: kli kjl klk kll
kie being the matrix of elementary stiffness of the quad i properly referred to its degrees of freedom considering the global identification of the nodes. The degrees of freedom are the smoothing coordinates, and they have been put in the order: ½x01 y10 x02 y20 x0m ym0 ;
ð20Þ
where m is the number of nodes. The assembled matrix is presented in Fig. 4. We have assumed that the numbering sequence of the nodes is i < k < l < j, in order to indicate the global assembling of the matrix elements. Actually, when assembling the matrix of global stiffness, we consider how each element affects the surrounding elements, thereby relating the degrees of freedom of each quad with its adjacent ones. Intuitively, it can be observed that the greater the number of quads sharing a particular node, the greater the stiffness of that node and, therefore, the greater the weight in the determination of the smoothed coordinates of the adjacent nodes. Once the global stiffness matrix of the system is assembled, we must solve the sparse system of linear algebraic equations defined by K~ c0 ¼ ~ b;
ð21Þ T ½x01 y10 x02 y20 x0m ym0
where ~ c0 ¼ the smoothed coordinates, right hand side, and K the dimension 2m 2m.
is the solution vector with T ~ b ¼ ½b1 b2 b3 b4 b2m the global stiffness matrix of
3.3. Assigning boundary conditions
Letting 2 D¼ 0
0 2
O¼
0 0 0 0
L¼
1 l
l ; 1
ð16Þ
The system of equation (Eq. (21)) would be indeterminate unless some of the nodes of the mesh were fixed.
we have kii ¼ kjj ¼ kkk ¼ kll ¼ D
kij ¼ kjk ¼ kkl ¼ kli ¼ L;
kik ¼ kki ¼ kjl ¼ klj ¼ O
kji ¼ kkj ¼ klk ¼ kil ¼ LT ;
x'i y'i
x'k y'k
x'l y'l x' j y' j
ð17Þ T
where L is the transpose matrix of L. This symmetric matrix k e can be expressed in terms of the nodal submatrices in the form: 2 3 D L 0 LT 6 LT D L 0 7 7: ð18Þ ke ¼ 6 40 LT D L 5 T L 0 L D
x'i y'i
ki i
ki k
ki l ki j
x'k y'k
kk i
kkk
kk l kk j
3.2. Assembly of the global stiffness matrix
x'l y' l x'j y'j
kl i
kl k
kl l kl j
kj i
k jk
kj l kj j
The matrix of global stiffness K for n quads will be assembled as X K¼ kie ; ð19Þ i¼1;n
Fig. 4. Assembling of the matrix of global stiffness K.
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
The fixed conditions are applied to the nodes of the boundary of the zone to be smoothed, or to any fixed inner nodes. If a node is free (non-fixed), the corresponding degrees of freedom in vector ~ b will be zero. Specific values will only appear if the node is fixed. The value of the latter can be physically interpreted as a force that prevents the point from moving due to the initial geometric distortion of the mesh. In order to handle the fixed value c0fix in the solution of the system of equations we proceed as in the methods of structural analysis; the stiffness in the main diagonal K½n; n of degree of freedom n is replaced by a very high fictitious stiffness (k1 1020 ), and simultaneously this fictitious stiffness multiplied by the constant value is introduced in the RHS of the system of equations (see Eq. (22)). b½n ¼ k1 c0fix :
K½n; n ¼ k1
ð22Þ
A fictitious stiffness of this value can be considered for numerical effects as infinity. When the system of equations is solved, the outcome of the value of the smoothed coordinate for that point is the same as the initial fixed value c0 ½n ¼ c0fix . In this way, the x or y coordinates of a node can be fixed since each degree of freedom is independent. In the case of the outer boundaries and inner boundaries of islands or isolated quads and segregated zones, the normal procedure is to fix both degrees of freedom of the node. A great advantage of this smoothing technique, is that boundary conditions are introduced only by altering the global stiffness of the nodes involved. This method is not limited to boundaries parallel to the x and y axes (see Fig. 7). If there is a fixed lineal dependency between the degrees of freedom n and g, given by the ratio rfix and the constant value sfix , see Eq. (23), we input this relationship using the fictitious stiffness given by Eq. (24) for K and Eq. (25) for ~ b. c0 ½n þ rfix c0 ½g ¼ sfix ; K½n; n ¼ k1 K½g; g ¼
2 rfix
ð23Þ
K½n; g ¼ K½g; n ¼ rfix k1 k1 ;
b½n ¼ sfix k1
b½g ¼ rfix sfix k1 :
ð24Þ ð25Þ
3.4. Resolution of the system of equations In Eq. (21) most of the terms of the matrix of global stiffness K are zero, since those degrees of freedom whose nodes do not share the same quad, are not bound by any stiffness, and therefore the input in the matrix of stiffness corresponding to those degrees of freedom is zero. A judicious reordering of the rows and columns for the coefficients matrix can lead to enormous savings in computer execution time and storage. Thus, we have
51
proposed a permutation that allows the rearrangement of the degrees of freedom, such that the same matrix can be stored in the most economic form. The system considered to implement the permutation has been the Mckee algorithm, while for the storage of the final matrix, the band and envelope method has been used [8]. Once the matrix K is stored, the system of equations is solved through the conventional Cholesky method, restoring the original equation order later, with the purpose of reassigning the values of smoothed coordinates to each node. If the matrix K is not symmetric and positive definite, some form of pivoting (row and/or column interchanges) is necessary to ensure the numerical stability of CholeskyÕs method. In order to avoid such difficulties, we have assumed k ¼ l and 0 < l 6 1; in this conditions K is symmetric and positive definite. 3.5. Two element mesh example In order to illustrate the method used in this work, we have considered the mesh of two quads of the Fig. 5. Assembly of the stiffness matrices for the two quads and using Eq. (22), leads to the following 12 equations: 2
4
6 6 0 6 6 2 6 6 6 0 6 6 0 6 6 6 0 6 6 1 6 6 6 l 6 6 1 6 6 6 l 6 6 0 4 0
0 4 0 2 0 0 l 1 l 1 0
2
0
0
0
1
0
2
0
0
l 1
l
1 l
0
1
0
0
3 2
x01
3
7 6 7 0 7 6 y10 7 7 6 7 6 07 4 0 1 l 0 0 0 0 1 l 7 7 6 x2 7 7 6 07 0 4 l 1 0 0 0 0 l 1 7 6 y2 7 7 6 7 6 07 1 l k1 0 1 l 0 0 0 0 7 7 6 x3 7 7 6 07 l 1 0 k1 l 1 0 0 0 0 7 6 y3 7 76 7 6 07 0 0 1 l k1 0 0 0 0 0 7 7 6 x4 7 7 6 07 6 7 0 0 l 1 0 k1 0 0 0 0 7 7 6 y4 7 6 07 0 0 0 0 0 0 k1 0 1 l 7 7 6 x5 7 7 6 07 6 7 0 0 0 0 0 0 0 k1 l 1 7 7 6 y5 7 7 0 7 1 l 0 0 0 0 1 l k1 0 5 6 4 x6 5 y60 l 1 0 0 0 0 l 1 0 k1 3
0 2 0 7 6 60 7 7 6 60 7 7 6 7 6 60 7 7 6 6 0k 7 6 17 7 6 6 1k1 7 7 ¼6 6 0k 7: 6 17 7 6 6 0k1 7 7 6 6 2k 7 6 17 7 6 6 0k1 7 7 6 6 2k 7 4 15 1k1
l
ð26Þ
Setting l ¼ 1 and k1 ¼ 1020 and solving the equations gives T ~ c0 ¼ ½0; 1; 1; 1; 0; 1; 0; 0; 2; 0; 2; 1 :
ð27Þ
52
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
Fig. 5. Initial mesh and smoothed mesh of two quads.
3.6. Degree of uniformity of the mesh Before and after the mesh has been smoothed, its degree of uniformity is calculated. To do this, it is important to select a suitable metric. When the verification of the value of the degree of uniformity of an element is above a certain threshold, the corresponding nodes are not included, and therefore, the size of the matrix of global stiffness is reduced. In this paper, we have taken a as the degree of uniformity, for a triangle with nodes Pi Pj Pk , as proposed by Lee and Lo [10] (see Fig. 6(a)), that is ! ! pffiffiffi kPk Pi Pk Pj k a ¼ 2 3 ! ð28Þ ! 2 ! 2 : 2 kPk Pi k þ kPi Pj k þ kPj Pk k To determine the degree of uniformity b of a quad (see Fig. 6(b) and (c)), it is divided into the two families of triangles bordered by their diagonals. The degree of triangular uniformity for each one of the four triangles is sorted in descending order of magnitude (a1 P a2 P a3 P a4 ). Finally, the value of the uniformity of the quad is a3 a4 b¼ : ð29Þ a1 a2
4. Mesh quality improvement results for the smoothing technique The scheme of smoothing developed here has been implemented in C++ within the Almec2D program, and it has been used to treat different meshes generated by this program. In Figs. 7–10 several tunnel models are illustrated indicating on the left the initial mesh, and on the right of the final smoothed mesh. As can be observed, the improvement in the mesh is very significant, and it tries to introduce the morphology of the tunnel cross-section very gradually. The minimum
(a)
(b) Fig. 7. Initial mesh (a) and smoothed mesh (b) of model A.
This degree of uniformity is 6 1 and approaches unity as the quad approaches a square.
(a) Fig. 6. Degree of uniformity a of the meshÕs elements.
(b) Fig. 8. Initial mesh (a) and smoothed mesh (b) of model B.
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
(a)
53
(b) Fig. 9. Initial mesh (a) and smoothed mesh (b) of model C.
Fig. 11. Degree of uniformity before smoothing (model D).
(a)
(b) Fig. 10. Initial mesh (a) and smoothed mesh (b) of model D.
and average values of b for each model are collected in Table 1, before and after making the smoothing. In the last column, we point out the increase of b in percentage. An increase of the average degree of uniformity is seen in all the models. In the best case, it is of the order of 4.5% (model A), although improvement values from 1.5% to 3.0% seem more representative (models C and D), as these are of the same order as improvements obtained with other smoothing algorithms [3,10]. The average value of the degree of uniformity does not clearly illustrate the improvement that takes place in the mesh. In order to overcome this difficulty we have to present the degree of local uniformity of each quadrilateral by mean of a map of isovalues. The zones are represented by colours ranging from red for near zero b values to blue for values close to unity. This legend and the degree of uniformity before smoothing are represented in Fig. 11, and after smoothing in Fig. 12. The improvement in the mesh can be observed clearly especially in the quadrilaterals adjacent to the boundaries of the sub-meshed into which the initial model was divided (see Fig. 1(a)).
Fig. 12. Degree of uniformity after smoothing (model D).
The results obtained for the rest of models are illustrated in Figs. 13–15.
5. Analysis of the shape factor In all the previous examples we have used a shape factor l defined by Eq. (32). Nevertheless, to complete the work we have to justify this choice. To do this, different l definitions have been studied in order to see
Table 1 Degree of uniformity of the meshes Model
A B C D
Number of nodes
440 1089 1790 1114
Number of quads
400 1000 1688 1038
Db (%)
Degree of uniformity b Before smoothing
After smoothing
Min
Average
Min
Average
0.334 0.042 0.254 0.391
0.699 0.728 0.884 0.858
0.363 0.050 0.278 0.356
0.744 0.737 0.905 0.878
4.5 0.9 2.1 2.8
54
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
Fig. 13. Degree of uniformity of initial mesh (a) and smoothed mesh (b) of model A.
Fig. 14. Degree of uniformity of initial mesh (a) and smoothed mesh (b) of model B.
Fig. 15. Degree of uniformity of initial mesh (a) and smoothed mesh (b) of model C.
which one gives the greatest degree of uniformity. In particular, four cases have been studied, as shown in Table 2. In case 1, l is constant for all quadrilaterals, in case 2 the relationship between the minimum and maxi-
mum diagonals is considered, in case 3 the relationship between the minimum and maximum side, and finally case 4, l is taken equal to the value of the degree of uniformity of the quadrilateral b given by Eq. (27).
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56 Table 2 Degree of uniformity and shape factor for each case Case
Maximum b
Shape factor l
1
0.772
l¼1
ð30Þ
2
0.795
minðPj Pl ; Pi Pk Þ l¼ maxðPj Pl ; Pi Pk Þ
ð31Þ
3
0.799
4
0.793
l¼
minðPi Pj ; Pj Pk ; Pk Pl ; Pl Pi Þ maxðPi Pj ; Pj Pk ; Pk Pl ; Pl Pi Þ
l¼b
ð32Þ ð33Þ
55
minimizes the user time required to produce an acceptable mesh in a simple way. It is a direct method that does not need iterative techniques of calculation and whenever the boundary conditions are well defined, the degree of uniformity improves remarkably. The shape factor helps to increase the degree of uniformity, especially when Eq. (32) is taken for its definition. The disadvantage is the high cost of memory when working with a matrix of global stiffness, but even so, this cost is less than that needed to make the final finite element analysis. The use of MckeeÕs method to solve the system of equations has been specially useful to mitigate this difficulty. This method is used within Almec2D to obtain more regular meshes, and it is one of the techniques for improving the numerical results of the calculation module. The possibility of extending this method to hexahedral meshes in 3D for the Almec3D program is being analyzed.
Acknowledgements Fig. 16. Degree of uniformity of cases 1–2–3–4.
In order to study several meshes, the simplest mesh has been taken as a reference, (the one corresponding to model A), and the density of nodes has been increased in the radial direction to obtain a series of eight meshes, identified by their number of nodes: 120, 200, 280, 360, 440, 520, 600 and 680 nodes. In Fig. 16 the values obtained for the degree of uniformity for eight consecutive meshes are represented. In the second column of Table 2, the maximum values reached by the degree of uniformity b are indicated. Initially b increases proportionally with the increase in the number of nodes, converging on a value of 0.8. However, it can be observed that from 500 nodes onwards, no further improvement is obtained in the degree of uniformity, regardless of the value l. Case 1 behaves better for meshes with a smaller number of nodes, but cases 2, 3 and 4 give the best results for meshes of more than 450 nodes. The greatest increase of average uniformity degree is obtained for case 3 (see Table 2), and for this reason we have used the value of l given by Eq. (32).
6. Conclusions An efficient 2D smoothing method has been developed. The method works with quadrilateral meshes by considering a criterion of global optimization that
The authors gratefully acknowledge the support of Robin Walker for the preparation of this paper in English.
References [1] Balendran B. A direct smoothing method for surface meshes. In: Proceedings of the Eighth International Meshing Roundtable, South Lake Tahoe, CA, USA 1999. p. 189–93. [2] Blacker TD, Stephenson MB. Paving: a new approach to automated quadrilateral mesh generation. Int J Numer Meth Eng 1991;32:811–47. [3] Canann Scott A, Tristano JR, Staten ML. An approach to combined Laplacian and optimization-based smoothing for triangular, quadrilateral, and quad-dominant meshes. In: Seventh International Meshing Roundtable. Sandia National Labs; 1998. p. 479–494. [4] Edelsbrunner H, Shah N. Incremental topological flipping works for regular triangulations. In: Proceedings of the Eighth ACM Symposium on Computational Geometry, 1992. p. 43–52. [5] Faux ID, Pratt MJ. Computational geometry for design and manufacture. Chichester, West Sussex, England: John Wiley & Sons; 1987. [6] Field D. Laplacian smoothing and Delaunay triangulations. Commun Numer Meth Eng 1988;4:709–12. [7] Freitag L. On combining Laplacian and optimization-based mesh smoothing techniques AMD trends in unstructured mesh generation. In: ASME, vol. 220. 1997. p. 37–43. [8] George A, Liu J. Computer solution of large sparse positive definite systems. New Jersey: Prentice-Hall; 1981. [9] George PL, Borouchaki H. Delaunay triangulation and meshing. Application to finite elements. Paris, France: Editions Hermes; 1998. p. 230–4.
56
A. Menendez Dıaz et al. / Computers and Geotechnics 31 (2004) 47–56
[10] Lee CK, Lo SH. A new scheme for the generation of a graded quadrilateral mesh. Comput Struct 1994;52(5):847–57. [11] Lohner R, Morgan K, Zienkiewicz OC. Adaptive grid refinement for the compressible Euler equations. In: Accuracy estimates and adaptive refinements in finite element computations. New York: Wiley; 1986. p. 281–97. [12] Menendez Dıaz A, Gonzalez Nicieza C, Alvarez Vigil A, Morıs Menendez G. Malladores Geometricos para el Disesimno Estructural por Diferencias Finitas, Disegno di Macchine e Progettazione Industriale, Junio 1996. p. 19–32.
[13] Parthasarathy V, Kodiyalam S. A constrained optimization approach to finite element mesh smoothing. Finite Elements Anal Des 1991;9:309–20. [14] Shimada K. Physically-based mesh generation: Automated triangulation of surfaces and volumes via bubble packing, Tesis Doctoral, ME Department, Massachusetts Institute of Technology, Cambridge, MA, 1993. [15] Staten ML, Canann SA. Post refinement element shape improvement for quadrilateral meshes, trends in unstructured mesh generation. ASME 1997;220:9–16.