New method for graded mesh generation of all hexahedral finite elements

Computers and Structures 76 (2000) 729±740

www.elsevier.com/locate/compstruc

New method for graded mesh generation of all hexahedral ®nite elements p Hua Li, Gengdong Cheng* State Key Laboratory of Structural Analysis of Industrial Equipment, Dalian University of Technology, Dalian 116023, People's Republic of China Received 25 June 1998; accepted 27 June 1999

Abstract Mapping method is widely applied by most of commercial mesh generators because of its eciency, mesh quality. One of the obstacles to apply the mapping method and generate a graded all hexahedral mesh of high quality in an arbitrarily three-dimensional domain is the generation of hexahedral parent elements on a super-element that allows for gradations in three co-ordinate directions. This paper presents a pattern module's method to generate the graded mesh of all hexahedral elements in a cube and thus improves the mapping method. The method requires few calculations. 7 2000 Elsevier Science Ltd. All rights reserved. Keywords: Hexahedral; Finite element; Mesh; Mapping; Pattern module; Conforming

1. Introduction The ®nite element method is a powerful and versatile numerical analysis tool to handle complicated structures, but its usefulness is hampered by the discretization of general geometry of the problem into a valid ®nite element mesh. The discretization process is tedious and error-prone if done manually. Furthermore, ®nite element analysis softwares have been becoming a part of integrated CAD software. In order to analyze a structure generated by CAD system, it is needed to achieve an automatic mesh generation that is integrated with an entity construction technology. High demand for automatic mesh generation is also from the area of non-linear structural analysis such as

The project (No. 19802005) is supported by the National Natural Science Foundation of China. * Corresponding author. p

metal forming, polymer manufacture process simulation or so. For di€erent physical problems, we are asked to divide these di€erent objects, such as plane domain, surface domain or three-dimensional (3D)object into ®nite element meshes. Because di€erent physical problems may need mesh of di€erent characteristics and di€erent quality, an excellent mesh generation method should adapt arbitrary boundary, nonconvex domain, and can control the mesh density, the direction of mesh lines and the location of speci®c mesh node. In recognition of these problems, a large number of methods have been proposed to automate the mesh generation task. In 2D situations this goal has been realized [1]. But it is much more dicult to handle 3D problem, especially in the generation of hexahedral ®nite element mesh. It is well known that hexahedral element is more accurate than the linear, constant-strain tetrahedral element and is thus more attractive for users of FEM. Lo [2] had reviewed the methods of tetrahedral mesh

0045-7949/00/$ - see front matter 7 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 9 ) 0 0 1 9 3 - 5

730

H. Li, G. Cheng / Computers and Structures 76 (2000) 729±740

generation. Compared with the tetrahedral element method, hexahedral mesh generation is more dicult. Landertshamer [3] gave an approach to generate 3D boundary-®tted hexahedral meshes, but it needs huge e€ort to adapt various corner points, boundaries, surfaces and the user has a high degree of in¯uence on the quality of the ®nal mesh. Shephard [4] and Perucchio [5] used the octree technique to divide the 3D domain into the mixed meshes of hexahedral and tetrahedral elements, but the shape of elements near the boundaries might be bad. Li and Armstrong [6] used the midpoint subdivision method and the integer programming to divide a simple subregion into hexahedral, but the subregion must satisfy certain simple topological condition. And also there are no e€ective methods to merge tetrahedral elements into hexahedral elements. Mapping element method is the mainstay of the existing commercial mesh generators. This method is fast. It can control mesh topology and density, and handle arbitrary boundary and easily meet various physical demands and the result elements have ®ne shape. The ®rst step of the mapping method for 3D cases is to divide a 3D solid into a number of subregions, which are hexahedral super-element de®ned by its 12 edges [7]. A typical 3D super-element in the physical space is shown in Fig. 1. By a trans®nite mapping Eq. (1), the super-element is transferred into a cube in the local x±Z±g co-ordinate. p…x, Z, g† ˆ …1 ÿ Z†…1 ÿ g†h1 …x† ‡ …1 ÿ Z†gh2 …x† ‡ Zgh3 …x† ‡ Z…1 ÿ g†h4 …x† ‡ …1 ÿ x†…1 ÿ g†h5 …Z † ‡ …1 ÿ x†gh6 …Z † ‡ xgh7 …Z † ‡ x…1 ÿ g†h8 …Z † ‡ …1 ÿ x†…1 ÿ Z†h9 …g † ‡ …1 ÿ x†Zh10 …g † ‡ xZh11 …g † ‡ x…1 ÿ Z†h12 …g † ‡ c…x, Z, g†

Fig. 1. Cartesian and body co-ordinate systems and boundary edge functions.

 c…x, Z, g† ˆ ÿ2 …1 ÿ x†…1 ÿ Z†…1 ÿ g†v…0, 0, 0† ‡ …1 ÿ x†…1 ÿ Z†gv…0, 0, 1† ‡ …1 ÿ x†Z…1 ÿ g†v…0, 1, 0† ‡ …1 ÿ x†Zgv…0, 1, 1†

…1†

‡ x…1 ÿ Z†…1 ÿ g†v…1, 0, 0† ‡ x…1 ÿ Z†gv…1, 0, 1† ‡ xZ…1 ÿ g†v…1, 1, 0† ‡ xZgv…1, 1, 1†



Where p…x, Z, g† ˆ fx…x; Z, g†, y…x, Z, g†, z…x, Z, g†g; hi …s† …i ˆ 1 to 12; s ˆ x, Z, g† are the boundary edge functions that specify x, y and z along the ith boundary edge. And a typical boundary edge function is a spline function; v…0, 0, 0†, v…0, 0, 1†, v…0, 1, 0†, v…0, 1, 1†, v…1, 0, 0†, v…1, 0, 1†, v…1, 1, 0†, v…1, 1, 1† are the eight nodes and v…i, j, k† ˆ fx…i; j, k†, y…i, j, k†, z…i, j, k†g …i, j, k ˆ 0 or 1†: The co-ordinates (x, y, z ) of a typical point within the super-element is related to the nodal co-ordinates on the 12 boundary edges and the coordinates …x, Z, g† of its corresponding point by Eq. (1). The second step of the mapping method is to generate mesh in the cube according to the requirement of the number of elements to be generated. And the last step is to transfer the mesh in the cube into the physical space by using Eq. (1). It is well recognised that for a physical problem, some areas of the structure may need mesh of high density, and others require low mesh density. Therefore, the mesh gradation problem, i.e., the gradual change of mesh density is one of the most important problems faced by all mesh generation techniques. Furthermore, the meshes of adjacent super-elements in the mapping method should be conforming in their interfaces where adjacent elements share a common side or a common face. In order to meet the requirement of both the mesh gradation and the conformability of the meshes between two adjacent super-elements, it becomes very critical for the success of the mapping method that one manages to generate the graded mesh in each super-element under very general conditions. Since the mesh in super-element is obtained by the mapping method from the mesh in a cube, the basic problem is to generate a graded mesh in a cube under the prescribed number of nodes on some sides and the element distribution on some faces of the cube. When the number of nodes are equal on opposite sides of the cube, mesh generation is trivial. When the number of nodes are di€erent on the opposite sides of the cube, various approaches to generate graded mesh composed of tetrahedral can be found in literatures, for examples, Kadivar and Shari® [8] applied the double mapping of isoparametric mesh generation to deal with the gradation problem. But the task is more dicult with a mesh of hexahedral elements. The present paper extends our previous study on mesh gener-

H. Li, G. Cheng / Computers and Structures 76 (2000) 729±740

731

ation in 2D domain [1] and presents a new pattern module's method that allows for the mesh gradation in three co-ordinate directions and generates all hexahedral elements. Moreover, we formulate an integer programming for ®nding the parameters of the pattern module's method and solve this problem by the linear integer goal programming technology. In this way we avoid the strict restriction on the distribution of the mesh density in the mapping method and improve the mapping element method. Fig. 2. A cube and its 12 sides.

2. All-hexahedral mesh generation of a cube Let us consider the basic problem of graded all-hexahedral mesh generation of a cube. The 12 sides of the cube are numbered as Fig. 2. The numbers of element divisions on 12 sides are Ni …i ˆ 1 to 12† respectively. Ni …i ˆ 1 to 12† must meet condition 1 in order to ensure that the elements on 6 faces are all quadrilateral element which is necessary to generate all hexahedral elements. Condition 1. N1 ‡ N2 ‡ N3 ‡ N4 ˆ even number N1 ‡ N5 ‡ N8 ‡ N9 ˆ even number N2 ‡ N5 ‡ N6 ‡ N10 ˆ even number N3 ‡ N6 ‡ N7 ‡ N11 ˆ even number N4 ‡ N7 ‡ N8 ‡ N12 ˆ even number N9 ‡ N10 ‡ N11 ‡ N12 ˆ even number

2.1. The pattern module

Fig. 4a). L pattern includes 12 junction sides. On the four parallel junction sides, the number of element divisions are equal to each other.

2.1.2. Corner patterns A corner pattern is applied to the eight corner regions. It consists of a corner element and a certain number of outer layers (see Fig. 5). The corner element is a cuboid and each outer layer includes three hexahedral elements. In order to describe the method easily, we name the three sides that are on the sides of the cube as junction sides of corner pattern (see the black and thick lines in Fig. 5b). In 12 sides of a corner pattern, three junction sides have the same number of element divisions and the number of element divisions of the nine others are all one. The corner pattern is the basic pattern of mesh subdivision of the corner region. Since the eight corner regions are di€erent only in the position and the direction, the pattern shown in Fig. 5a, which is for the region B, can be applied to any of the eight corner regions through translation and rotation. The union of the junction sides of the corner patterns in the eight corner regions and the L pattern in the L region constitutes the 12 sides of the cube.

Let us divide the cube into nine 3D regions: A, B, C, D, E, F, G, H and L, see Fig. 3. Among the 9 regions, A, B, C, D, E, F, G, H are corner region and L is central region. Further, let us design three patterns that are applied to di€erent region. 2.1.1. L pattern L pattern is applied to the region L. It only includes cuboid elements (see Fig. 4a) and its projections on three co-ordinate planes are similar (see Fig. 4b). In order to describe the method easily, we name the sides as junction sides of L pattern that are on the boundaries of the cube (see the black and thick lines in

Fig. 3. Nine regions in a cube.

732

H. Li, G. Cheng / Computers and Structures 76 (2000) 729±740

Fig. 4. The L pattern and its projection outline on the co-ordinate planes.

2.1.3. Side patterns A side pattern consists of an inner element and a certain number of outer layers (see Fig. 6). The inner element is a cuboid element and each outer layer includes four hexahedral elements. In 12 sides of a side pattern, the numbers of element divisions of 11 sides are one. And there is only one side for which number of element division is larger than 1 and we name the side as junction side of side pattern (see the black and thick lines in Fig. 6). The side pattern is to be imbedded in the corner region and to be applied in the corner region. 2.1.4. The embedding rule of the side patterns In the present method, the side patterns are imbedded in the corner element of corner patterns, and its junction side is located on one of the junction side of the corner pattern. Depending on the situation of

the result mesh, the side pattern may be located on any junction sides of the corner pattern. Because the number of corner regions is 8, and each corner pattern in each corner region has three junction sides, the side pattern could be imbedded in the 3  8 = 24 possible positions. Fig. 7 shows the three cases in which the side pattern is located on the three junction sides of the corner pattern of A region corresponding to the 1, 4 and 8 sides of the cube (see the black and thick lines in Fig. 7), respectively. In this paper, the embedding rule of the side patterns is that one corner region at most includes one side pattern. In summary, we apply the L pattern in the L region, and a proper combination of the side pattern and the corner pattern in the eight corner region with the limitation that one corner region at most includes one side pattern.

Fig. 5. A corner pattern.

H. Li, G. Cheng / Computers and Structures 76 (2000) 729±740

733

structed in a cube for a given set of parameters of a pattern module. De®nition 4. Full degree side is a side which dose not include any side pattern.

2.2. Determination of parameters of a pattern modules

De®nition 5. A side is called as the shortest side if the degree number of two corner patterns which are applied to the two corner region and lie on the same side of the cube both are 1, and this side has no side pattern.

2.2.1. Terminology and de®nitions

Obviously, the shortest side is full degree side.

Fig. 6. A side pattern.

De®nition 1. Degree number of L pattern is the number of element divisions on the junction side of L pattern. If the numbers of element divisions on the junction side of L pattern in three co-ordinate directions are LX , LY , LZ , respectively, then the degree number of L pattern is LX , LY , LZ : De®nition 2. Degree number of corner pattern is the number of element divisions on the junction side of corner pattern. According to the region which the corner pattern is applied to, its corresponding degree number is denoted as A, B, C, D, E, F, G and H, respectively. De®nition 3. Degree number of side pattern is the number of its outer layers. The degree numbers of side patterns are denoted as Ei , in which i 2 f1; 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12g, and i refers to the side number of the cube which the junction side of the side pattern lie on. For example, in Fig. 7a, A = 2; E1 ˆ 1: The parameters of a pattern modules are LX , LY , LZ , A, B, C, D, E, F, G, H, and Ei for i 2 f1; 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12g: Following the above description, a mesh of all hexahedral elements may be con-

2.2.2. Basic relations of the degree number of patterns The result meshes depend on the degree number of patterns applied to each region. From the topologies of the 9 regions and the embedding rule, the degree numbers of patterns and Ni always meet the following basic relations LX ‡ A ‡ B ‡ 2E1 ˆ N1 LY ‡ B ‡ C ‡ 2E2 ˆ N2 LX ‡ C ‡ D ‡ 2E3 ˆ N3 LY ‡ A ‡ D ‡ 2E4 ˆ N4 LZ ‡ B ‡ F ‡ 2E5 ˆ N5 LZ ‡ C ‡ G ‡ 2E6 ˆ N6 LZ ‡ D ‡ H ‡ 2E7 ˆ N7 LZ ‡ A ‡ E ‡ 2E8 ˆ N8

Fig. 7. The side patterns in A pattern.

734

H. Li, G. Cheng / Computers and Structures 76 (2000) 729±740

LX ‡ E ‡ F ‡ 2E9 ˆ N9 LY ‡ F ‡ G ‡ 2E10 ˆ N10 LX ‡ G ‡ H ‡ 2E11 ˆ N11 LY ‡ E ‡ H ‡ 2E12 ˆ N12

…2†

The above relations (2) are very important in the present method. And depending on the requirement of mesh generation, the relations (2) can be further rearranged. 2.2.3. Mathematical formulation of di€erent mesh generation requirement For the mapping method, we must consider the relations between the two adjacent super-elements. If we design the mesh in each super-element without consideration of the adjacent super-elements, the mesh may not be conforming to the adjacent ones. Fig. 8 gives an example of two nonconforming meshes. The reason of nonconforming is that the degree numbers of patterns on the adjacent faces of two super-elements are di€erent. In super-element 1, A = 2, B = 1, C = 1, D = 1, LX ˆ 1, LY ˆ 1; in super-element 2, A 0 ˆ 1, B 0 ˆ 1, C 0 ˆ 1, D 0 ˆ 1, LX ˆ 2, LY ˆ 2: To describe the possible mesh generation requirement in a cube, four notions about sides and faces of a super-element ®rstly are de®ned.

Fig. 8. The nonconforming mesh where adjacent elements do not share a whole face.

Depending on the state of the super-element to be meshed, the basic problem for the cube may be di€erent. We will present the mathematical formulation for three cases. Case 1. If 12 sides of a super-element are all ®xed side and there is no conforming face, then A, B, C, D, E, F, G, H, LX , LY , LZ , Ei …i ˆ 1, 2, . . ., 12† are the unknowns, and Ni …i ˆ 1, 2, . . ., 12† are the known numbers. The mathematical formulation for determining the parameters of patterns module is as follows: To ®nd positive integer A, B, C, D, E, F, G, H and non-negative integer LX , LY , LZ , Ei …i ˆ 1, 2, . . ., 12† such that LX ‡ A ‡ B ‡ 2E1 ˆ N1

De®nition 6 (Free side). If a side of a super-element is not adjacent to any other super-elements, and the number of element division on the side is not ®xed, then the side is called free side.

LY ‡ B ‡ C ‡ 2E2 ˆ N2

De®nition 7 (Fixed side). If a side of a super-element is not adjacent to any other super-elements, and the number of element division on a side is given, then the side is called ®xed side

LY ‡ A ‡ D ‡ 2E4 ˆ N4

De®nition 8 (Conforming side). If the number of element division and the degree number of patterns on a side of a super-element have been de®ned by the adjacent super-elements and does not permit to be changed then the side is called conforming side.

LX ‡ C ‡ D ‡ 2E3 ˆ N3

LZ ‡ B ‡ F ‡ 2E5 ˆ N5 LZ ‡ C ‡ G ‡ 2E6 ˆ N6 LZ ‡ D ‡ H ‡ 2E7 ˆ N7

De®nition 9. Conforming face: If the number of element division and the degree number of patterns on all four sides of a face have been given by the adjacent superelement, then the face is called conforming face.

LZ ‡ A ‡ E ‡ 2E8 ˆ N8

Obviously, all sides located on a conforming face are conforming sides. If two super-elements only share a side and have not any common faces, this shared side is ®xed side and not conforming side.

LY ‡ F ‡ G ‡ 2E10 ˆ N10

LX ‡ E ‡ F ‡ 2E9 ˆ N9

LX ‡ G ‡ H ‡ 2E11 ˆ N11

H. Li, G. Cheng / Computers and Structures 76 (2000) 729±740

LY ‡ E ‡ H ‡ 2E12 ˆ N12

…3†

LZ ‡ B ‡ 2E5 ˆ N5 ÿ F

Case 2. If the face enclosed by side 1, 2, 3, 4 in Fig. 2 is conforming face, then the side 1, 2, 3, 4 are conforming side. Let the other sides be ®xed side. In this case, the known number are Ni …i ˆ 1, 2, . . ., 12†, A, B, C, D, LX , LY , and Ei …i ˆ 1, 2, 3, 4†, the unknown number are E, F, G, H, LZ , Ei …i ˆ 5, 6, . . ., 12†: The mathematical formulation for determining the pattern modules is as follows: To ®nd positive integer E, F, G, H and non-negative integer LZ , Ei …i ˆ 5, 6, . . ., 12† such that

LZ ‡ C ‡ 2E6 ˆ N6 ÿ G

735

LZ ‡ D ‡ 2E7 ˆ N7 ÿ H LZ ‡ A ‡ 2E8 ˆ N8 ÿ E H ‡ 2E11 rN11 ÿ LX ÿ G H ‡ 2E12 RN12 ÿ LY ÿ E

…5†

LZ ‡ F ‡ 2E5 ˆ N5 ÿ B Depending on the state of the super-elements, di€erent gradation and conformation requirement may exist. We may set up the mathematical formulation for determining the unknowns based on the Eq. (2) in the way similar to the above modi®cation. Solving the unknowns, we obtain the parameters of the pattern module and construct the mesh in the cube.

LZ ‡ G ‡ 2E6 ˆ N6 ÿ C LZ ‡ H ‡ 2E7 ˆ N7 ÿ D LZ ‡ E ‡ 2E8 ˆ N8 ÿ A E ‡ F ‡ 2E9 ˆ N9 ÿ LX F ‡ G ‡ 2E10 ˆ N10 ÿ LY G ‡ H ‡ 2E11 ˆ N11 ÿ LX E ‡ H ‡ 2E12 ˆ N12 ÿ LY

…4†

Case 3. If the side 1±8 are all ®xed side, side 9 and side 10 are conforming side, the side 11 and side 12 are free sides and the number of element division of side 11 is required to be not less than N11 and the number of element division of side 12 is required to be not large than N12, then A, B, C, D, H, LZ , Ei …i ˆ 1, 2, . . ., 8, 11, 12† are the unknowns, and E, F, G, LX , LY , E9, E10, Ni …i ˆ 1, 2, . . ., 12† are the known numbers. The mathematical formulation for determining the pattern modules is as follows: To ®nd positive integer A, B, C, D, H and non-negative integer LZ , Ei …i ˆ 1, 2, . . ., 8, 11, 12† such that A ‡ B ‡ 2E1 ˆ N1 ÿ LX B ‡ C ‡ 2E2 ˆ N2 ÿ LY C ‡ D ‡ 2E3 ˆ N3 ÿ LX A ‡ D ‡ 2E4 ˆ N4 ÿ LY

2.2.4. Linear integer goal programming (LIGP) model of pattern module's method Let x ˆ …x 1 , x 2 , . . ., x n †T be the unknowns. The superscript T denotes vector transposition. The Eqs. (3)±(5), after transforming the positive integer to nonnegative integer, can be described generally as follows: Find xÅ ˆ …x 1 , x 2 , . . ., x n † such that 8 n X > > > …xÅ † ˆ cij x j Rbi i ˆ 1,


, p f > i > > > jˆ1 > > > > n < X fi …xÅ † ˆ cij x j ˆ bi i ˆ p ‡ 1, . . . , q …6† > > jˆ1 > > > > n > X > > > …Åx † ˆ f cij x j rbi i ˆ q ‡ 1, . . . , m i > : jˆ1

and x j …j ˆ 1 to n† is non-negative integer. Where bi is the right-hand-side value of Eqs. (3), (4) or (5); ci, j is the coecient associated with variable j in the ith equation of Eqs. (3), (4) or (5). To solve the set equalities or inequalities, we apply the Linear Integer Goal Programming technique. We consider each Eq. (6) as an objective and introduce unknown negative deviation ni and positive deviation pi …i ˆ 1, 2, . . ., m†, which re¯ect the underachievement or overachievement of the ith objective. Eq. (6) can be reformulated as a Linear Integer Goal Programming (LIGP) model [9] as follows:  Find xˆ…x 1 , x 2 , . . ., x n † so as to ! p q m X X X min …7† pi ‡ ‡ n n p ‡ … i i† i iˆ1

iˆp‡1

iˆq‡1

736

H. Li, G. Cheng / Computers and Structures 76 (2000) 729±740

subject to

details, the reader is referred to [9]. If given N1±N12 are 5, 4, 7, 10, 7, 8, 6, 2, 6, 7, 5, 6 respectively and the degree number of patterns meet Eq. (3), then solution are as follows:

n X ci, j x j ‡ ni ÿ pi ˆ bi jˆ1

LX ˆ 1, LY ˆ 2, LZ ˆ 0

x j , pi , ni r0 and integer

A ˆ 1, B ˆ 1, C ˆ 1, D ˆ 3, E ˆ 1, F ˆ 4, G ˆ 1, H ˆ 3

i ˆ 1, 2, . . . , m Therefore, for example, the LIGP model of, for  …x 1 , x 2 , . . ., x 23 † so as to example, Eq. (3) is Find xˆ 12 X min …pi ‡ ni †

…8†

iˆ1

Subject to x 1 ‡ x 4 ‡ x 5 ‡ 2x 12 ‡ n1 ÿ p1 ˆ N1 ÿ 2 x 2 ‡ x 5 ‡ x 6 ‡ 2x 1 ‡ n2 ÿ p2 ˆ N2 ÿ 2 x 1 ‡ x 6 ‡ x 7 ‡ 2x 14 ‡ n3 ÿ p3 ˆ N3 ÿ 2 x 2 ‡ x 4 ‡ x 7 ‡ 2x 15 ‡ n4 ÿ p4 ˆ N4 ÿ 2 x 3 ‡ x 5 ‡ x 9 ‡ 2x 16 ‡ n5 ÿ p5 ˆ N5 ÿ 2 x 3 ‡ x 6 ‡ x 10 ‡ 2x 17 ‡ n6 ÿ p6 ˆ N6 ÿ 2 x 3 ‡ x 7 ‡ x 11 ‡ 2x 18 ‡ n7 ÿ p7 ˆ N7 ÿ 2 x 3 ‡ x 4 ‡ x 8 ‡ 2x 19 ‡ n8 ÿ p8 ˆ N8 ÿ 2 x 1 ‡ x 8 ‡ x 9 ‡ 2x 20 ‡ n9 ÿ p9 ˆ N9 ÿ 2 x 2 ‡ x 9 ‡ x 10 ‡ 2x 21 ‡ n10 ÿ p10 ˆ N10 ÿ 2 x 1 ‡ x 10 ‡ x 11 ‡ 2x 22 ‡ n11 ÿ p11 ˆ N11 ÿ 2 x 2 ‡ x 8 ‡ x 11 ‡ 2x 23 ‡ n12 ÿ p12 ˆ N12 ÿ 2 x i , pi , ni …i ˆ 1, 2, . . . , 12† are non ÿ negative integer number: The methodology for solving linear goal programs can not, in general, be used to solve linear integer goal programs since it allows the value of the decision variables to be fractional. All Integer Cutting Plane Algorithm will be used to solve the LIGP model. For the

E1 ˆ 1, E2 ˆ 0, E3 ˆ 1, E4 ˆ 2, E5 ˆ 1, E6 ˆ 3, E7 ˆ 0, E8 ˆ 0, E9 ˆ 0, E10 ˆ 0, E11 ˆ 0, E12 ˆ 0 For Eqs. (4) and (5), we can build two new LIGP models in a similar way and obtain the solutions. 2.2.5. The number of elements in the cube and each pattern The total number of elements in the cube and every pattern can be calculated from the degree numbers of each pattern. Denote the symbol j  j the number of elements in pattern , we have jL j ˆ LX
LY
LZ ‡ 4…LX ‡ LY ‡ LZ † jP j ˆ 3…P ÿ 1 † ‡ 1 jEi j ˆ 4…Ei ÿ 1† ‡ 1 In which, P $ { A, B, C, D, E, F, G, H} is the name of corner pattern; P 2 fA; B, C, D, E, F, G, H g is the degree number of corner pattern; Ei is the name of side pattern and Ei is the degree number of side pattern, i 2 f1; 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12g: The number of elements in the cube is jcube j ˆ LX
LY
LZ ‡ 4…LX ‡ LY ‡ LZ † ‡

X …3P ÿ 2 † P2fA, B, C, D, E, F, G, Hg

‡

12 X …4Ei ÿ 3† iˆ1

Fig. 9 shows three cubes. The numbers of nodes on each side of the three cubes are speci®ed in advance. The graded meshes are constructed by using the above pattern module method. In Fig. 9(a), the numbers of N1±N12 are 13, 5, 8, 4, 10, 9, 8, 7, 4, 8, 11, 9, respectively. The degree number of L pattern is: LX ˆ 0, LY ˆ 1, LZ ˆ 2; the degree number of each corner patterns are A = 2, B = 3, C = 1, D = 1, E = 3, F = 1, G = 6, H = 5, respectively; the degree number of each side patterns are E1 = 4, E3 = 3, E5 = 2, the others are 0; side pattern E1, E3, E5 are embedded in corner pattern A, C and B, respectively; the number of nodes and elements in the result mesh of the cube are 318 and 172, respectively.

Fig. 9. Several cubes and their mesh.

738

H. Li, G. Cheng / Computers and Structures 76 (2000) 729±740

Fig. 10. Neighbourhood of an internal node i.

Fig. 11. The smoothing mesh of a cube.

Fig. 12. Hexahedral mesh of a general 3D object (1).

H. Li, G. Cheng / Computers and Structures 76 (2000) 729±740

739

Fig. 13. Hexahedral mesh of a general 3D object (2).

In Fig. 9(b), the numbers of N1±N12 are 6, 9, 9, 4, 5, 5, 4, 6, 15, 9, 6, 10, respectively; In Fig. 9(c), the numbers of N1±N12 are 10, 9, 8, 7, 5, 5, 4, 8, 9, 7, 7, 9, respectively.

3. The mesh smoothing The elements generated in Section 2 are all hexahedral, but the elements are not well distributed. We apply the Eq. (9) to reposition an internal node i in the cube [10].

8 Ni > 1 X > > xi ˆ > …x na ‡ x nb ‡ x nc ÿ x nd ÿ x ne ÿ x nf ‡ x ng † > > Ni nˆ1 > > > > > Ni < 1 X yi ˆ …yna ‡ ynb ‡ ync ÿ ynd ÿ yne ÿ ynf ‡ yng † > Ni nˆ1 > > > > > Ni > > 1 X > > z ˆ > …zna ‡ znb ‡ znc ÿ znd ÿ zne ÿ znf ‡ zng † i : N i nˆ1

…9† Where Ni is the number of elements around the node i. the subscripts na, nb, nc, nd, ne, nf, ng refer to the neighbouring nodes, see Fig. 10. The nodes on 6 faces of a cube will be in¯uenced

740

H. Li, G. Cheng / Computers and Structures 76 (2000) 729±740

only by these nodes which are located on the same faces. Fig. 11 is the smoothing mesh of Fig. 9(a).

4. Algorithm and examples From the above study, an algorithm for a graded hexahedral mesh generation in an arbitrary domain is as follows: (a) Divided the entire 3D structure into a number of curved hexahedral subdomains and specify the number and position of nodes on each side of the subdomains. (b) For each subdomain, apply the method in Section 2 to generate a graded all hexahedral mesh in its parametric domain Ð a cube. (b.1) According to the state of a cube, building the LIGP model. (b.2) Using the all integer cutting plane algorithm obtain the degree number of each patterns. (b.3) Using embedding rule de®ne the embedded position of each side patterns. (b.4) Mesh L area and eight corner areas respectively. (b.5) Delete the repeat nodes and obtain the result mesh in a cube. (c) According to the position of nodes in the sides of subdomains to rearrange the position of corresponding node in cube, and smoothing the mesh in parametric domain. (d) Map the mesh back to the corresponding physical subdomains. Figs. 12 and 13 show two 3D complicated shapes, that are meshed by the above method.

5. Conclusion This present paper tries to improve the mapping method for mesh generation of 3D domain. The basic

strategy of the method is to develop a pattern module for generating a graded all-hexahedral mesh in a cube. The graded mesh in the parametric domain then is mapped to the physical domain. It has been shown by numerical experience that with the pattern module method the mesh gradation is very ¯exible and can be adjusted according to the practical requirement. In comparison with other approach, the present method needs much less calculation. In this way, the method keeps all merits of the mapping methods. However, the automation of ®rst step of the mapping method, i.e., to divide a 3D solid into a number of subregions that are hexahedral super-elements remains challenging.

References [1] Gengdong Cheng, Huo Li. New method for graded mesh generation of all quadrilateral ®nite element. Comp and Stru 1996;59:823±9. [2] Lo SH. Volume discretization into tetrahedra, Part II: 3D triangulation by advancing front approach. Comp and Stru 1991;39:501±11. [3] Landertshamer F. Method to generate complex computational meshes eciently. Communications in numerical methods in engineering 1994;10:373±84. [4] Shephard MS, Georges MK. Automatic three-dimensional mesh generation by the ®nite octree technique. SCOREC Report no. 1, 1991. [5] Perucchio R, Saxena M, Kela A. Automatic mesh generation from solid models based on recursive spatial decomposition. Int J Num Meth Engng 1988;24:2135±59. [6] Li CG, McKeag RM, Armstrong CG. Hexahedral meshing using midpoint subdivision and integer programing. Computer methods in applied mechanics and engineering 1995;124:171±93. [7] Cook WA. Body oriented (natural) co-ordinates for generating three dimensional meshes. Int J Num Meth Engng 1974;8:27±43. [8] Kadivar MH, Shari® H. Double mapping of isoparametric mesh generation. Comp and Stru 1996;59:471±7. [9] Ignizio JP. Goal programming and extensions, Lexington Books, 1976. [10] Herrmann LR. Laplacian-isoparametric grid generation scheme. J Engrg Mech Div 1976;102(EMS 10):749±56.