A

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng (2008) Published online in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/nme.2470

Octree-based reasonable-quality hexahedral mesh generation using a new set of refinement templates Yasushi Ito∗, †, ‡ , Alan M. Shih§ and Bharat K. Soni¶ Enabling Technology Laboratory (ETLab), Department of Mechanical Engineering, University of Alabama at Birmingham, U.S.A.

SUMMARY An octree-based mesh generation method is proposed to create reasonable-quality, geometry-adapted unstructured hexahedral meshes automatically from triangulated surface models without any sharp geometrical features. A new, easy-to-implement, easy-to-understand set of refinement templates is developed to perform local mesh refinement efficiently even for concave refinement domains without creating hanging nodes. A buffer layer is inserted on an octree core mesh to improve the mesh quality significantly. Laplacian-like smoothing, angle-based smoothing and local optimization-based untangling methods are used with certain restrictions to further improve the mesh quality. Several examples are shown to demonstrate the capability of our hexahedral mesh generation method for complex geometries. Copyright 2008 John Wiley & Sons, Ltd. Received 29 November 2007; Revised 1 August 2008; Accepted 16 August 2008 KEY WORDS:

unstructured hexahedral mesh; refinement templates; CSM; biomedical applications; mesh generation; grid generation

1. INTRODUCTION Computational simulations are widely used to predict or evaluate the performance of mechanical systems. The biomedical application of computational simulations aids in the understanding of complex human biological systems. One of the issues that needs to be resolved is how to represent ∗ Correspondence

to: Yasushi Ito, Enabling Technology Laboratory (ETLab), 1530 3rd Avenue South, BEC 356B, Birmingham, AL 35294-4461, U.S.A. † E-mail: [email protected] ‡ Research Assistant Professor. § Research Associate Professor, Director of ETLab. ¶ Professor and Chair. Contract/grant sponsor: Southern Consortium for Injury Biomechanics (SCIB); contract/grant number: DTNH22-01H-07551

Copyright

2008 John Wiley & Sons, Ltd.

Y. ITO, A. M. SHIH AND B. K. SONI

complex, realistic three-dimensional (3D) geometries given as biomedical image data, such as computed tomography (CT) or magnetic resonance imaging data. Patient-specific computational modeling is desired to correctly understand specific requirements. Triangular and tetrahedral meshes are widely used for computational simulations because of their flexibility for complex geometries [1]. However, triangles and tetrahedra are usually considered to be constant strain elements in finite element methods, i.e. only one value of strain is assigned for each element without any gradient. Many simplicial elements are needed to achieve the same accuracy as bilinear elements. Thus, quadrilaterals and hexahedra are preferred in computational structural mechanics (CSM) simulations. Automatic hexahedral mesh generation for arbitrary geometries is still a challenging problem. The mesh generation group at Sandia National Laboratories has been developing unconstrained paving and plastering methods [2]. They have demonstrated high-quality hexahedral mesh generation for mechanical objects. The methods, however, cannot be applied to general biomedical geometries because of their complexity. Octree-based mesh generation methods are a promising alternative to create hexahedral meshes automatically even for complex geometries. Zhang et al. proposed an isosurface extraction method to automatically generate hexahedral meshes for geometrically complex structures using an octree-based method [3] and a smoothing method using geometric flow [4]. In these approaches, a bottom-up surface topology preserving octree-based algorithm was applied to select a starting octree level. Then, a dual contouring method was used to extract a preliminary uniform hexahedral mesh. The mesh could be refined based on a feature-sensitive error function calculated from the region that users are interested in or based on computational results using refinement templates for hexahedra by Schneiders et al. [5]. Finally, a smoothing method was applied to the mesh to improve the quality. Boyd and M¨uller proposed a similar, but simpler, octree-based mesh generation method without any refinement [6]. Mesh refinement for hexahedral elements without any hanging nodes is another challenging problem. Since computational resources are always limited for simulations, meshes with as few nodes as possible are preferred. Unlike tetrahedral elements, it is not easy to maintain connectivity of hexahedral elements during the mesh refinement. Refinement template-based methods are widely used for octree-based hexahedral mesh generation to maintain the connectivity of hexahedra easily. Schneiders et al. proposed the first template-based method using four templates for node, edge, face and volume refinement as shown in Figure 1 [5]. The black points in the figure represent the positions to be refined. The first three templates produce transition elements. The volume refinement template subdivides a hexahedron into 27 similar hexahedra. In their approach, hexahedra in the region of interest are subdivided using the volume refinement template first. The neighboring hexahedra are then subdivided using all the four templates so that the resulting mesh has no hanging nodes. This subdivision is based on the nodes marked as needed to be refined, corresponded to the black points in Figure 1. If a hexahedron has a different pattern of the marked nodes, it must be converted to one of the four templates as shown in Table I. The table shows that most of the refinement patterns must be subdivided using the volume refinement template. As a result, the four templates allow the refinement of only convex domains. The entire domain or most parts of the domain sometimes need to be refined to achieve conformal refinement. A new set of templates is needed to refine meshes more efficiently. Zhang and Zhao proposed a hexahedral mesh refinement method using six templates [7]. The templates are for node, edge, face, volume, three-node (three nodes of one of the six quadrilaterals) and two-face (two quadrilaterals sharing an edge) refinement. The templates for node and volume refinement are the same as those developed by Schneiders et al. The idea of the remaining templates Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

Figure 1. Four refinement templates by Schneiders et al. for: (a) node; (b) edge; (c) face; and (d) volume refinement.

Table I. Possible refinement patterns and corresponding refinement templates by Schneiders et al.

is from conformal refinement using a shrink-and-connect strategy, known as pillowing [8–10]. Interestingly, the structure of the face refinement template is the same as the upper part of the face refinement template used by Schneiders et al. Although Zhang and Zhao showed interesting examples, their paper has three major problems. First, they did not show the internal structures of their templates. Thus, no one can reproduce their work easily. Second, not all the faces of the last two templates for three-node and two-face refinement have line symmetry based on Table I in Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

their paper; hence, each of the hexahedra assigned the three-node or two-face refinement template cannot be subdivided independently. The absence of discussion on how the neighboring hexahedra are subdivided also makes it impossible to replicate and truly understand their process. Third, it is not clear how the templates for node and three-node refinement are used. Parrish et al. proposed a selective approach to conformal refinement of existing unstructured hexahedral meshes [11]. Although their refinement approach is similar to the one by Zhang and Zhao [7], the discussion is clearer. Their refinement process consists of two methods: (A) elementby-element refinement using six templates and (B) directional refinement. Although the idea of directional refinement is from pillowing, elements are independently refined using a method, the so-called ranking system, to improve the efficiency of the algorithm. However, there is no need to apply that extra, still less efficient refinement method if the templates used in method A can be applied to concave refinement domains. This makes their refinement system difficult to implement. Moreover, the discussion is limited to one level of local refinement. If more than one level of local refinement is needed, transition elements can be refined more than once. The quality of such elements is usually quite poor. For these reasons, we propose a new, better set of templates and discuss how to achieve multiple levels of refinement without refining transition elements in this paper. The octree-based hexahedral mesh generation methods cannot create excellent-quality, reasonably-coarse hexahedral meshes, especially when mesh refinement is applied. Although it is almost impossible to create such meshes for any geometry using an octree, we need a method that produces better-quality meshes than the existing methods. There are no demonstrated methods for creating octree-based, geometry-adapted, good-quality hexahedral meshes. Excellent-quality meshes cannot be created using octrees because there are many irregular nodes on the boundary surface. Irregular nodes are defined as nodes that have too many or too few neighboring hexahedra. If interior nodes of a hexahedral mesh are considered, irregular nodes can be defined as nodes connected to more or less than six edges. The biggest disadvantage of the octree-based hexahedral mesh generation methods is the quality of resulting meshes near the boundary surfaces, where significant changes occur during computational simulations. Suppose there is a hexahedron and five of its eight nodes need to be projected to a surface in order to create a boundary conforming mesh. If the surface is flat, the shape of the hexahedron becomes quite poor. An efficient method is required to create good-quality hexahedral elements around the boundary surfaces to improve the overall mesh quality. In this paper, we propose an octree-based mesh generation method to automatically create hexahedral meshes for complex biomedical geometries with reasonable mesh quality. The starting point of the mesh generation can be either original medical image data or triangulated surface models extracted from the medical image data. Since the latter approach is more flexible to control noise in the medical image data, we use the triangulated surface models as input. A new set of refinement templates is proposed to create geometry- or solution-adaptive hexahedral meshes without producing many extra elements. To improve the final mesh quality, a buffer layer is inserted on an octree core mesh and node smoothing is applied. Although the node smoothing method is simple, it is practically important to know how it should be applied to octree-based hexahedral meshes. The remainder of this paper is organized as follows. Section 2 describes our mesh generation method. In Section 3, several biomedical geometries are shown to demonstrate and discuss the capability of our octree-based mesh generation approach. Finally, conclusions are given in Section 4. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

2. MESH GENERATION METHOD Figures 2 and 3 show the mesh generation process for a human pediatric liver model. There are four major steps: (1) the generation of a hexahedral core mesh (Figure 2(b)) inside an input triangulated surface model (Figure 2(a)) using five refinement templates (Section 2.1); (2) the removal of special hexahedra that become poor-quality elements in the later steps (Section 2.2); (3) the addition of a buffer layer (Figure 2(c)) on the boundary surface to improve the quality of the final mesh (Section 2.3); and (4) the application of node smoothing and boundary projection methods to create a boundary conforming mesh (Figure 3(a); Section 2.4). We assume in this paper that the input surface model does not have any sharp geometrical features. The recovery of topology, such as vertexes and curves, is not considered. The scaled Jacobian metric is used to check the quality of hexahedra or quadrilaterals during the mesh generation [12]. Since a hexahedron has eight nodes, the minimum of eight scaled Jacobian values, jmin, is used to check its quality. Negative Jacobian elements are defined as elements whose minimum scaled Jacobian values are less than zero. 2.1. Refinement templates The existing sets of refinement templates cannot provide flexibility for complex geometries and ease of implementation as briefly discussed in Section 1. A new set of five refinement templates is proposed here to resolve this issue as shown in Figure 4. Templates for edge (Figure 4(a)), face (Figure 4(b)) and volume (Figure 4(c)) refinement are from pillowing. Zhang and Zhao [7] and Parrish et al. [11] used the same templates as part of their template sets. In this paper, two new templates are proposed for three-node and two-face refinement. Before these new templates are introduced, let us prepare symbols. Each of the templates has eight base nodes, A, B, C, D, E, F, G and H . The orientation of the templates is shown in Figure 4. The templates for edge, face, volume, three-node and two-face refinement are represented as T2 (ABCDEFGH), T4 (ABCDEFGH), T8 (ABCDEFGH), T3 (ABCDEFGH) and T6 (ABCDEFGH), respectively. The underlines indicate nodes to be refined, which correspond to the black points with white letters in Figure 4. Figure 4(d, e) shows the structures of the two new templates and does not contain all the subdivided hexahedra to simplify the drawing. Nodes I –P in each of the new templates are added to create base hexahedra. For T3 , the hexahedron is divided into six base hexahedra first, as shown in Figure 4(d). T2 is then applied to two hexahedra: T2 (ABJIEFNM) and T2 (DAILHEMP). The other four hexahedra, CDLKGHPO, BCKJFGON, IJKLMNOP and MNOPEFGH, are used as is. For T6 , the hexahedron is divided into five base hexahedra first, as shown in Figure 4(e). T2 is applied to two hexahedra: T2 (BFGCJNOK) and T2 (HDCGPLKO). T4 is applied to two other hexahedra: T4 (AEFBIMNJ) and T4 (DHEALPMI). The hexahedron in the center, IJKLMNOP, is used as is. To understand the five refinement templates better, Figure 5 shows three types of their sides, TQ2 (ABCD), TQ3 (ABCD) and TQ4 (ABCD). Note that these are created by default if the five refinement templates for hexahedra are used. TQ3 is used for Quadrilateral ABCD in Figure 4(d) as TQ3 (ABCD) and Quadrilaterals ABCD and EFGH in Figure 4(e) as TQ3 (ABCD) and TQ3 (EFGH), respectively. TQ2 , TQ3 and TQ4 have line symmetry for easy implementation. For example, Line AC of TQ3 is the line of symmetry. This great feature allows refining hexahedral elements independently without worrying about the conformity problem discussed in [11]. Table II shows the number of hexahedra in each of the five templates. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

Figure 2. Hexahedral mesh generation process for a human pediatric liver model: (a) original triangulated surface model; (b) hexahedral bumpy core mesh; and (c) core mesh after applied smoothing and added a buffer layer. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

Figure 3. Hexahedral meshes for a human pediatric liver model: (a) final mesh after conforming to the boundary of the geometry; (b) mesh showing a cross section; (c) mesh without adaptation; and (d) mesh using the templates by Schneiders et al.

Table III shows possible refinement patterns and corresponding refinement templates using our approach (Candidate 1). A node refinement template is not used. The two new templates for threenode, T3 , and two-face refinement, T6 , allow much more flexibility for adaptive mesh generation because many of the refinement patterns assigned to T8 in [5] (Table I) are reassigned to T6 . As a result, concave domains can be refined effectively. Table IV shows another set of possible refinement patterns and corresponding refinement templates (Candidate 2). Although Candidate 2 seems to create meshes with fewer elements, Candidates 1 and 2 practically create almost the same number of elements while holding all other conditions constant. Moreover, Candidate 1 often creates the same or slightly better-quality meshes. Thus, we use Candidate 1 to assign the refinement templates to hexahedra. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

Figure 4. New set of refinement templates for: (a) edge (T2 ); (b) face (T4 ); (c) volume (T8 ); (d) three-node (T3 ); and (e) two-face refinement (T6 ). Note that (d) and (e) only show the main structures of our new templates.

Figure 5. Quadrilateral refinement templates for: (a) edge (TQ2 ); (b) three-node (TQ3 ); and (c) full refinement (TQ4 ).

Four steps are needed to create an octree core mesh using the refinement templates: (A) uniform refinement; (B) local refinement based on any sensor values for geometry-adaptive mesh generation; (C) assignment of a refinement template to each cube; and (D) removal of hexahedra outside Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

Table II. Number of hexahedra in each of the five refinement templates. Templates

T2

T4

T8

T3

T6

No. of hexahedra

5

13

27

14

37

Table III. Possible refinement patterns and corresponding refinement templates (Candidate 1).

the meshing domain. Figure 6 shows an example in 2D. Although the refinement templates can easily be applied to existing hexahedral meshes based on solution data, we do not discuss the solution-adaptive mesh generation in this paper. In Step (A), an octree is constructed from a root node—a cube that covers the entire triangulated surface model; ‘cube’ is used to represent an octree node hereafter to avoid confusion—by recursively subdividing each leaf cube into 27 sub-cubes [5]. This operation is equivalent to applying T8 to the entire domain a certain number of times. If the size of the leaf cubes becomes smaller than a user-specified value, L, the construction is stopped. In Step (B), each leaf cube that contains or intersects triangles of the initial surface model is checked. If the normals of any two of the triangles make an angle of more than !, T8 is applied Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

Table IV. Possible refinement patterns and corresponding refinement templates (Candidate 2).

to the leaf cube. ! is a user-specified parameter to detect geometrical features (or in other words, to select leaf cubes to be refined) and is used as a main sensor in this paper. The default value of ! is 60◦ . Surface curvature can be used as an additional sensor if the initial surface mesh is of good quality. l levels of local refinement are allowed. The default value of l is 1. To avoid creating low-quality elements, transition elements created by T2 , T4 , T3 and T6 are not refined when l!2. Suppose each leaf cube knows how many times it is refined, lc (0"lc "l). At node i, we can find one of the leaf cubes connecting to it that has the maximum lc . ldi = max (lc j ) c j ∈E i

(1)

where E i is the leaf cubes connecting to the node. Each leaf cube has eight or more nodes in this stage. For example, the lower left square in Figure 6(a) has 10 nodes even in 2D. At leaf cube i, the following value can be calculated: "ci = lci max −lci min lci max = max (ld j ) d j ∈Ni

(2)

lci min = min (ld j ) d j ∈Ni

Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

Figure 6. Octree core mesh generation—the numbers and shaded area represent ld and part of an input geometry, respectively: (a) uniform (bold lines) and local refinement; (b) building conforming closure; (c) assignment of the refinement templates; and (d) removal of hexahedra outside the meshing domain.

where Ni is the nodes that leaf cube i has. "ci should be 0 or 1 to apply the refinement templates. If not, the leaf cube is refined using T8 and lci +1 is assigned to the 27 new cubes as their refinement level. This process is called building conforming closure in [5] and is continued until all the leaf cubes satisfy this condition (Figure 6(b)). In Step (C), leaf cubes are subdivided using the five templates (Figure 6(c)). The pseudocode is as follows: For i = 1 to l { Select leaf cubes with lc =l −i Do { For j = 1 to 6 For k = 1 to the number of leaf cubes with lc =l −i { Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

Calculate "c of leaf cube k If "c >1 { Build conforming closure } Else if "c = 1 { Select an appropriate refinement template if cube k has one of Patterns (8− j ) to 8. }

}

} } } while any leaf cubes are modified.

where Patterns 1–8 are shown in Table III. The subdivision is first applied to leaf cubes with lc =l −1. If such a leaf cube has "c >1, the process for building conforming closure is applied. If "c = 1, an appropriate template is selected based on Table III. Black points in Table III correspond to nodes with ld =lc +1. The others correspond to nodes with ld =lc . T3 and T6 can be applied to certain refinement patterns in several ways. For example, Pattern 2B can be represented as P2B (ABCDEFGH). T3 can be applied to P2B as T3 (BCDAFGHE) or T3 (DABCHEFG), and node B or D needs to be refined accordingly. T6 can be applied to Patten 2C, P2C (ABCDEFGH ), in six ways. The six configurations are T6 (BAEFCDH G), T6 (BCDAFGH E), T6 (CBFGDAEH), T6 (DABCHEFG), T6 (EADHFBCG) and T6 (EFBAHGCD). A refinement probe, #, is used to select the best configuration. To obtain ldi (ld at node i), lc is used in Equation (1). Instead of using lc , the maximum refinement level in neighboring elements, lc max in Equation (2), can be used. ld% i = max (lc j max ) c j ∈E i

#=

!

di ∈R

ld% i

(3)

where R is nodes that may be additionally refined. The configuration with the largest # is accepted. In the P2B case, ld% B and ld% D are compared. If the refinement probe returns the same number for all the configurations of the selected pattern and template, any one of them can be selected. To avoid selecting many extra nodes to be refined, cubes with Pattern 7 or 8 are assigned T8 first, as shown in the pseudocode. Cubes with patterns that have fewer refined nodes are then assigned appropriate templates. ld at each of its eight nodes is updated if the refinement pattern of the leaf cube does not exactly match one of the five templates. For example, if T3 (BCDAFGHE) is applied to P2B , then ld at node B needs to be updated. These are repeated until all the leaf cubes with lc =l −1 have "c "1 and those with "c = 1 are assigned appropriate refinement templates. The same procedure is applied to leaf cubes with lc =l −2, lc =l −3, . . ., 0. In Step (D), hexahedra outside or intersecting the input triangulated surface model are removed to create an octree core mesh (Figures 2(b) and 6(d)). Note that leaf cubes must be small enough to capture important small features of the input geometry by specifying a small L and/or a large l. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

2.2. Removal of special hexahedra Based on the templates assigned, hexahedral elements (leaf cubes and transition elements in the refinement templates) are created inside the input triangulated surface model. To create good-quality elements in the later stages, certain elements need to be eliminated beforehand. Two important criteria are the number of quadrilaterals of a hexahedron that are on the boundary, n qh , and the number of hexahedra at a node, n hn . To identify unwanted hexahedra, edges on the boundary surface are prepared. Each of the boundary edges can have two quadrilaterals by our definition. If four quadrilaterals connect to an edge as shown in Figure 7, two edges are created, and they are considered to be duplicated. At each node on the boundary, the number of boundary edges, n en , the number of sets of duplicated boundary edges, n dn , and the number of boundary quadrilaterals, n qn , are counted. A hexahedron is removed if n qh = 4 and the hexahedron is not equilateral (i.e. a transition element in one of the refinement templates; Figure 8(a)), or if n qh = 5 (Figure 8(b)) or 6. n qh = 6 means that the hexahedron is isolated. At a node on the boundary, the hexahedra connecting to it are selected. One of them with the largest n qh is removed (A) if n hn = 2, n en !5 and n dn !1, (B) if n hn = 2 and n qn = 6, (C) if n hn = 3, n en !7 and n dn !1, (D) if n hn = 4, n en !6, n en & = 7 and n dn = 1, (H) if n hn = 5, n en !9 and n dn = 1, or (I) if n hn = 6, n en !6, n en & = 7 and n dn = 1. Figure 9 shows examples at nodes A, B, C and D. One of the hexahedra connected to node A, B or C is removed because of Condition A (n hn = 2, n en = 6 and n dn = 1) at node A, Condition B (n hn = 2 and n qn = 6) at node B, and Condition C (n hn = 3, n en = 7 and n dn = 1) at node C. The four hexahedra at node D are kept because n hn = 4, n en = 8 and n dn = 2. 2.3. Insertion of a buffer layer The next step is to insert a buffer layer on the boundary surface of the octree core mesh (Figure 10(c)). This idea comes from a paper by Mitchell and Tautges on the removal of special irregular nodes, called doublets, to improve the mesh quality significantly [8]. In 2D, a double is two quadrilaterals sharing two edges and a surface doublet is a single quadrilateral where two of whose edges are constrained to lie on a straight or near-straight curve bounding the meshed object (see Figure 10(b)). This approach suggests inserting a buffer layer under the boundary surface (Figure 10(c)), which can improve the quality of the octree-based hexahedral meshes rather than the single application of any node smoothing method. Shepherd et al. [13, 14] and Merkley et al. [15] recently applied a similar method to hexahedral meshes to improve the mesh quality. Before adding a buffer layer, a shape smoothing method is applied to the boundary surface of the octree core mesh to eliminate humps [16]. Single marching directions can be calculated at

Figure 7. Two cubes sharing edge AB. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

Figure 8. Removal of special hexahedra based on n qh : (a) two transition hexahedra with n qh = 4 and (b) hexahedron with n qh = 5.

Figure 9. Removal of special hexahedra based on n hn : bold lines represent duplicated edges.

nodes on the boundary except singular points using a standard technique used in near-field mesh generation for viscous flow simulations [17]. A singular point is a point where a single marching direction cannot be obtained without creating negative Jacobian elements. There are two types of singular points: (A) those where two local two-manifold surfaces can be defined (node D in Figure 9 and the black point in the center in Figure 11) and (B) those where they cannot be defined. Two marching directions are prepared at each of type A singular points. Examples of type B singular points are shown in [18, 19]. Since no single marching direction can be calculated at those singular points, the buffer layer is not added locally. Figures 11 and 12 show examples Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

Figure 10. Projection of a boundary surface: (a) original boundary and boundary surface of octants; (b) and (c) projection of boundary nodes without and with a buffer layer, respectively.

Figure 11. Addition of a buffer layer: suppose the black points are singular points.

of Type B singular points in 2D (imitation) and 3D, respectively, and the resulting buffer layers. Although the configuration at node A looks simple in Figure 12 (n hn = 4, n en = 6 and n qn = 6), no single marching direction can be calculated. The height of the buffer layer, h, is calculated based on L and l, h = 0.1L/3l , to avoid creating overlapping elements. The constant height is used because the quality of elements in the buffer layer is not very important at this stage unless they are negative Jacobian elements. Their quality is improved using the methods in the next section. 2.4. Smoothing and projection Smoothing methods can also improve the final mesh quality. Laplacian smoothing generally works well for octree-based hexahedral meshes as mentioned in [5]. This is sometimes true for meshes without adaptation, but not for meshes with adaptation. Negative Jacobian elements can be easily created during the smoothing process. The quality of tetrahedral meshes can be improved significantly not only with node smoothing but also with face swapping. However, face swapping is difficult for hexahedral meshes. Better smoothing methods are needed. In addition, the boundary surfaces of hexahedral meshes must be recovered at the same time. In our approach, the following set of smoothing and projection methods is applied to hexahedral meshes several times: (A) surface mesh smoothing and projection, (B) volume mesh smoothing while holding the surface fixed and (C) local surface/volume mesh smoothing without projection. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

Figure 12. Buffer layer and a Type B singular point, node A: (a) surface mesh around node A, which has four neighboring hexahedra and (b) suppressed buffer layer at node A.

Figure 13. Element quality distribution by the scaled Jacobian metric for the pediatric liver meshes.

The pseudocode for the smoothing and surface projection is as follows: For i = 1 to the number of iterations, $, { Apply Step A with each of three smoothing methods (3 iterations each) Apply Step B with each of three smoothing methods (3 iterations each) If 10.6, where %1 is the maximum principal curvature). The three levels of local refinement allow representing small features precisely, while keeping the number of elements small. The number of hexahedra with jmin"0.2 is 340. The lowest minimum scaled Jacobian is 0.0158. The smoothing method may need to be improved to create better-quality meshes with more than one level of local refinement. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

Figure 19. Hexahedral mesh for a human torso (side view): (a) surface mesh and (b) cross section through the center.

4. CONCLUSIONS An octree-based mesh generator has been developed to automatically create hexahedral meshes for complex geometries suitable for CSM simulations. Two new refinement templates are proposed in this paper and three other existing templates are used to create geometry-adaptive meshes without adding too many extra elements. The addition of a buffer layer on the boundary surface and the node smoothing methods with certain restrictions improve the resulting mesh quality significantly. The capability of our mesh generation method is demonstrated using four triangulated surface models. The resulting meshes do not have any negative Jacobian elements even for complex geometries, Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

and the quality of the meshes is acceptable. The five refinement templates allow the creation of adaptive meshes efficiently. ACKNOWLEDGEMENTS

The initial triangulated surface models of the liver and pelvis are provided courtesy of Dr King Yang and Ms. Christina Huber of Wayne State University, Detroit, MI, U.S.A. The human torso model is provided courtesy of Alexandre Olivier-Mangon and George Drettakis of INRIA by the [email protected] Shape Repository. The authors would like thank the reviewers for their constructive comments. This research is supported in part by the Southern Consortium for Injury Biomechanics (SCIB) No. DTNH22-01-H-07551. REFERENCES 1. Ito Y, Shum PC, Shih AM, Soni BK, Nakahashi K. Robust generation of high-quality unstructured meshes on realistic biomedical geometry. International Journal for Numerical Methods in Engineering 2006; 65(6):943–973. DOI: 10.1002/nme.1482. 2. Staten ML, Kerr RA, Owen SJ, Blacker TD. Unconstrained paving and plastering: progress update. Proceedings of the 15th International Meshing Roundtable, Birmingham, AL, 2006; 469–486. DOI: 10.1007/978-3-54034958-7 27. 3. Zhang Y, Bajaj C. Adaptive and quality quadrilateral/hexahedral meshing from volumetric data. Computer Methods in Applied Mechanics and Engineering 2006; 195:942–960. DOI: 10.1016/j.cma.2005.02.016. 4. Zhang Y, Bajaj C, Xu G. Surface smoothing and quality improvement of quadrilateral/hexahedral meshes with geometric flow. Proceedings of the 14th International Meshing Roundtable, San Diego, CA, 2005; 449–468. 5. Schneiders R, Schindler R, Weiler F. Octree-based generation of hexahedral element meshes. Proceedings of the 5th International Meshing Roundtable, Sandia National Laboratories, 1996; 205–216. 6. Boyd SK, M¨uller R. Smooth surface meshing for automated finite element model generation from 3D image data. Journal of Biomechanics 2006; 39:1287–1295. DOI: 10.1016/j.jbiomech.2005.03.006. 7. Zhang H, Zhao G. Adaptive hexahedral mesh generation based on local domain curvature and thickness using a modified grid-based method. Finite Elements in Analysis and Design 2007; 43:691–704. DOI: 10.1016/ j.finel.2007.03.001. 8. Mitchell SA, Tautges TJ. Pillowing doublets: refining a mesh to ensure that faces share at most one edge. Proceedings of the 4th International Meshing Roundtable, Sandia National Laboratories, 1995; 231–240. 9. Tchon K, Dompierre J, Camarero R. Automated refinement of conformal quadrilateral and hexahedral meshes. International Journal for Numerical Methods in Engineering 2004; 59:1539–1562. DOI: 10.1002/nme.926. 10. Harris N, Benzley S, Owen S. Conformal refinement of all-hexahedral element meshes based on multiple twist plane insertion. Proceedings of the 13th International Meshing Roundtable, Sandia National Laboratories, Williamsburg, VA, 2004; 157–168. SAND #2004-3765C. 11. Parrish M, Borden M, Staten M, Benzley S. A selective approach to conformal refinement of unstructured hexahedral finite element meshes. Proceedings of the 16th International Meshing Roundtable, Seattle, WA, 2007; 251–268. DOI: 10.1007/978-3-540-75103-8 15. 12. Knupp PM. Achieving finite element mesh quality via optimization of the Jacobian matrix norm and associated quantities. Part II-A framework for volume mesh optimization and the condition number of the Jacobian matrix. International Journal for Numerical Methods in Engineering 2000; 48:1165–1185. DOI: 10.1002/(SICI)10970207(20000720)48:83.0.CO;2-Y. 13. Shepherd JF. Topologic and geometric constraint-based hexahedral mesh generation. Ph.D. Dissertation, School of Computing, University of Utah, 2007. 14. Shepherd JF, Zhang Y, Tuttle CJ, Silva CT. Quality improvement and Boolean-like cutting operations in hexahedral meshes. Proceedings of the 10th International Conference on Numerical Grid Generation in Computational Field Simulations, Crete, Greece, 2007 (CD-ROM). 15. Merkley K, Ernst C, Shepherd J, Borden MJ. Methods and applications of generalized sheet insertion for hexahedral meshing. Proceedings of the 16th International Meshing Roundtable, Seattle, WA, 2007; 233–250. DOI: 10.1007/978-3-540-75103-8 14. 16. Taubin G. A signal processing approach to fair surface design. Proceedings of ACM SIGGRAPH 95, 1995; 351–358. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION 17. Kallinderis Y, Ward S. Prismatic grid generation for three-dimensional complex geometries. AIAA Journal 1993; 31(10):1850–1856. 18. Ito Y, Shih AM, Soni BK, Nakahashi K. Multiple marching direction approach to generate high quality hybrid meshes. AIAA Journal 2007; 45(1):162–167. DOI: 10.2514/1.23260. 19. Ito Y, Shih AM, Soni BK. Hybrid mesh generation with embedded surfaces. Proceedings of the 10th International Conference on Numerical Grid Generation in Computational Field Simulations, Crete, Greece, September 2007 (CD-ROM). 20. Zhou T, Shimada K. An angle-based approach to two-dimensional mesh smoothing. Proceedings of the 9th International Meshing Roundtable, New Orleans, LA, 2000; 373–384. 21. Ito Y, Nakahashi K. Improvements in the reliability and quality of unstructured hybrid mesh generation. International Journal for Numerical Methods in Fluids 2004; 45(1):79–108. DOI: 10.1002/fld.669. 22. Freitag LA, Plassmann PE. Local optimization-based untangling algorithms for quadrilateral meshes. Proceedings of the 10th International Meshing Roundtable, Newport Beach, CA, 2001; 397–406.

Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Octree-based reasonable-quality hexahedral mesh generation using a new set of refinement templates Yasushi Ito∗, †, ‡ , Alan M. Shih§ and Bharat K. Soni¶ Enabling Technology Laboratory (ETLab), Department of Mechanical Engineering, University of Alabama at Birmingham, U.S.A.

SUMMARY An octree-based mesh generation method is proposed to create reasonable-quality, geometry-adapted unstructured hexahedral meshes automatically from triangulated surface models without any sharp geometrical features. A new, easy-to-implement, easy-to-understand set of refinement templates is developed to perform local mesh refinement efficiently even for concave refinement domains without creating hanging nodes. A buffer layer is inserted on an octree core mesh to improve the mesh quality significantly. Laplacian-like smoothing, angle-based smoothing and local optimization-based untangling methods are used with certain restrictions to further improve the mesh quality. Several examples are shown to demonstrate the capability of our hexahedral mesh generation method for complex geometries. Copyright 2008 John Wiley & Sons, Ltd. Received 29 November 2007; Revised 1 August 2008; Accepted 16 August 2008 KEY WORDS:

unstructured hexahedral mesh; refinement templates; CSM; biomedical applications; mesh generation; grid generation

1. INTRODUCTION Computational simulations are widely used to predict or evaluate the performance of mechanical systems. The biomedical application of computational simulations aids in the understanding of complex human biological systems. One of the issues that needs to be resolved is how to represent ∗ Correspondence

to: Yasushi Ito, Enabling Technology Laboratory (ETLab), 1530 3rd Avenue South, BEC 356B, Birmingham, AL 35294-4461, U.S.A. † E-mail: [email protected] ‡ Research Assistant Professor. § Research Associate Professor, Director of ETLab. ¶ Professor and Chair. Contract/grant sponsor: Southern Consortium for Injury Biomechanics (SCIB); contract/grant number: DTNH22-01H-07551

Copyright

2008 John Wiley & Sons, Ltd.

Y. ITO, A. M. SHIH AND B. K. SONI

complex, realistic three-dimensional (3D) geometries given as biomedical image data, such as computed tomography (CT) or magnetic resonance imaging data. Patient-specific computational modeling is desired to correctly understand specific requirements. Triangular and tetrahedral meshes are widely used for computational simulations because of their flexibility for complex geometries [1]. However, triangles and tetrahedra are usually considered to be constant strain elements in finite element methods, i.e. only one value of strain is assigned for each element without any gradient. Many simplicial elements are needed to achieve the same accuracy as bilinear elements. Thus, quadrilaterals and hexahedra are preferred in computational structural mechanics (CSM) simulations. Automatic hexahedral mesh generation for arbitrary geometries is still a challenging problem. The mesh generation group at Sandia National Laboratories has been developing unconstrained paving and plastering methods [2]. They have demonstrated high-quality hexahedral mesh generation for mechanical objects. The methods, however, cannot be applied to general biomedical geometries because of their complexity. Octree-based mesh generation methods are a promising alternative to create hexahedral meshes automatically even for complex geometries. Zhang et al. proposed an isosurface extraction method to automatically generate hexahedral meshes for geometrically complex structures using an octree-based method [3] and a smoothing method using geometric flow [4]. In these approaches, a bottom-up surface topology preserving octree-based algorithm was applied to select a starting octree level. Then, a dual contouring method was used to extract a preliminary uniform hexahedral mesh. The mesh could be refined based on a feature-sensitive error function calculated from the region that users are interested in or based on computational results using refinement templates for hexahedra by Schneiders et al. [5]. Finally, a smoothing method was applied to the mesh to improve the quality. Boyd and M¨uller proposed a similar, but simpler, octree-based mesh generation method without any refinement [6]. Mesh refinement for hexahedral elements without any hanging nodes is another challenging problem. Since computational resources are always limited for simulations, meshes with as few nodes as possible are preferred. Unlike tetrahedral elements, it is not easy to maintain connectivity of hexahedral elements during the mesh refinement. Refinement template-based methods are widely used for octree-based hexahedral mesh generation to maintain the connectivity of hexahedra easily. Schneiders et al. proposed the first template-based method using four templates for node, edge, face and volume refinement as shown in Figure 1 [5]. The black points in the figure represent the positions to be refined. The first three templates produce transition elements. The volume refinement template subdivides a hexahedron into 27 similar hexahedra. In their approach, hexahedra in the region of interest are subdivided using the volume refinement template first. The neighboring hexahedra are then subdivided using all the four templates so that the resulting mesh has no hanging nodes. This subdivision is based on the nodes marked as needed to be refined, corresponded to the black points in Figure 1. If a hexahedron has a different pattern of the marked nodes, it must be converted to one of the four templates as shown in Table I. The table shows that most of the refinement patterns must be subdivided using the volume refinement template. As a result, the four templates allow the refinement of only convex domains. The entire domain or most parts of the domain sometimes need to be refined to achieve conformal refinement. A new set of templates is needed to refine meshes more efficiently. Zhang and Zhao proposed a hexahedral mesh refinement method using six templates [7]. The templates are for node, edge, face, volume, three-node (three nodes of one of the six quadrilaterals) and two-face (two quadrilaterals sharing an edge) refinement. The templates for node and volume refinement are the same as those developed by Schneiders et al. The idea of the remaining templates Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

Figure 1. Four refinement templates by Schneiders et al. for: (a) node; (b) edge; (c) face; and (d) volume refinement.

Table I. Possible refinement patterns and corresponding refinement templates by Schneiders et al.

is from conformal refinement using a shrink-and-connect strategy, known as pillowing [8–10]. Interestingly, the structure of the face refinement template is the same as the upper part of the face refinement template used by Schneiders et al. Although Zhang and Zhao showed interesting examples, their paper has three major problems. First, they did not show the internal structures of their templates. Thus, no one can reproduce their work easily. Second, not all the faces of the last two templates for three-node and two-face refinement have line symmetry based on Table I in Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

their paper; hence, each of the hexahedra assigned the three-node or two-face refinement template cannot be subdivided independently. The absence of discussion on how the neighboring hexahedra are subdivided also makes it impossible to replicate and truly understand their process. Third, it is not clear how the templates for node and three-node refinement are used. Parrish et al. proposed a selective approach to conformal refinement of existing unstructured hexahedral meshes [11]. Although their refinement approach is similar to the one by Zhang and Zhao [7], the discussion is clearer. Their refinement process consists of two methods: (A) elementby-element refinement using six templates and (B) directional refinement. Although the idea of directional refinement is from pillowing, elements are independently refined using a method, the so-called ranking system, to improve the efficiency of the algorithm. However, there is no need to apply that extra, still less efficient refinement method if the templates used in method A can be applied to concave refinement domains. This makes their refinement system difficult to implement. Moreover, the discussion is limited to one level of local refinement. If more than one level of local refinement is needed, transition elements can be refined more than once. The quality of such elements is usually quite poor. For these reasons, we propose a new, better set of templates and discuss how to achieve multiple levels of refinement without refining transition elements in this paper. The octree-based hexahedral mesh generation methods cannot create excellent-quality, reasonably-coarse hexahedral meshes, especially when mesh refinement is applied. Although it is almost impossible to create such meshes for any geometry using an octree, we need a method that produces better-quality meshes than the existing methods. There are no demonstrated methods for creating octree-based, geometry-adapted, good-quality hexahedral meshes. Excellent-quality meshes cannot be created using octrees because there are many irregular nodes on the boundary surface. Irregular nodes are defined as nodes that have too many or too few neighboring hexahedra. If interior nodes of a hexahedral mesh are considered, irregular nodes can be defined as nodes connected to more or less than six edges. The biggest disadvantage of the octree-based hexahedral mesh generation methods is the quality of resulting meshes near the boundary surfaces, where significant changes occur during computational simulations. Suppose there is a hexahedron and five of its eight nodes need to be projected to a surface in order to create a boundary conforming mesh. If the surface is flat, the shape of the hexahedron becomes quite poor. An efficient method is required to create good-quality hexahedral elements around the boundary surfaces to improve the overall mesh quality. In this paper, we propose an octree-based mesh generation method to automatically create hexahedral meshes for complex biomedical geometries with reasonable mesh quality. The starting point of the mesh generation can be either original medical image data or triangulated surface models extracted from the medical image data. Since the latter approach is more flexible to control noise in the medical image data, we use the triangulated surface models as input. A new set of refinement templates is proposed to create geometry- or solution-adaptive hexahedral meshes without producing many extra elements. To improve the final mesh quality, a buffer layer is inserted on an octree core mesh and node smoothing is applied. Although the node smoothing method is simple, it is practically important to know how it should be applied to octree-based hexahedral meshes. The remainder of this paper is organized as follows. Section 2 describes our mesh generation method. In Section 3, several biomedical geometries are shown to demonstrate and discuss the capability of our octree-based mesh generation approach. Finally, conclusions are given in Section 4. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

2. MESH GENERATION METHOD Figures 2 and 3 show the mesh generation process for a human pediatric liver model. There are four major steps: (1) the generation of a hexahedral core mesh (Figure 2(b)) inside an input triangulated surface model (Figure 2(a)) using five refinement templates (Section 2.1); (2) the removal of special hexahedra that become poor-quality elements in the later steps (Section 2.2); (3) the addition of a buffer layer (Figure 2(c)) on the boundary surface to improve the quality of the final mesh (Section 2.3); and (4) the application of node smoothing and boundary projection methods to create a boundary conforming mesh (Figure 3(a); Section 2.4). We assume in this paper that the input surface model does not have any sharp geometrical features. The recovery of topology, such as vertexes and curves, is not considered. The scaled Jacobian metric is used to check the quality of hexahedra or quadrilaterals during the mesh generation [12]. Since a hexahedron has eight nodes, the minimum of eight scaled Jacobian values, jmin, is used to check its quality. Negative Jacobian elements are defined as elements whose minimum scaled Jacobian values are less than zero. 2.1. Refinement templates The existing sets of refinement templates cannot provide flexibility for complex geometries and ease of implementation as briefly discussed in Section 1. A new set of five refinement templates is proposed here to resolve this issue as shown in Figure 4. Templates for edge (Figure 4(a)), face (Figure 4(b)) and volume (Figure 4(c)) refinement are from pillowing. Zhang and Zhao [7] and Parrish et al. [11] used the same templates as part of their template sets. In this paper, two new templates are proposed for three-node and two-face refinement. Before these new templates are introduced, let us prepare symbols. Each of the templates has eight base nodes, A, B, C, D, E, F, G and H . The orientation of the templates is shown in Figure 4. The templates for edge, face, volume, three-node and two-face refinement are represented as T2 (ABCDEFGH), T4 (ABCDEFGH), T8 (ABCDEFGH), T3 (ABCDEFGH) and T6 (ABCDEFGH), respectively. The underlines indicate nodes to be refined, which correspond to the black points with white letters in Figure 4. Figure 4(d, e) shows the structures of the two new templates and does not contain all the subdivided hexahedra to simplify the drawing. Nodes I –P in each of the new templates are added to create base hexahedra. For T3 , the hexahedron is divided into six base hexahedra first, as shown in Figure 4(d). T2 is then applied to two hexahedra: T2 (ABJIEFNM) and T2 (DAILHEMP). The other four hexahedra, CDLKGHPO, BCKJFGON, IJKLMNOP and MNOPEFGH, are used as is. For T6 , the hexahedron is divided into five base hexahedra first, as shown in Figure 4(e). T2 is applied to two hexahedra: T2 (BFGCJNOK) and T2 (HDCGPLKO). T4 is applied to two other hexahedra: T4 (AEFBIMNJ) and T4 (DHEALPMI). The hexahedron in the center, IJKLMNOP, is used as is. To understand the five refinement templates better, Figure 5 shows three types of their sides, TQ2 (ABCD), TQ3 (ABCD) and TQ4 (ABCD). Note that these are created by default if the five refinement templates for hexahedra are used. TQ3 is used for Quadrilateral ABCD in Figure 4(d) as TQ3 (ABCD) and Quadrilaterals ABCD and EFGH in Figure 4(e) as TQ3 (ABCD) and TQ3 (EFGH), respectively. TQ2 , TQ3 and TQ4 have line symmetry for easy implementation. For example, Line AC of TQ3 is the line of symmetry. This great feature allows refining hexahedral elements independently without worrying about the conformity problem discussed in [11]. Table II shows the number of hexahedra in each of the five templates. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

Figure 2. Hexahedral mesh generation process for a human pediatric liver model: (a) original triangulated surface model; (b) hexahedral bumpy core mesh; and (c) core mesh after applied smoothing and added a buffer layer. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

Figure 3. Hexahedral meshes for a human pediatric liver model: (a) final mesh after conforming to the boundary of the geometry; (b) mesh showing a cross section; (c) mesh without adaptation; and (d) mesh using the templates by Schneiders et al.

Table III shows possible refinement patterns and corresponding refinement templates using our approach (Candidate 1). A node refinement template is not used. The two new templates for threenode, T3 , and two-face refinement, T6 , allow much more flexibility for adaptive mesh generation because many of the refinement patterns assigned to T8 in [5] (Table I) are reassigned to T6 . As a result, concave domains can be refined effectively. Table IV shows another set of possible refinement patterns and corresponding refinement templates (Candidate 2). Although Candidate 2 seems to create meshes with fewer elements, Candidates 1 and 2 practically create almost the same number of elements while holding all other conditions constant. Moreover, Candidate 1 often creates the same or slightly better-quality meshes. Thus, we use Candidate 1 to assign the refinement templates to hexahedra. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

Figure 4. New set of refinement templates for: (a) edge (T2 ); (b) face (T4 ); (c) volume (T8 ); (d) three-node (T3 ); and (e) two-face refinement (T6 ). Note that (d) and (e) only show the main structures of our new templates.

Figure 5. Quadrilateral refinement templates for: (a) edge (TQ2 ); (b) three-node (TQ3 ); and (c) full refinement (TQ4 ).

Four steps are needed to create an octree core mesh using the refinement templates: (A) uniform refinement; (B) local refinement based on any sensor values for geometry-adaptive mesh generation; (C) assignment of a refinement template to each cube; and (D) removal of hexahedra outside Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

Table II. Number of hexahedra in each of the five refinement templates. Templates

T2

T4

T8

T3

T6

No. of hexahedra

5

13

27

14

37

Table III. Possible refinement patterns and corresponding refinement templates (Candidate 1).

the meshing domain. Figure 6 shows an example in 2D. Although the refinement templates can easily be applied to existing hexahedral meshes based on solution data, we do not discuss the solution-adaptive mesh generation in this paper. In Step (A), an octree is constructed from a root node—a cube that covers the entire triangulated surface model; ‘cube’ is used to represent an octree node hereafter to avoid confusion—by recursively subdividing each leaf cube into 27 sub-cubes [5]. This operation is equivalent to applying T8 to the entire domain a certain number of times. If the size of the leaf cubes becomes smaller than a user-specified value, L, the construction is stopped. In Step (B), each leaf cube that contains or intersects triangles of the initial surface model is checked. If the normals of any two of the triangles make an angle of more than !, T8 is applied Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

Table IV. Possible refinement patterns and corresponding refinement templates (Candidate 2).

to the leaf cube. ! is a user-specified parameter to detect geometrical features (or in other words, to select leaf cubes to be refined) and is used as a main sensor in this paper. The default value of ! is 60◦ . Surface curvature can be used as an additional sensor if the initial surface mesh is of good quality. l levels of local refinement are allowed. The default value of l is 1. To avoid creating low-quality elements, transition elements created by T2 , T4 , T3 and T6 are not refined when l!2. Suppose each leaf cube knows how many times it is refined, lc (0"lc "l). At node i, we can find one of the leaf cubes connecting to it that has the maximum lc . ldi = max (lc j ) c j ∈E i

(1)

where E i is the leaf cubes connecting to the node. Each leaf cube has eight or more nodes in this stage. For example, the lower left square in Figure 6(a) has 10 nodes even in 2D. At leaf cube i, the following value can be calculated: "ci = lci max −lci min lci max = max (ld j ) d j ∈Ni

(2)

lci min = min (ld j ) d j ∈Ni

Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

Figure 6. Octree core mesh generation—the numbers and shaded area represent ld and part of an input geometry, respectively: (a) uniform (bold lines) and local refinement; (b) building conforming closure; (c) assignment of the refinement templates; and (d) removal of hexahedra outside the meshing domain.

where Ni is the nodes that leaf cube i has. "ci should be 0 or 1 to apply the refinement templates. If not, the leaf cube is refined using T8 and lci +1 is assigned to the 27 new cubes as their refinement level. This process is called building conforming closure in [5] and is continued until all the leaf cubes satisfy this condition (Figure 6(b)). In Step (C), leaf cubes are subdivided using the five templates (Figure 6(c)). The pseudocode is as follows: For i = 1 to l { Select leaf cubes with lc =l −i Do { For j = 1 to 6 For k = 1 to the number of leaf cubes with lc =l −i { Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

Calculate "c of leaf cube k If "c >1 { Build conforming closure } Else if "c = 1 { Select an appropriate refinement template if cube k has one of Patterns (8− j ) to 8. }

}

} } } while any leaf cubes are modified.

where Patterns 1–8 are shown in Table III. The subdivision is first applied to leaf cubes with lc =l −1. If such a leaf cube has "c >1, the process for building conforming closure is applied. If "c = 1, an appropriate template is selected based on Table III. Black points in Table III correspond to nodes with ld =lc +1. The others correspond to nodes with ld =lc . T3 and T6 can be applied to certain refinement patterns in several ways. For example, Pattern 2B can be represented as P2B (ABCDEFGH). T3 can be applied to P2B as T3 (BCDAFGHE) or T3 (DABCHEFG), and node B or D needs to be refined accordingly. T6 can be applied to Patten 2C, P2C (ABCDEFGH ), in six ways. The six configurations are T6 (BAEFCDH G), T6 (BCDAFGH E), T6 (CBFGDAEH), T6 (DABCHEFG), T6 (EADHFBCG) and T6 (EFBAHGCD). A refinement probe, #, is used to select the best configuration. To obtain ldi (ld at node i), lc is used in Equation (1). Instead of using lc , the maximum refinement level in neighboring elements, lc max in Equation (2), can be used. ld% i = max (lc j max ) c j ∈E i

#=

!

di ∈R

ld% i

(3)

where R is nodes that may be additionally refined. The configuration with the largest # is accepted. In the P2B case, ld% B and ld% D are compared. If the refinement probe returns the same number for all the configurations of the selected pattern and template, any one of them can be selected. To avoid selecting many extra nodes to be refined, cubes with Pattern 7 or 8 are assigned T8 first, as shown in the pseudocode. Cubes with patterns that have fewer refined nodes are then assigned appropriate templates. ld at each of its eight nodes is updated if the refinement pattern of the leaf cube does not exactly match one of the five templates. For example, if T3 (BCDAFGHE) is applied to P2B , then ld at node B needs to be updated. These are repeated until all the leaf cubes with lc =l −1 have "c "1 and those with "c = 1 are assigned appropriate refinement templates. The same procedure is applied to leaf cubes with lc =l −2, lc =l −3, . . ., 0. In Step (D), hexahedra outside or intersecting the input triangulated surface model are removed to create an octree core mesh (Figures 2(b) and 6(d)). Note that leaf cubes must be small enough to capture important small features of the input geometry by specifying a small L and/or a large l. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

2.2. Removal of special hexahedra Based on the templates assigned, hexahedral elements (leaf cubes and transition elements in the refinement templates) are created inside the input triangulated surface model. To create good-quality elements in the later stages, certain elements need to be eliminated beforehand. Two important criteria are the number of quadrilaterals of a hexahedron that are on the boundary, n qh , and the number of hexahedra at a node, n hn . To identify unwanted hexahedra, edges on the boundary surface are prepared. Each of the boundary edges can have two quadrilaterals by our definition. If four quadrilaterals connect to an edge as shown in Figure 7, two edges are created, and they are considered to be duplicated. At each node on the boundary, the number of boundary edges, n en , the number of sets of duplicated boundary edges, n dn , and the number of boundary quadrilaterals, n qn , are counted. A hexahedron is removed if n qh = 4 and the hexahedron is not equilateral (i.e. a transition element in one of the refinement templates; Figure 8(a)), or if n qh = 5 (Figure 8(b)) or 6. n qh = 6 means that the hexahedron is isolated. At a node on the boundary, the hexahedra connecting to it are selected. One of them with the largest n qh is removed (A) if n hn = 2, n en !5 and n dn !1, (B) if n hn = 2 and n qn = 6, (C) if n hn = 3, n en !7 and n dn !1, (D) if n hn = 4, n en !6, n en & = 7 and n dn = 1, (H) if n hn = 5, n en !9 and n dn = 1, or (I) if n hn = 6, n en !6, n en & = 7 and n dn = 1. Figure 9 shows examples at nodes A, B, C and D. One of the hexahedra connected to node A, B or C is removed because of Condition A (n hn = 2, n en = 6 and n dn = 1) at node A, Condition B (n hn = 2 and n qn = 6) at node B, and Condition C (n hn = 3, n en = 7 and n dn = 1) at node C. The four hexahedra at node D are kept because n hn = 4, n en = 8 and n dn = 2. 2.3. Insertion of a buffer layer The next step is to insert a buffer layer on the boundary surface of the octree core mesh (Figure 10(c)). This idea comes from a paper by Mitchell and Tautges on the removal of special irregular nodes, called doublets, to improve the mesh quality significantly [8]. In 2D, a double is two quadrilaterals sharing two edges and a surface doublet is a single quadrilateral where two of whose edges are constrained to lie on a straight or near-straight curve bounding the meshed object (see Figure 10(b)). This approach suggests inserting a buffer layer under the boundary surface (Figure 10(c)), which can improve the quality of the octree-based hexahedral meshes rather than the single application of any node smoothing method. Shepherd et al. [13, 14] and Merkley et al. [15] recently applied a similar method to hexahedral meshes to improve the mesh quality. Before adding a buffer layer, a shape smoothing method is applied to the boundary surface of the octree core mesh to eliminate humps [16]. Single marching directions can be calculated at

Figure 7. Two cubes sharing edge AB. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

Figure 8. Removal of special hexahedra based on n qh : (a) two transition hexahedra with n qh = 4 and (b) hexahedron with n qh = 5.

Figure 9. Removal of special hexahedra based on n hn : bold lines represent duplicated edges.

nodes on the boundary except singular points using a standard technique used in near-field mesh generation for viscous flow simulations [17]. A singular point is a point where a single marching direction cannot be obtained without creating negative Jacobian elements. There are two types of singular points: (A) those where two local two-manifold surfaces can be defined (node D in Figure 9 and the black point in the center in Figure 11) and (B) those where they cannot be defined. Two marching directions are prepared at each of type A singular points. Examples of type B singular points are shown in [18, 19]. Since no single marching direction can be calculated at those singular points, the buffer layer is not added locally. Figures 11 and 12 show examples Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

Figure 10. Projection of a boundary surface: (a) original boundary and boundary surface of octants; (b) and (c) projection of boundary nodes without and with a buffer layer, respectively.

Figure 11. Addition of a buffer layer: suppose the black points are singular points.

of Type B singular points in 2D (imitation) and 3D, respectively, and the resulting buffer layers. Although the configuration at node A looks simple in Figure 12 (n hn = 4, n en = 6 and n qn = 6), no single marching direction can be calculated. The height of the buffer layer, h, is calculated based on L and l, h = 0.1L/3l , to avoid creating overlapping elements. The constant height is used because the quality of elements in the buffer layer is not very important at this stage unless they are negative Jacobian elements. Their quality is improved using the methods in the next section. 2.4. Smoothing and projection Smoothing methods can also improve the final mesh quality. Laplacian smoothing generally works well for octree-based hexahedral meshes as mentioned in [5]. This is sometimes true for meshes without adaptation, but not for meshes with adaptation. Negative Jacobian elements can be easily created during the smoothing process. The quality of tetrahedral meshes can be improved significantly not only with node smoothing but also with face swapping. However, face swapping is difficult for hexahedral meshes. Better smoothing methods are needed. In addition, the boundary surfaces of hexahedral meshes must be recovered at the same time. In our approach, the following set of smoothing and projection methods is applied to hexahedral meshes several times: (A) surface mesh smoothing and projection, (B) volume mesh smoothing while holding the surface fixed and (C) local surface/volume mesh smoothing without projection. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

Figure 12. Buffer layer and a Type B singular point, node A: (a) surface mesh around node A, which has four neighboring hexahedra and (b) suppressed buffer layer at node A.

Figure 13. Element quality distribution by the scaled Jacobian metric for the pediatric liver meshes.

The pseudocode for the smoothing and surface projection is as follows: For i = 1 to the number of iterations, $, { Apply Step A with each of three smoothing methods (3 iterations each) Apply Step B with each of three smoothing methods (3 iterations each) If 10.6, where %1 is the maximum principal curvature). The three levels of local refinement allow representing small features precisely, while keeping the number of elements small. The number of hexahedra with jmin"0.2 is 340. The lowest minimum scaled Jacobian is 0.0158. The smoothing method may need to be improved to create better-quality meshes with more than one level of local refinement. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION

Figure 19. Hexahedral mesh for a human torso (side view): (a) surface mesh and (b) cross section through the center.

4. CONCLUSIONS An octree-based mesh generator has been developed to automatically create hexahedral meshes for complex geometries suitable for CSM simulations. Two new refinement templates are proposed in this paper and three other existing templates are used to create geometry-adaptive meshes without adding too many extra elements. The addition of a buffer layer on the boundary surface and the node smoothing methods with certain restrictions improve the resulting mesh quality significantly. The capability of our mesh generation method is demonstrated using four triangulated surface models. The resulting meshes do not have any negative Jacobian elements even for complex geometries, Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

Y. ITO, A. M. SHIH AND B. K. SONI

and the quality of the meshes is acceptable. The five refinement templates allow the creation of adaptive meshes efficiently. ACKNOWLEDGEMENTS

The initial triangulated surface models of the liver and pelvis are provided courtesy of Dr King Yang and Ms. Christina Huber of Wayne State University, Detroit, MI, U.S.A. The human torso model is provided courtesy of Alexandre Olivier-Mangon and George Drettakis of INRIA by the [email protected] Shape Repository. The authors would like thank the reviewers for their constructive comments. This research is supported in part by the Southern Consortium for Injury Biomechanics (SCIB) No. DTNH22-01-H-07551. REFERENCES 1. Ito Y, Shum PC, Shih AM, Soni BK, Nakahashi K. Robust generation of high-quality unstructured meshes on realistic biomedical geometry. International Journal for Numerical Methods in Engineering 2006; 65(6):943–973. DOI: 10.1002/nme.1482. 2. Staten ML, Kerr RA, Owen SJ, Blacker TD. Unconstrained paving and plastering: progress update. Proceedings of the 15th International Meshing Roundtable, Birmingham, AL, 2006; 469–486. DOI: 10.1007/978-3-54034958-7 27. 3. Zhang Y, Bajaj C. Adaptive and quality quadrilateral/hexahedral meshing from volumetric data. Computer Methods in Applied Mechanics and Engineering 2006; 195:942–960. DOI: 10.1016/j.cma.2005.02.016. 4. Zhang Y, Bajaj C, Xu G. Surface smoothing and quality improvement of quadrilateral/hexahedral meshes with geometric flow. Proceedings of the 14th International Meshing Roundtable, San Diego, CA, 2005; 449–468. 5. Schneiders R, Schindler R, Weiler F. Octree-based generation of hexahedral element meshes. Proceedings of the 5th International Meshing Roundtable, Sandia National Laboratories, 1996; 205–216. 6. Boyd SK, M¨uller R. Smooth surface meshing for automated finite element model generation from 3D image data. Journal of Biomechanics 2006; 39:1287–1295. DOI: 10.1016/j.jbiomech.2005.03.006. 7. Zhang H, Zhao G. Adaptive hexahedral mesh generation based on local domain curvature and thickness using a modified grid-based method. Finite Elements in Analysis and Design 2007; 43:691–704. DOI: 10.1016/ j.finel.2007.03.001. 8. Mitchell SA, Tautges TJ. Pillowing doublets: refining a mesh to ensure that faces share at most one edge. Proceedings of the 4th International Meshing Roundtable, Sandia National Laboratories, 1995; 231–240. 9. Tchon K, Dompierre J, Camarero R. Automated refinement of conformal quadrilateral and hexahedral meshes. International Journal for Numerical Methods in Engineering 2004; 59:1539–1562. DOI: 10.1002/nme.926. 10. Harris N, Benzley S, Owen S. Conformal refinement of all-hexahedral element meshes based on multiple twist plane insertion. Proceedings of the 13th International Meshing Roundtable, Sandia National Laboratories, Williamsburg, VA, 2004; 157–168. SAND #2004-3765C. 11. Parrish M, Borden M, Staten M, Benzley S. A selective approach to conformal refinement of unstructured hexahedral finite element meshes. Proceedings of the 16th International Meshing Roundtable, Seattle, WA, 2007; 251–268. DOI: 10.1007/978-3-540-75103-8 15. 12. Knupp PM. Achieving finite element mesh quality via optimization of the Jacobian matrix norm and associated quantities. Part II-A framework for volume mesh optimization and the condition number of the Jacobian matrix. International Journal for Numerical Methods in Engineering 2000; 48:1165–1185. DOI: 10.1002/(SICI)10970207(20000720)48:83.0.CO;2-Y. 13. Shepherd JF. Topologic and geometric constraint-based hexahedral mesh generation. Ph.D. Dissertation, School of Computing, University of Utah, 2007. 14. Shepherd JF, Zhang Y, Tuttle CJ, Silva CT. Quality improvement and Boolean-like cutting operations in hexahedral meshes. Proceedings of the 10th International Conference on Numerical Grid Generation in Computational Field Simulations, Crete, Greece, 2007 (CD-ROM). 15. Merkley K, Ernst C, Shepherd J, Borden MJ. Methods and applications of generalized sheet insertion for hexahedral meshing. Proceedings of the 16th International Meshing Roundtable, Seattle, WA, 2007; 233–250. DOI: 10.1007/978-3-540-75103-8 14. 16. Taubin G. A signal processing approach to fair surface design. Proceedings of ACM SIGGRAPH 95, 1995; 351–358. Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme

OCTREE-BASED REASONABLE-QUALITY HEXAHEDRAL MESH GENERATION 17. Kallinderis Y, Ward S. Prismatic grid generation for three-dimensional complex geometries. AIAA Journal 1993; 31(10):1850–1856. 18. Ito Y, Shih AM, Soni BK, Nakahashi K. Multiple marching direction approach to generate high quality hybrid meshes. AIAA Journal 2007; 45(1):162–167. DOI: 10.2514/1.23260. 19. Ito Y, Shih AM, Soni BK. Hybrid mesh generation with embedded surfaces. Proceedings of the 10th International Conference on Numerical Grid Generation in Computational Field Simulations, Crete, Greece, September 2007 (CD-ROM). 20. Zhou T, Shimada K. An angle-based approach to two-dimensional mesh smoothing. Proceedings of the 9th International Meshing Roundtable, New Orleans, LA, 2000; 373–384. 21. Ito Y, Nakahashi K. Improvements in the reliability and quality of unstructured hybrid mesh generation. International Journal for Numerical Methods in Fluids 2004; 45(1):79–108. DOI: 10.1002/fld.669. 22. Freitag LA, Plassmann PE. Local optimization-based untangling algorithms for quadrilateral meshes. Proceedings of the 10th International Meshing Roundtable, Newport Beach, CA, 2001; 397–406.

Copyright

2008 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2008) DOI: 10.1002/nme