- Email: [email protected]
Resolving topology ambiguity for multiple-material domains
Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
Contents lists available at SciVerse ScienceDirect
Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma
Resolving topology ambiguity for multiple-material domains Yongjie Zhang ⇑, Jin Qian Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, PA 15213, USA
a r t i c l e
i n f o
Article history: Received 1 April 2012 Received in revised form 18 July 2012 Accepted 31 July 2012 Available online 24 August 2012 Keywords: Dual contouring Topology ambiguity Multiple-material domain Hybrid octree Triangular mesh Tetrahedral mesh
a b s t r a c t This paper describes a novel approach to resolve topology ambiguity in generating quality triangular and tetrahedral meshes for multiple-material domains. In previous works, we developed an octree-based dual contouring method to construct surface and volumetric meshes for complicated domains with topology ambiguity resolved. However, this method is based on enumeration and cannot handle domains with multiple materials. As a follow up, in this study we develop an automatic and efficient approach to resolve topology ambiguity for multiple-material domains. We first use an indicator variable to identify the topology inside each octree leaf cell; whenever one cell is ambiguous, we split the cell into 12 tetrahedra. In the following we analyze each material-change edge to create triangles or boundary tetrahedra with surrounding minimizers. For tetrahedral meshes, in addition to material-change edges we need to analyze each interior edge to create interior tetrahedra. With the ambiguity-free mesh ready, we are able to apply quality improvement techniques to attain good mesh quality. Our main contribution lies in resolving topology ambiguity by analyzing a hybrid octree with both cubic and tetrahedral leaf cells, and this method works for both homogeneous and multiple-material domains. We have applied our algorithm to three complicated material datasets and obtained good results. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction With the development of scanned techniques and super computers, the capability of collecting and storing large datasets has been advanced rapidly, and the scanned domain such as polycrystalline materials is usually heterogeneous and with complicated topology. In recent years, various ways have been developed to model these multiple-material domains [30,36,10]. However in comprehensive modeling and prediction of material behavior, there still lacks automatic and robust geometric modeling and mesh generation techniques, which are able to generate quality meshes for all material regions with correct topology. Resolving topology ambiguity is a challenging problem especially for multiple-material domains. In previous works, we have developed an octree-based Dual Contouring (DC) method to construct triangular and tetrahedral meshes for homogeneous domains with correct topology [33], by categorizing all possible ambiguity cases into 31 groups and modifying the dual mesh accordingly based on enumeration. However, this method only works for domains with one single material, and cannot handle complicated multiple-material domains like the polycrystalline material in Fig. 1. To address this deficiency, in this study we develop a novel approach to automatically and efficiently handle topology ambiguity for multiple material domains. ⇑ Corresponding author. Tel.: +1 412 268 5332; fax: +1 412 268 3348. E-mail address: [email protected] (Y. Zhang). 0045-7825/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2012.07.022
To distinguish topology ambiguities, we need to know the function property inside each cubic cell. Unlike single-material domains, where we can use a tri-linear interpolation function to model the interior domain, there is no such modeling function for multiple material cells. Instead, here we choose an indicator variable to examine each composing material and see whether topology ambiguity arises. Once one ambiguity is found, we split the cubic cell into 12 tetrahedral cells. The tessellation by tetrahedra guarantees to eliminate ambiguities in the resulting mesh. Till now we are able to modify the standard DC method based on such a hybrid octree consisting of both cubic and tetrahedral leaf cells. For each leaf cell, we insert one minimizer by either minimizing a pre-defined error function or simply choosing the cell center. Then, we classify all edges into cell edges, internal edges and face diagonals. In the following, we analyze each material-change edge, whose two end points lie in two different materials, and use its surrounding minimizers to form a polygon in the triangular mesh, or a polyhedron in the tetrahedral mesh. In the case of tetrahedralization, we also need to study each interior edge, whose two end points lie in the same material, to construct interior tetrahedra. Till now a mesh with no topology ambiguity is formed. However, we still need to improve the mesh quality via a combination of vertex classification, face swapping, edge removal and geometric flowbased smoothing. Our algorithm generates a conformal mesh and attains good quality with all topology ambiguity resolved. We have applied our algorithm to several complicated material datasets and achieved good results.
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
167
Fig. 1. (a) Adaptive triangular and tetrahedral meshes of a 92-grain polycrystalline material (multi-material domain). (b) Adaptive mesh between two interior grains; (c) Zoom-in picture of the planar surface. (d, e) zoom-in details of an interior local region with topology ambiguity (the red circle). The green line in (d) denotes one non-manifold edge, which is resolved in (e). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
The remainder of this paper is organized as follows: Section 2 reviews related previous work. Section 3 talks about the DC methods. Sections 4 and 5 discuss resolving 2D and 3D ambiguities. Section 6 presents how to improve the mesh quality. Section 7 illustrates results, and finally Section 8 draws conclusions. 2. Previous work 2.1. Isocontouring Isocontouring is the process of generating a piece-wise linear approximation to a given zero-surface of an implicit function. Of the various contouring methods, the Marching Cubes (MC) and DC are of most importance. The MC technique [15] is able to generate a closed, manifold triangular mesh for any scalar domain. To capture sharp features, additional information such as surface normals are used to handle sampling alias [8]. In addition, various techniques were developed to triangulate over cells with different levels. However, there are cracks when the resolution levels of adjacent cells are different, and a fan of triangles needs to be inserted around the gravity center of the coarse triangles [29]. Moreover, people developed ways to generate tetrahedral meshes involved with the body centered cubic lattice [9]. Topology ambiguity issues were awared of and studied [1,18,14]. Differently, the DC technique [7,24] uses an adaptive octree and produces water-tight contours over adaptive octree grids. However, it fails to recognize the ambiguous situations in the octree cells.
be in two separate pieces or there may be a single ‘‘tunnel’’ piece [18,1]. These ambiguities can be tracked using the collapsing rule [5], which repeatedly collapses the edges of the cell whose end points lie on the same side of the iso-surface. If the collapsing result is not an edge or point, the cell is ambiguous. To study possible topologies, the interior part of the cube is usually modeled as trilinear interpolations. After discussing the property of such tri-linear interpolation function, we obtain a classification of 31 cases [1]. With such classification, multiple attempts to solve the ambiguity problem were made [18,2,14]. These techniques are all based on the MC algorithm, and they model the iso-surface cell by cell. These methods provide a comprehensive solution to the ambiguity problem, however, they only work for single-material domains. 2.3. Image-based modeling for microstructures Image-based modeling is an important tool to study the properties and to predict the response of microstructures. With the development of automated collection and advanced image processing techniques, characterization of microstructure materials becomes tractable [6,22,26]. The collected datasets often provide statistically important representations of the microstructure. Moreover, they can be incorporated into the modeling and simulation for the design and prediction of material behaviors, which has been a cutting-edge research area in recent years [12,17,16]. People have developed various ways to model the microstructures [30,36,10]. However, an efficient and robust meshing technique that can handle topology ambiguity is still lacking.
2.2. Resolving topology ambiguity 3. A review of dual contouring methods Topology ambiguities arise when there are multiple choices of topology configurations for the given iso-surface [28]. It includes face ambiguity, when the four edges of a cube face intersect the iso-surface [20]; and body ambiguity, where the iso-surface may
Since our algorithm is based on the DC method, let us first review it. DC [7] is a grid-based method for contouring a scalar field. With a uniform grid, the memory consumption on complex models
168
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
with small features can be very high. Therefore we choose an adaptive octree, which has denser cells near the boundary or feature, and coarser cells inside the volume. There are three steps in our algorithm: construct an adaptive octree, compute a Quadratic Error Function (QEF) [3] for each leaf, and then generate dual meshes for this octree. The adaptive octree is achieved by subdividing the volume into octree cells recursively until the desired surface accuracy is achieved. The subdivision procedure begins with one octree cell that contains the whole domain, and it ends when the octree has achieved a certain adaptation level. The octree adaptation can be controlled with a feature sensitive function [34], which measures the iso-contour difference between two neighboring octree levels. The feature sensitive function is defined as P F ¼ ðjf iþ1 f i j=jrf i jÞ, where f i is a tri-linear interpolation function inside the cell at Level i [34]. The subdivision is carried out until the error is less than a pre-defined threshold. In order to obtain triangles with good quality, we restrict the neighboring octree level difference to be less than or equal to two. The next step is to construct a predefined QEF and to obtain a minimizer for each leaf cell by minimizing this function. The QEF is defined P as QEFðxÞ ¼ ðni ðx pi ÞÞ2 , where pi and ni are the position and unit normal vectors at the intersection points between the isosurface and octree cell edges [3,4]. The final step is to analyze each sign-change edge in manifold domains, which is defined as one edge whose two end points lie on different sides of the isocontour. Since the adaptive octree consists of leaves with different resolutions, each leaf may have neighbors at different levels. In order to avoid generating redundant elements, we always choose the minimal edge that does not contain an edge of the neighboring leaf. In this adaptive octree, the minimal edge is shared by either four or three cells, and we connect the minimizers of these cells to obtain a hybrid mesh consisting of quadrilateral (quad) and triangular elements. Finally these quads are divided to form an adaptive triangular mesh. This DC method has been extended to tetrahedral [34] and hexahedral [32] mesh generation. In tetrahedral mesh, interior edges, whose two end points lie inside the volume, need to be analyzed in addition to sign-change edges. For each sign-change edge, its three or four surrounding minimizers and its interior end point construct a tetrahedron or a pyramid. For each interior edge, its three or four surrounding minimizers and its two interior end points construct a diamond or pyramid. These pyramids and diamonds are then split into tetrahedra. In hexahedral meshes, we analyze each interior grid point instead. In a uniform case, each grid point is shared by eight octree cells and we construct a hexahedron from the minimizers of these cells. Several variants of the DC method have been developed as well, which introduce better topology and mesh quality [19,23,24,31]. Parallel algorithms were also introduced so that DC can be performed directly on modern Graphics Processing Unit architectures [25]. Furthermore, the DC method has also been extended to meshing a domain with multiple materials [36]. Based on the DC methodology, the material-change edge is used to identify the interface between material regions. In the meshing process, both materialchange edges and interior edges are analyzed to construct tetrahedral meshes, and interior grid points are analyzed for hexahedral mesh construction. Despite being adaptive and feature-preserving, one drawback of the above DC methods is that it fails to handle ambiguity in the topology. We have developed a technique that is able to distinguish and resolve all topology ambiguities in homogeneous materials with enumeration [33]. However, resolving topology ambiguity in multiple-material domains still remains an open problem. We will address this issue in the following sections.
4. Resolving topology ambiguity in 2D Ambiguity arises when the topology in each cell is not uniquely defined. In homogeneous domains, the input data contains a scalar field. If the function value at a grid point is greater than a given isovalue, then it has a positive sign; otherwise it is assigned negative. In multiple-material domains, the input data contains the material ID for each grid point. These values form the sign/material configuration of each cell. Topology ambiguity is one situation where the same sign/material configuration may lead to different topologies in representing the iso-surface. One way for checking is to use the collapsing rule [5], which repeatedly collapses the edges of one quad/ cubic cell, whose corners have the same sign/material ID, to a single vertex. The quad contains no ambiguity if and only if the result of this reduction is a single edge or a single vertex. To resolve the ambiguity, we need to distinguish the ambiguous situations and then make a conformal mesh accordingly. In this section, we will discuss 2D ambiguity problems, and then talk about complex 3D problems in the next section. 4.1. Identifying 2D ambiguity To begin with, we now study the 2D ambiguities. Since the material IDs are only given at the four corners of each quadtree cell, here we choose an indicator variable to decide the material ID of interior nodes. In one quad cell, suppose there are n different materials, and their corresponding material IDs are ID1 ; ID2 ; . . . ; IDn , where nranges from 1 to 4. For each of the material i (i = 1–n) and grid point j (j = 0–3) of the quad cell, we assign the following indicator variable:
vi;j ¼
1 If IDi ¼ material ID of Grid Point j 0
Otherwise
To illustrate, we use the example in Fig. 2. In this quad cell, there are two materials, red (ID1 ) and green (ID2 ), so n = 2. For material ID1 , grid points P 1 and P 3 belong to it; therefore their indicator variables for ID1 are 1. On the other hand, the indicator variables of grid points P0 and P2 are 0, as shown in (b). Similarly, for material ID2 , grid points P0 and P2 belong to it; therefore their indicator variables for ID2 are 1. On the other hand, the indicator variables of grid points P1 and P3 are 0, as shown in (c). For each material i (i = 1 to n), we choose a bi-linear interpolation, Fðn; gÞ, to represent the indicator variable in the interior domain. The bi-linear interpolation can be written in terms of the function values at four corners F ij ði; j ¼ 0; 1Þ:
Fðn; gÞ ¼ F 00 ð1 nÞð1 gÞ þ F 01 ð1 nÞg þ F 10 nð1 gÞ þ F 11 ng:
ð1Þ
Given an interior point P, suppose there are n materials in the cell, and the material ID is simply the material ID whose indicator variable is the maximum. Depending on the material ID of the quad center, we summarize all the 2D configurations into 8 unique cases with the consideration of symmetry and complementary. In order to differentiate distinct situations, we use the labelling notations as follows: the first number denotes the number of materials in the cube; the second label distinguishes the material ID of the four nodes; and the third label denotes the resolving of topology ambiguities. The quad cell may only contain one material (Case 1); two materials (Cases 2-1, 2-2-1, 2-2-2); or three materials (Cases 3-1, 3-2-1, and 3-22); otherwise the cell has four distinct materials (Case 4). From Fig. 3, we can tell that ambiguity arises when one of its materials has topology ambiguity, i.e, when one of its materials does not satisfy the collapsing rule. To be specific, there are ambiguities in Cases 2-2-1, 2-2-2, 3-2-1, 3-2-2, where two nodes of the
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
169
Fig. 2. An example of indicator variables: (a) contains two materials and their materials IDs are ID1 (red), ID2 (green); (b-c) show the indicator variables of material ID1 and ID2 , respectively, where a solid circle denotes vi;j ¼ 1, and an empty circle denotes vi;j ¼ 0. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
same material occupy the cell diagonal. After resolving these ambiguities, there are a total of 8 cases for a 2D cell (see Fig. 3). Simple configurations: Cases 1, 2-1, 3-1, and 4 have no ambiguity. These cases can be meshed using the standard DC. Configuration 2-2: There are two possible cases: Case 2-2-1: Green grid points are joined inside the cell; and the green material spreads across the cell; Case 2-2-2: Green grid points are separated; and the green isosurface breaks into two pieces inside this cell. On the other hand, red grid points are joined. Configuration 3-2: There are two possible cases: Case 3-2-1: Red grid points are joined inside the cell; and the red material spreads across the cell; Case 3-2-2: Red grid points are separated; and the green isosurface breaks into two pieces inside this cell. 4.2. Handling 2D ambiguity in triangulation After distinguishing topology ambiguities, now we mesh the 2D multiple-material domain. For those simple cases, the standard DC is sufficient to represent the correct topology. However, in order to handle the ambiguous cases, modifications to the standard DC is necessary. As discussed above, ambiguities are due to the alternative choice of topology for a quad cell containing two opposite corners with the same material ID. In fact, these ambiguities can be removed entirely by tessellating the ambiguous quad cell with triangles [28].
The disadvantage of this scheme is that it creates more triangles in the final mesh. However, considering the small portion of ambiguous cells in the quadtree, the element number increase is negligible. Following the above principle, for those ambiguous cells we choose a quad diagonal along which the material ID does not change, and use this diagonal to split the quad into two triangles. If no such diagonal exists (Case 3-2-2), we choose either one of the two diagonals. These triangles are treated as ‘‘triangular cells’’, and therefore the original quadtree has become a hybrid tree consisting of quad and triangular cells. Till now the tessellation is done. Now, we construct meshes with respect to the hybrid tree. Material-change edges and interior edges are analyzed to construct triangular meshes for each material region. Material-change edge: Each material-change edge is shared by two boundary cells in a uniform case, and one minimizer is computed for each of them. Note that these cells may be quads or triangles. These two minimizers and each end point of the material-change edge construct a triangle. If one material is the background, then no triangle is created from it. Therefore, a total of one or two triangles are obtained. In adaptive case, we always consider the minimal edges. With such restriction, each material-change edge is also shared by two boundary cells, but these two cells may have different sizes. From these two cells we can also get two minimizers to construct two triangles.
Fig. 3. 2D iso-surface topology possibilities. The cell has up to four materials. The iso-surface (red, green, blue and yellow) intersects the cell. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
170
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
Fig. 4. The meshing of 2D topology ambiguities. The dash lines form quad/triangular cells, the grey nodes are minimizers, while the colored nodes are grid points. In (a, b), the red nodes are connected in the middle while in (c, d) the red nodes are separated. In (a) and (c) the red, blue and green color represents the material ID of each triangle; in (b) and (d) the color denotes the type of edges around which triangles are built. The yellow triangles are generated by analyzing material-change edges; and the tan triangles are generated by interior edges (cell diagonal). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 5. An adaptive local region with ambiguity. The rendering color represents different materials. (a) The red material spreads the central domain; and (b) the red domain is separated in the central domain. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Interior edge: Each interior edge is shared by two quad or triangle cells. For a boundary cell, we determine the minimizer using the QEF. Otherwise the center of the cell is chosen as the minimizer. This interior edge and each of the minimizer point construct a triangle, therefore two triangles are obtained. Uniform and adaptive grids only differ in cell size, and the same meshing method can be applied. We demonstrate the resolving and meshing of topology ambiguities for Case 3-2-1 (a, b) and Case 3-2-2 (c, d) in Fig. 4. Note in (a) and (c) the rendering color represents the material ID while in (b) and (d) the color denotes the edge type. Fig. 5 shows the meshing example of an adaptive ambiguous region. To make the drawing clear, we do not draw all the triangle edges. 5. Resolving topology ambiguity in 3D Resolving 3D topology ambiguity is much more complicated. Any grid point can belong to eight possible materials, yielding a total of 88 cases. Even after considering symmetry and complementary, there are still more than 300 unique cases, making enumeration impracticable. Therefore, a general approach is needed to solve this problem automatically without considering too much tedious details. Here we first introduce the categorization of topology ambiguities, then present how to resolve them, and finally discuss how to integrate these topologies into mesh generation. Topology ambiguity is caused by tessellating the space using cubes. There can be multiple choices of topologies for a cube with any face or body diagonal containing the same material ID. Again we can use the same collapsing rule [5] to find out the topology ambiguities. We can repeatedly collapse edges of the cubic cell whose two end points have the same material ID to a single vertex. If the result of this reduction turns out to be an edge or a vertex, the cubic cell is free of ambiguity.
Topology ambiguity can be categorized into two groups: face ambiguity and interior ambiguity. Face ambiguity arises when two grid points with the same material ID occupy one cubic face diagonal, as shown in Fig. 6(a) and (b). The red grid points belong to Material 1, while the green ones belong to Material 2. The four edges of this face are all material-change edges, and we need to decide how to triangulate this local region. Meanwhile, there are additional interior ambiguities inside the cubic cell, as shown in Fig. 6(c) and (d). In these situations the iso-surface of Material 1 can either be two separated surfaces or a tunnel piece passing through the cubic cell. These interior structures decide what the topology looks like inside this cell. 5.1. Identifying 3D ambiguity In 3D cells, the material IDs are only given at the eight corner points. In order to resolve ambiguity, we need to know the material ID of interior points. Once again we use an indicator variable. In one cubic cell, suppose there are n different materials, and their corresponding material IDs are ID1 ; ID2 ; . . . ; IDn , where n ranges from 1 to 8. For each of the material i (i = 1 to n) and grid point j (j = 0 to 7) of the cubic cell, we assign the following indicator variable:
vi;j ¼
1 If IDi ¼ material ID of Grid Point j 0
Otherwise
To illustrate, we use the example in Fig. 7. In this cubic cell, there are three materials, red (ID1 ), green (ID2 ) and blue (ID3 ), so n = 3. For material ID1 , grid points P 0 ; P3 and P5 belong to it; therefore their indicator variables for ID1 are 1. On the other hand, the indicator variables of all the other grid points are 0, as shown in (b). Similarly, for material ID2 , grid points P2 and P 6 belong to it; therefore their indicator variables for ID2 are 1. On the other hand, the indicator variables of all the other grid points are 0, as shown in (c). Finally, for material ID3 , grid points P1 P4 and P7 belong to it; therefore
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
171
Fig. 6. Examples of 3D face and interior ambiguities. The cubic cell has two grid points (red) belonging to Material 1 and all the other grid points (green) belonging to Material 2. The iso-surface (red) intersects with this cell. (a, b) One face ambiguity; and (c, d) one interior ambiguity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 7. An example of indicator variables: (a) contains three materials and their material IDs are ID1 (red), ID2 (green) and ID3 (blue); (b–d) show the indicator variables of material material ID1 ; ID2 and ID3 , respectively, where a solid circle denotes vi;j ¼ 1, and an empty circle denotes vi;j ¼ 0. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
their indicator variables are 1. On the other hand, the indicator variables of all the other grid points are 0, as shown in (d). The indicator variable is tri-linearly interpolated in the interior domain. The tri-linear interpolation, Fðn; g; fÞ, is written in terms of the function values at the eight grid points, F ijk ði; j; k ¼ 0; 1Þ:
Fðn; g; fÞ ¼ F 000 ð1 nÞð1 gÞð1 fÞ þ F 001 ð1 nÞð1 gÞf þ F 010 ð1 nÞgð1 fÞ þ F 011 ð1 nÞgf þ F 100 nð1 gÞð1 fÞ þ F 101 nð1 gÞf þ F 110 ngð1 fÞ þ F 111 ngf:
ð2Þ
Given an interior point P, suppose there are n materials in the cell, the material ID is simply the material ID whose indicator variable is the maximum. For each material the cell contains, we study the topology classification according to the criteria used in the singlematerial domain [33]. The only difference is that instead of using the function value of the face/body saddle points, we now study the material IDs at the face and body centers of each cell. If any of the materials in the cube has topology ambiguity, the cell is considered ambiguous. With the topology classification known, we are able to generate the dual mesh and incorporate the specific topology into mesh generation. 5.2. Handling 3D ambiguity in tetrahedralization It is obvious that resolving 3D ambiguity is much more complicated. Since enumeration is impossible, we develop a general approach similar to 2D domains. This approach involves with splitting ambiguous cells and mesh generation. In fact, topology ambiguities can be removed entirely by tessellating the space using tetrahedra rather than cubes [28,27]. There are 24 ¼ 16 cases and after considering symmetry there are only five unique ones, see Fig. 8. The main disadvantage of this tessellating scheme is that more elements are generated in the final mesh.
To solve this problem, we only split the ambiguous cubic cells into tetrahedra. Because the number of ambiguous cells is small compared with the total cell number, this increase can be neglected. Based on this idea, we propose the following four steps to handle topology ambiguities. These four steps include: (1) find the ambiguous octree cells and split each of them into 12 tetrahedra, forming a hybrid tree consisting of cubes and tetrahedra; (2) calculate one minimizer for each cube and choose the center of each tetrahedron as the minimizer; (3) analyze each material-change edge and interior edge to build a volumetric mesh. Details of each step are described as follows. 5.2.1. Splitting ambiguous cells To eliminate the ambiguities, we first split each ambiguous cell into tetrahedra. Moreover, each of the six faces is split into two triangles using the same rule as in 2D ambiguous faces. The cell center and each of these triangles construct a tetrahedron. This scheme is able to resolve all the topology ambiguities. Face ambiguities are resolved automatically because the 2D scheme applied on cubic faces distinguishes 2D ambiguities. The body ambiguities are also resolved because the material ID of the body center decides either a ‘‘tunnel’’ or two separated pieces exist. To illustrate, we use the cube in Fig. 9 to show the splitting schemes. For better visualization, only the front face is split. If the front face belongs to 2D Case 3-2-1 (see Fig. 3), then we choose the diagonal shown in (a). Otherwise we choose the diagonal shown in (b). 5.2.2. Computing one minimizer for each hybrid cell With all topology ambiguities identified, we are able to compute one minimizer for each of the hybrid cells. There are two types of octree cells: boundary cells and interior cells. Boundary cells include cubic cells and tetrahedral cells. For cubic and tetrahedral cells, we follow the standard DC method to compute one
172
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
Fig. 8. Five possible configurations of one tetrahedral cell.
minimizers to form a volumetric mesh. For each edge type, we have different treatments as follows.
Fig. 9. Cube splitting schemes when red nodes are joined (a) or separated (b) on the front face.
minimizer by minimizing the predefined QEF. For other interior cells we also choose the center as the minimizer. Our approach assumes that only one minimizer is generated in each cell, and since all topology ambiguities are resolved beforehand, our approach is free of any topological mis-matches. 5.2.3. Building tetrahedral mesh by analyzing edges To construct tetrahedral meshes from the hybrid octree, we need to analyze each material-change edge as well as interior edges. All the material-change edges and interior edges belong to the cubic cells or the tetrahedral cells, and they can be cell edges, internal edges and face diagonals, as shown in Fig. 10. Cell edges (Fig. 10(a)) are surrounded by a combination of cubic and tetrahedral cells. There can be 4 to 8 surrounding minimizers, corresponding to the situations where the edge is shared by 4 cubic cells, and where the edge is shared by 8 tetrahedral cells, respectively. The internal edges (Fig. 10(b)) are those which connect the center to cell corners. These edges are surrounded by 3 to 6 tetrahedra. Face diagonals (Fig. 10(c)) can be shared by either 1 cubic cell and 2 tetrahedral cells, or 4 tetrahedral cells. Therefore there can be 3 or 4 surrounding minimizers. For adaptive octree, we follow the standard DC method and always use the minimal edges. We analyze each material-change edge and interior edge in the analysis domain, connect the QEF
Material-change edge: Each material-change edge is shared by several cubic or tetrahedral cells, and one minimizer is computed for each cell. We use these minimizers to form a convex polygon. Depending on the position of the edge end points, we now have the following two tetrahedralization schemes. (1) Material-change edge with one end point in the background material: For these edges, we first triangulate the convex polygon, and construct one tetrahedron by connecting each triangle with the edge end point not in the background, as shown in Fig. 11(a). The nodes V 1 V 3 V 2 V 0 form one tetrahedron, and V 1 V 4 V 3 V 0 form the other. (2) Other material-change edges: For these edges, we construct a polyhedron using the polygon and two edge end points, and then split it into tetrahedra. Note that the tetrahedra in both sides of the polygon belong to different materials. If we choose the tetrahedralization scheme in Fig. 11(a), we may experience inverted elements in adaptive octree. Therefore we split the polyhedron by inserting one edge center node, as shown in Fig. 11(b). We construct a tetrahedron using this center node and each triangle in the polyhedral surface. In Fig. 11(b), the node V 6 is the edge center node, and we connect it with every polyhedral surface triangle, for instance, triangle V 1 V 2 V 5 to construct a tetrahedron. Interior edges: For these edges, we construct a polyhedron using the polygon and two edge end points, and then split it into tetrahedra. Again in order to avoid inverted elements, we construct the tetrahedra using the two edge end points and one edge of the polygon. In Fig. 11(c), four tetrahedra are created, and they are V 5 V 0 V 2 V 1 , V 5 V 0 V 3 V 2 , V 5 V 0 V 4 V 3 and V 5 V 0 V 1 V 4 . Fig. 12 shows the resulting meshes of five local regions before and after resolving topology ambiguity. To show the connecting
Fig. 10. Three kinds of edges (red): cell edges (a), internal edges (b) and face diagonals (c). Only the adjacent tetrahedra are drawn. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
173
Fig. 11. The tetrahedralization of three types of edges: material-change edges with one end point lying in background (a), other material-change edges (b), and interior edges (c). Different meshing schemes are shown, and the material ID of each tetrahedron is rendered in color. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 12. Meshes of five regions before (first row) and after (second row) the re-connection. The red circles indicate the locations of the ambiguity, green lines or the blue dot indicates non-manifold situations, and red dashed lines denote the axes of a tunnel topology. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
schemes more clearly, we only show the mesh of one material. The green lines or the blue node denotes non-manifold situations, where one edge is shared by more than two triangles on the surface, or one vertex whose neighborhood does not form a disc. Note in (c), along the red dashed axes, a hollow tunnel topology is missing, and our re-connection algorithm recovers it. In these five figures, it is obvious that the topology ambiguities are resolved and appropriate meshes are built. 6. Quality improvement Till now the generated mesh has correct topology and the element quality is good especially at the core; however, it may have poor-quality elements around the material boundary that influence the convergence and stability of the numerical simulation. Therefore, we need to take an additional step to improve the mesh quality. Our quality improvement scheme includes vertex classification, face swapping, edge removal, and geometric flow-based smoothing. To begin with, we need metrics to measure the quality of the triangular and tetrahedral meshes. For triangular meshes, the metric we choose is the aspect ratio, defined as 2 r in =r out , where rin is the radius of the inscribed circle, and rout is the radius of the circumcircle. For tetrahedral meshes, we choose the edge-ratio, a minimum volume bound and the Joe-Liu parameter [13] to measure the mesh quality. The edge-ratio is the ratio of the longest
edge length over the shortest edge length of one tetrahedron; the minimum volume bound is the minimum volume of all the tetrahedra in the mesh; the Joe-Liu parameter is defined as follows:
24=3 3 jsj2=3 Fðs; dÞ ¼ P ; 2 06i
ð3Þ
where jsj denotes the volume of one tetrahedron, V i (i = 0 to 3) are the four vertices, and eij denotes the edge connecting V i to V j . With these measures, the mesh quality can be judged by observing the worst element quality, and the distribution of elements in terms of their quality values. For microstructures, the datasets are defined in a cubic representative volume element (RVE). We need to preserve the RVE boundary during quality improvement. According to the relative positions to the RVE and to each grain, all the vertices can be categorized into seven distinct groups, as shown in Fig. 13 and each group needs different improving techniques. Group 1: The eight corners of the RVE box, which are fixed during the improvement in order to keep the RVE boundary. Group 2: Vertices on the twelve edges of the RVE box, which only move along the edge. If they lie inside one grain (2a), they are smoothed only along the edge; if they are shared by two or more than two grains (2b), they are fixed during the improvement. Group 3: Vertices on the six faces of the RVE box, which are smoothed only on the plane. If they are shared by more than two grains (3a), they are planar non-manifold points and are fixed dur-
174
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
Fig. 13. Vertex classification. (a) Groups 1–3 on the RVE boundary; and (b) Groups 4–7 inside the RVE volume.
ing the improvement; If they are on the planar curve shared by two grains (3b), they are smoothed along the tangent direction of the planar curve during the improvement; finally if they are the vertices of one grain (3c), they are only moved on the plane. Group 4: Vertices inside one grain, which will be improved using the weighted-averaging method. These vertices do not exist in triangular meshes. Group 5: Vertices located on the grain boundary surface patches shared by two grains, which are smoothed on the tangent plane of the grain boundary during the improvement. Group 6: Vertices located on interior curves, which can only move along the tangent direction of the interior curve. Group 7: Interior non-manifold vertices shared by more than two interior curves, which are fixed during the improvement. During the improvement process, the vertex classification and the corresponding moving method need to be followed. For instance, when we relocate a vertex in Group 3, the movement should be projected to the RVE plane. For planar or interior curves, we use a curve fairing method based on the mean curvature flow to improve them. We fix the starting and ending non-manifold points (Groups 3a, 7), and move the interior vertices (Groups 3b, 6) along the tangent direction using a B-spline interpolation and re-sampling technique [21]. For vertices on grain boundary surface patches (Group 5), we improve them using a surface diffusion flow developed in [35]. Finally for the interior vertices (Group 4), we relocate them to the mass center of their neighboring elements. In many cases, elements with bad aspect ratios may still remain in the mesh after smoothing. Therefore, we consider modifying the mesh topology and use it together with smoothing to improve the quality. In adaptive meshes, there might be vertices embedded in a triangle with an element valence number of 3 (triangular mesh), or
a tetrahedron with an element valence of 4 (tetrahedral mesh). Such vertices are first removed. In addition, if one tetrahedron has two edges lying on the same boundary curve, these two edges form a big angle and the local quality is usually bad. To handle them, we change the local connectivities to make sure that each tetrahedron has at most one edge on the boundary curve. Furthermore, face swapping and edge removal are also used to improve the local mesh topology [11]. The face swapping operation swaps the faces of several neighboring tetrahedra. There are four basic templates including 2–3 flip, 3–2 flip, 2–2 flip and 4–4 flip, as shown in Fig. 14. One swapping operation is legal only when it does not generate inverted elements. The algorithm begins by checking the Joe-Liu parameter of all the tetrahedra. When the quality is less than a given threshold, we perform the face swapping technique. For each tetrahedron, we need to check all its 4 faces and 6 edges. The following procedures are taken: (1) Check 4 faces: We examine each face of the tetrahedron. Whenever a 2–3 flip is legal, we save its resulting local qualities. (2) Check 6 edges: We analyze each edge of the tetrahedron. Depending on the element valence number of that edge, we can use 3–2 flip (the valence number equals to 3), 2–2 flip (the valence number equals to 2 and the edge lies on the RVE boundary) or 4–4 flip (the valence number equals to 4). Once a flip is legal, we also save its resulting local qualities. (3) Carry out the flip: After checking all its faces and edges, we pick the one that leads to the best local quality and carry out the face swapping operation according to the templates.
Fig. 14. Face swapping operations: (a, b) 2–3 and 3–2 flip; (b) 2–2 flip; and (c) 4–4 flip.
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
175
Fig. 15. An edge removal operation. (a) The original polyhedron with the red center edge; (b) the center edge and the surrounding minimizers; and (c) the resulting mesh after edge removal. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 16. Zoom-in details of a 308-grain polycrystalline material mesh (multi-material domain) before (a) and after (b) the quality improvement.
Edge swapping is another efficient approach in eliminating sliver tetrahedral elements. We use an edge removal operation to swap interior edges. These edges are shared by multiple tetrahedra. In Fig. 15, we show the removal of one red edge. In this process the middle polygon needs to be triangulated in the most optimized way. Note that to preserve boundary features of the multiple-material domain, the face swapping and edge removal operations cannot be carried out across grain boundaries. Fig. 16 shows the effects of our quality improvement. 7. Results and discussion Our developed algorithm is automatic and robust in resolving topology ambiguity for both homogeneous and multiple-material domains. We have applied our algorithm to two multiple-material datasets (Figs. 1 and 18) as well as one homogeneous dataset (Fig. 17). The computations were carried out on a PC with Intel Q9600 CPU and 4 GB DDR-II memories. Tables 1 and 2 show statistics of the generated triangular and tetrahedral meshes. Fig. 19 shows the histogram of aspect ratios for triangular meshes and Joe-Liu parameters for tetrahedral meshes. The minimum– maximum octree level of the Foam, 92-grain and 308-grain models are 5–8, 4–7 and 4–7, respectively. The mesh size is increased after resolving the ambiguities, since one cubic leaf cell may be split into 12 tetrahedral cells. The mesh quality is improved greatly after quality improvement. Note that the computation time includes the time for the mesh generation, resolving topology ambiguities, together with quality improvement, which takes most of the time. Our technique combines the idea of DC and Marching Tetrahedra. The spatial grid is a hybrid tree structure consisting of cubic and tet-
rahedral cells. Since the percentage of ambiguous cells is small (usually 1 in 1000 cells), such combination solves the ambiguity problem without increasing the number of vertices too much. Furthermore, compared with the topology preservation technique developed for homogeneous materials [36], which categorizes all topologies into 31 cases, our new technique has two advantages. First, it is able to handle both homogeneous and multiple-material domains. In this algorithm, the homogeneous material is treated as one special case of multiple-material domains. Second, instead of based on enumeration, our new algorithm is more general. This makes it much easier to implement and extend to tetrahedral mesh generation. For our octree-based method, it comes with several inherent limitations. Bad quality elements are generated along the boundary although good quality elements are present at the core, and the obtained mesh quality may be inferior to Delaunay meshes. For complex models with small local features, the memory consumption can be very high. Therefore it is critical to choose an adaptive octree and generate adaptive meshes. During the octree subdivision, the maximum depth is determined by surface accuracy. The maximum number of materials is limited to four for 2D domains and eight for 3D domains. This can be overcome by refining the cells locally. 8. Conclusion In this study, we have developed an efficient and automatic algorithm to resolve topology ambiguity in generating triangular and tetrahedral meshes for multiple-material domains. We begin with an indicator variable to decide the topology inside each octree leaf cell. When one ambiguity is detected, we split the cell into 12
176
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
Fig. 17. (a) Adaptive triangular and tetrahedral meshes of a foam (homogeneous) material. (b) shows one local region of the tetrahedral mesh with some elements removed; (c) shows a zoom-in picture of one local boundary; (d) and (e) show zoom-in details of an interior local region with topology ambiguity (the red circle). The green line in (d) denotes one non-manifold edge, which is resolved in (e). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 18. (a) Adaptive triangular and tetrahedral meshes of a 308-grain polycrystalline material (multi-material domains). (b) shows the adaptive mesh between two interior grains; (c) shows a zoom-in picture of the planar surface. (d) and (e) show zoom-in details of an interior local region with topology ambiguity (the red circle). The green line in (d) denotes one non-manifold edge, which is resolved in (e). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
177
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178 Table 1 Statistics of the generated triangle meshes. Model
Mesh
Mesh size (Vert #, Elem #)
Aspect ratio (Worst, Best)
Ambiguity Cell #
Time (s)
Foam (Fig. 17)
Original Improved
(188,352, 378,755) (189,423, 380,412)
(0.01, 1.0) (0.10, 1.0)
214 214
135 161
92 Grain data (Fig. 1)
Original Improved
(15,789, 34,918) (16,000, 36,286)
(0.01, 1.0) (0.13, 1.0)
30 30
14 21
308 Grain data (Fig. 18)
Original Improved
(17,092, 40,715) (18,713, 42,092)
(0.04, 1.0) (0.14, 1.0)
14 14
23 30
Table 2 Statistics of the generated tetrahedral meshes. Model
Mesh
Mesh size (Vert #, Elem #)
Edge ratio (Worst, Best)
Volume (Min, Max)
Joe–Liu parameter (Worst, Best)
Time (s)
Foam (Fig. 17)
Original Improved
(553,612, 2,712,452) (554,706, 2,721,620)
(411.5, 1.03) (11.50, 1.01)
(1.0e5, 0.87) (6.2e2, 0.67)
(2.0e6, 1.0) (0.1, 1.0)
258 395
92 Grain data (Fig. 1)
Original Improved
(76,765, 407,525) (77,333, 408,664)
(312.32, 1.01) (13.21, 1.0)
(1.2e4, 1.52) (2.7e-3, 1.20)
(2.5e4, 0.98) (0.10, 1.0)
37 101
308 Grain data (Fig. 18)
Original Improved
(92,110, 494,455) (92,769, 495,353)
(631.44, 1.01) (13.12, 1.02)
(1.0e6, 1.30) (4.6e3, 0.84)
(1.0e6, 1.0) (0.10, 1.0)
65 123
Fig. 19. The histogram of aspect ratios for triangular meshes (a) and Joe-Liu Parameters for tetrahedral meshes (b).
tetrahedra. In the following we study each material-change edge, and use its surrounding minimizers to create triangles or boundary tetrahedra. In addition, we analyze each interior edge to create interior tetrahedra. Finally, we apply quality improvement techniques to attain a quality mesh. In this work we build up a hybrid octree and generate quality meshes with all topology ambiguities resolved. We have applied our algorithm to three complicated datasets and obtained good results. In the future we would like to test more multiple material datasets in various application fields. Acknowledgement The research was supported in part by an AFOSR grant FA955011-1-0346 and a NSF/DoD-MRSEC seed grant. References [1] E. Chernyaev, Marching cubes 33: construction of topologically correct isosurfaces, Technical Report CH/95-17, CERN, 1995. [2] P. Cignoni, F. Ganovelli, C. Montani, R. Scopigno, Reconstruction of topologically correct and adaptive trilinear isosurfaces, Computers and Graphics 24 (2000) 399–418.
[3] M. Garland, P. Heckbert, Simplifying surfaces with color and texture using quadric error metrics, in: IEEE Visualization’98, 1998, pp. 263–270. [4] M. Garland, E. Shaffer, A multiphase approach to efficient surface simplification, in: IEEE Visualization’02, 2002, pp. 117–124. [5] T. Gerstner, R. Paarola, Topology preserving and controlled topology simplifying multiresolution isosurface extraction, in: IEEE Visualization’00, 2000, pp. 259–266. [6] M.A. Groeber, B.K. Haley, M.D. Uchic, D.M. Dimiduk, S. Ghosh, 3D reconstruction and characterization of polycrystalline microstructures using a FIBCSEM system, Materials Characterization 57 (4–5) (2006) 399–418. [7] T. Ju, F. Losaasso, S. Schaefer, J. Warren, Dual contouring of Hermite data, ACM Transactions on Graphics 21 (2002) 339–346. [8] L. Kobbelt, M. Botsch, U. Schwanecke, Feature sensitive surface extraction from volume data, in: SIGGRAPH’01, 2001, pp. 57–66. [9] F. Labelle, J. Shewchuk, Isosurface stuffing: fast tetrahedral meshes with good dihedral angles, ACM Transactions on Graphics 26 (3) (2007) 57.1– 57.10. [10] S.A. Langer, E.R. Fuller Jr., W.C. Carter, OOF: an image-based finite-element analysis of material microstructures, Computing in Science Engineering 3 (3) (2001) 15–23. [11] J. Leng, Y. Zhang, G. Xu, A novel geometric flow-driven approach for quality improvement of segmented tetrahedral meshes, in: International Meshing Roundtable, 2011, pp. 347–364. [12] A.C. Lewis, A.B. Geltmacher, Image-based modeling of the response of experimental 3D microstructures to mechanical loading, Scripta Materialia 55 (1) (2006) 81–85. [13] A. Liu, B. Joe, Relationship between tetrahedron shape measures, BIT Numerical Mathematics 34 (2) (1994) 268–287.
178
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
[14] A. Lopes, K. Brodlie, Improving the robustness and accuracy of the marching cubes algorithm for isosurfacing, IEEE Transactions on Visualization and Computer Graphics 9 (1) (2003) 16–29. [15] W. Lorensen, H. Cline, Marching cubes: a high resolution 3D surface construction algorithm, in: SIGGRAPH’87, vol. 21, 1987, pp. 163–169. [16] J. MacSleynea, M.D. Uchicb, J.P. Simmonsb, M. De Graefa, Three-dimensional analysis of secondary precipitates in René-88 dt and UMF-20 superalloys, Acta Materialia 57 (20) (2009) 6251–6267. [17] J. Madisona, J. Spowartb, D. Rowenhorstc, L.K. Aagesend, K. Thorntona, T.M. Pollocka, Modeling fluid flow in three-dimensional single crystal dendritic structures, Acta Materialia 58 (8) (2010) 2864–2875. [18] B. Natarajan, On generating topologically correct isosurfaces from uniform samples, The Visual Computer 11 (1994) 52–62. [19] G. Nielson, Dual marching cubes, in: IEEE Visualization’04, 2004, pp. 489– 496. [20] G. Nielson, B. Hamann, The asymptotic decider: resolving the ambiguity in marching cubes, in: IEEE Visualization’92, 1992, pp. 83–91. [21] J. Qian, Y. Zhang, W. Wang, A.C. Lewis, M.A.S. Qidwai, A.B. Geltmacher, Quality improvement of non-manifold hexahedral meshes for critical feature determination of microstructure materials, International Journal for Numerical Methods in Engineering 82 (11) (2010) 1406–1423. [22] D.J. Rowenhorst, A.C. Lewis, G. Spanos, Three-dimensional analysis of grain topology and interface curvature in a b-titanium alloy, Acta Materialia 58 (16) (2010) 5511–5519. [23] S. Schaefer, J. Warren, Dual marching cubes: primal contouring of dual grids, in: Computer Graphics and Applications’04, 2004, pp. 70–76. [24] S. Schaefer, T. Ju, J. Warren, Manifold dual contouring, IEEE Transactions on Visualization and Computer Graphics 13 (3) (2007) 610–619. [25] L. Schmitz, L.F. Scheidegger, D.K. Osmari, C.A. Dietrich, J.L.D. Comba, Efficient and quality contouring algorithms on the GPU, Computer Graphics Forum 29 (8) (2010) 2569–2578.
[26] Z. Shan, A.M. Gokhale, Micromechanics of complex three-dimensional microstructures, Acta Materialia 49 (11) (2001) 2001–2015. [27] B.S. Sohn, Topology preserving tetrahedral decomposition applied to trilinear interval volume tetrahedrization, KSII Transaction one Internet and Information Systems 3 (6) (2009) 667–681. [28] G.M. Treece, R.W. Prager, A.H. Gee, Regularised marching tetrahedra: improved iso-surface extraction, Computers and Graphics 23 (1998) 583–598. [29] J.R. Westermann, L. Kobbelt, T. Ertl, Real-time exploration of regular volume data by adaptive reconstruction of isosurfaces, The Visual Computer 15 (2) (1999) 100–111. [30] Z. Wu, J.M. Sullivan Jr., Multiple material marching cubes algorithm, International Journal for Numerical Methods in Engineering 58 (2) (2003) 189–207. [31] N. Zhang, W. Hong, A. Kaufman, Dual contouring with topology preserving simplification using enhanced cell representation, in: IEEE Visualization’04, 2004, pp. 505–512. [32] Y. Zhang, C. Bajaj, Adaptive and quality quadrilateral/hexahedral meshing from volumetric data, Computer Methods in Applied Mechanics and Engineering 195 (9–12) (2006) 942–960. [33] Y. Zhang, J. Qian, Dual contouring for domains with topology ambiguity, Computer Methods in Applied Mechanics and Engineering 217–220 (2012) 34–45. [34] Y. Zhang, C. Bajaj, B.-S. Sohn, 3D finite element meshing from imaging data, Computer Methods in Applied Mechanics and Engineering 194 (2005) 5083– 5106. [35] Y. Zhang, C. Bajaj, G. Xu, Surface smoothing and quality improvement of quadrilateral/hexahedral meshes with geometric flow, Communications in Numerical Methods in Engineering 25 (1) (2009) 1–18. [36] Y. Zhang, T. Hughes, C. Bajaj, An automatic 3D mesh generation method for domains with multiple materials, Computer Methods in Applied Mechanics and Engineering 199 (5–8) (2010) 405–415.
Contents lists available at SciVerse ScienceDirect
Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma
Resolving topology ambiguity for multiple-material domains Yongjie Zhang ⇑, Jin Qian Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, PA 15213, USA
a r t i c l e
i n f o
Article history: Received 1 April 2012 Received in revised form 18 July 2012 Accepted 31 July 2012 Available online 24 August 2012 Keywords: Dual contouring Topology ambiguity Multiple-material domain Hybrid octree Triangular mesh Tetrahedral mesh
a b s t r a c t This paper describes a novel approach to resolve topology ambiguity in generating quality triangular and tetrahedral meshes for multiple-material domains. In previous works, we developed an octree-based dual contouring method to construct surface and volumetric meshes for complicated domains with topology ambiguity resolved. However, this method is based on enumeration and cannot handle domains with multiple materials. As a follow up, in this study we develop an automatic and efficient approach to resolve topology ambiguity for multiple-material domains. We first use an indicator variable to identify the topology inside each octree leaf cell; whenever one cell is ambiguous, we split the cell into 12 tetrahedra. In the following we analyze each material-change edge to create triangles or boundary tetrahedra with surrounding minimizers. For tetrahedral meshes, in addition to material-change edges we need to analyze each interior edge to create interior tetrahedra. With the ambiguity-free mesh ready, we are able to apply quality improvement techniques to attain good mesh quality. Our main contribution lies in resolving topology ambiguity by analyzing a hybrid octree with both cubic and tetrahedral leaf cells, and this method works for both homogeneous and multiple-material domains. We have applied our algorithm to three complicated material datasets and obtained good results. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction With the development of scanned techniques and super computers, the capability of collecting and storing large datasets has been advanced rapidly, and the scanned domain such as polycrystalline materials is usually heterogeneous and with complicated topology. In recent years, various ways have been developed to model these multiple-material domains [30,36,10]. However in comprehensive modeling and prediction of material behavior, there still lacks automatic and robust geometric modeling and mesh generation techniques, which are able to generate quality meshes for all material regions with correct topology. Resolving topology ambiguity is a challenging problem especially for multiple-material domains. In previous works, we have developed an octree-based Dual Contouring (DC) method to construct triangular and tetrahedral meshes for homogeneous domains with correct topology [33], by categorizing all possible ambiguity cases into 31 groups and modifying the dual mesh accordingly based on enumeration. However, this method only works for domains with one single material, and cannot handle complicated multiple-material domains like the polycrystalline material in Fig. 1. To address this deficiency, in this study we develop a novel approach to automatically and efficiently handle topology ambiguity for multiple material domains. ⇑ Corresponding author. Tel.: +1 412 268 5332; fax: +1 412 268 3348. E-mail address: [email protected] (Y. Zhang). 0045-7825/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2012.07.022
To distinguish topology ambiguities, we need to know the function property inside each cubic cell. Unlike single-material domains, where we can use a tri-linear interpolation function to model the interior domain, there is no such modeling function for multiple material cells. Instead, here we choose an indicator variable to examine each composing material and see whether topology ambiguity arises. Once one ambiguity is found, we split the cubic cell into 12 tetrahedral cells. The tessellation by tetrahedra guarantees to eliminate ambiguities in the resulting mesh. Till now we are able to modify the standard DC method based on such a hybrid octree consisting of both cubic and tetrahedral leaf cells. For each leaf cell, we insert one minimizer by either minimizing a pre-defined error function or simply choosing the cell center. Then, we classify all edges into cell edges, internal edges and face diagonals. In the following, we analyze each material-change edge, whose two end points lie in two different materials, and use its surrounding minimizers to form a polygon in the triangular mesh, or a polyhedron in the tetrahedral mesh. In the case of tetrahedralization, we also need to study each interior edge, whose two end points lie in the same material, to construct interior tetrahedra. Till now a mesh with no topology ambiguity is formed. However, we still need to improve the mesh quality via a combination of vertex classification, face swapping, edge removal and geometric flowbased smoothing. Our algorithm generates a conformal mesh and attains good quality with all topology ambiguity resolved. We have applied our algorithm to several complicated material datasets and achieved good results.
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
167
Fig. 1. (a) Adaptive triangular and tetrahedral meshes of a 92-grain polycrystalline material (multi-material domain). (b) Adaptive mesh between two interior grains; (c) Zoom-in picture of the planar surface. (d, e) zoom-in details of an interior local region with topology ambiguity (the red circle). The green line in (d) denotes one non-manifold edge, which is resolved in (e). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
The remainder of this paper is organized as follows: Section 2 reviews related previous work. Section 3 talks about the DC methods. Sections 4 and 5 discuss resolving 2D and 3D ambiguities. Section 6 presents how to improve the mesh quality. Section 7 illustrates results, and finally Section 8 draws conclusions. 2. Previous work 2.1. Isocontouring Isocontouring is the process of generating a piece-wise linear approximation to a given zero-surface of an implicit function. Of the various contouring methods, the Marching Cubes (MC) and DC are of most importance. The MC technique [15] is able to generate a closed, manifold triangular mesh for any scalar domain. To capture sharp features, additional information such as surface normals are used to handle sampling alias [8]. In addition, various techniques were developed to triangulate over cells with different levels. However, there are cracks when the resolution levels of adjacent cells are different, and a fan of triangles needs to be inserted around the gravity center of the coarse triangles [29]. Moreover, people developed ways to generate tetrahedral meshes involved with the body centered cubic lattice [9]. Topology ambiguity issues were awared of and studied [1,18,14]. Differently, the DC technique [7,24] uses an adaptive octree and produces water-tight contours over adaptive octree grids. However, it fails to recognize the ambiguous situations in the octree cells.
be in two separate pieces or there may be a single ‘‘tunnel’’ piece [18,1]. These ambiguities can be tracked using the collapsing rule [5], which repeatedly collapses the edges of the cell whose end points lie on the same side of the iso-surface. If the collapsing result is not an edge or point, the cell is ambiguous. To study possible topologies, the interior part of the cube is usually modeled as trilinear interpolations. After discussing the property of such tri-linear interpolation function, we obtain a classification of 31 cases [1]. With such classification, multiple attempts to solve the ambiguity problem were made [18,2,14]. These techniques are all based on the MC algorithm, and they model the iso-surface cell by cell. These methods provide a comprehensive solution to the ambiguity problem, however, they only work for single-material domains. 2.3. Image-based modeling for microstructures Image-based modeling is an important tool to study the properties and to predict the response of microstructures. With the development of automated collection and advanced image processing techniques, characterization of microstructure materials becomes tractable [6,22,26]. The collected datasets often provide statistically important representations of the microstructure. Moreover, they can be incorporated into the modeling and simulation for the design and prediction of material behaviors, which has been a cutting-edge research area in recent years [12,17,16]. People have developed various ways to model the microstructures [30,36,10]. However, an efficient and robust meshing technique that can handle topology ambiguity is still lacking.
2.2. Resolving topology ambiguity 3. A review of dual contouring methods Topology ambiguities arise when there are multiple choices of topology configurations for the given iso-surface [28]. It includes face ambiguity, when the four edges of a cube face intersect the iso-surface [20]; and body ambiguity, where the iso-surface may
Since our algorithm is based on the DC method, let us first review it. DC [7] is a grid-based method for contouring a scalar field. With a uniform grid, the memory consumption on complex models
168
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
with small features can be very high. Therefore we choose an adaptive octree, which has denser cells near the boundary or feature, and coarser cells inside the volume. There are three steps in our algorithm: construct an adaptive octree, compute a Quadratic Error Function (QEF) [3] for each leaf, and then generate dual meshes for this octree. The adaptive octree is achieved by subdividing the volume into octree cells recursively until the desired surface accuracy is achieved. The subdivision procedure begins with one octree cell that contains the whole domain, and it ends when the octree has achieved a certain adaptation level. The octree adaptation can be controlled with a feature sensitive function [34], which measures the iso-contour difference between two neighboring octree levels. The feature sensitive function is defined as P F ¼ ðjf iþ1 f i j=jrf i jÞ, where f i is a tri-linear interpolation function inside the cell at Level i [34]. The subdivision is carried out until the error is less than a pre-defined threshold. In order to obtain triangles with good quality, we restrict the neighboring octree level difference to be less than or equal to two. The next step is to construct a predefined QEF and to obtain a minimizer for each leaf cell by minimizing this function. The QEF is defined P as QEFðxÞ ¼ ðni ðx pi ÞÞ2 , where pi and ni are the position and unit normal vectors at the intersection points between the isosurface and octree cell edges [3,4]. The final step is to analyze each sign-change edge in manifold domains, which is defined as one edge whose two end points lie on different sides of the isocontour. Since the adaptive octree consists of leaves with different resolutions, each leaf may have neighbors at different levels. In order to avoid generating redundant elements, we always choose the minimal edge that does not contain an edge of the neighboring leaf. In this adaptive octree, the minimal edge is shared by either four or three cells, and we connect the minimizers of these cells to obtain a hybrid mesh consisting of quadrilateral (quad) and triangular elements. Finally these quads are divided to form an adaptive triangular mesh. This DC method has been extended to tetrahedral [34] and hexahedral [32] mesh generation. In tetrahedral mesh, interior edges, whose two end points lie inside the volume, need to be analyzed in addition to sign-change edges. For each sign-change edge, its three or four surrounding minimizers and its interior end point construct a tetrahedron or a pyramid. For each interior edge, its three or four surrounding minimizers and its two interior end points construct a diamond or pyramid. These pyramids and diamonds are then split into tetrahedra. In hexahedral meshes, we analyze each interior grid point instead. In a uniform case, each grid point is shared by eight octree cells and we construct a hexahedron from the minimizers of these cells. Several variants of the DC method have been developed as well, which introduce better topology and mesh quality [19,23,24,31]. Parallel algorithms were also introduced so that DC can be performed directly on modern Graphics Processing Unit architectures [25]. Furthermore, the DC method has also been extended to meshing a domain with multiple materials [36]. Based on the DC methodology, the material-change edge is used to identify the interface between material regions. In the meshing process, both materialchange edges and interior edges are analyzed to construct tetrahedral meshes, and interior grid points are analyzed for hexahedral mesh construction. Despite being adaptive and feature-preserving, one drawback of the above DC methods is that it fails to handle ambiguity in the topology. We have developed a technique that is able to distinguish and resolve all topology ambiguities in homogeneous materials with enumeration [33]. However, resolving topology ambiguity in multiple-material domains still remains an open problem. We will address this issue in the following sections.
4. Resolving topology ambiguity in 2D Ambiguity arises when the topology in each cell is not uniquely defined. In homogeneous domains, the input data contains a scalar field. If the function value at a grid point is greater than a given isovalue, then it has a positive sign; otherwise it is assigned negative. In multiple-material domains, the input data contains the material ID for each grid point. These values form the sign/material configuration of each cell. Topology ambiguity is one situation where the same sign/material configuration may lead to different topologies in representing the iso-surface. One way for checking is to use the collapsing rule [5], which repeatedly collapses the edges of one quad/ cubic cell, whose corners have the same sign/material ID, to a single vertex. The quad contains no ambiguity if and only if the result of this reduction is a single edge or a single vertex. To resolve the ambiguity, we need to distinguish the ambiguous situations and then make a conformal mesh accordingly. In this section, we will discuss 2D ambiguity problems, and then talk about complex 3D problems in the next section. 4.1. Identifying 2D ambiguity To begin with, we now study the 2D ambiguities. Since the material IDs are only given at the four corners of each quadtree cell, here we choose an indicator variable to decide the material ID of interior nodes. In one quad cell, suppose there are n different materials, and their corresponding material IDs are ID1 ; ID2 ; . . . ; IDn , where nranges from 1 to 4. For each of the material i (i = 1–n) and grid point j (j = 0–3) of the quad cell, we assign the following indicator variable:
vi;j ¼
1 If IDi ¼ material ID of Grid Point j 0
Otherwise
To illustrate, we use the example in Fig. 2. In this quad cell, there are two materials, red (ID1 ) and green (ID2 ), so n = 2. For material ID1 , grid points P 1 and P 3 belong to it; therefore their indicator variables for ID1 are 1. On the other hand, the indicator variables of grid points P0 and P2 are 0, as shown in (b). Similarly, for material ID2 , grid points P0 and P2 belong to it; therefore their indicator variables for ID2 are 1. On the other hand, the indicator variables of grid points P1 and P3 are 0, as shown in (c). For each material i (i = 1 to n), we choose a bi-linear interpolation, Fðn; gÞ, to represent the indicator variable in the interior domain. The bi-linear interpolation can be written in terms of the function values at four corners F ij ði; j ¼ 0; 1Þ:
Fðn; gÞ ¼ F 00 ð1 nÞð1 gÞ þ F 01 ð1 nÞg þ F 10 nð1 gÞ þ F 11 ng:
ð1Þ
Given an interior point P, suppose there are n materials in the cell, and the material ID is simply the material ID whose indicator variable is the maximum. Depending on the material ID of the quad center, we summarize all the 2D configurations into 8 unique cases with the consideration of symmetry and complementary. In order to differentiate distinct situations, we use the labelling notations as follows: the first number denotes the number of materials in the cube; the second label distinguishes the material ID of the four nodes; and the third label denotes the resolving of topology ambiguities. The quad cell may only contain one material (Case 1); two materials (Cases 2-1, 2-2-1, 2-2-2); or three materials (Cases 3-1, 3-2-1, and 3-22); otherwise the cell has four distinct materials (Case 4). From Fig. 3, we can tell that ambiguity arises when one of its materials has topology ambiguity, i.e, when one of its materials does not satisfy the collapsing rule. To be specific, there are ambiguities in Cases 2-2-1, 2-2-2, 3-2-1, 3-2-2, where two nodes of the
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
169
Fig. 2. An example of indicator variables: (a) contains two materials and their materials IDs are ID1 (red), ID2 (green); (b-c) show the indicator variables of material ID1 and ID2 , respectively, where a solid circle denotes vi;j ¼ 1, and an empty circle denotes vi;j ¼ 0. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
same material occupy the cell diagonal. After resolving these ambiguities, there are a total of 8 cases for a 2D cell (see Fig. 3). Simple configurations: Cases 1, 2-1, 3-1, and 4 have no ambiguity. These cases can be meshed using the standard DC. Configuration 2-2: There are two possible cases: Case 2-2-1: Green grid points are joined inside the cell; and the green material spreads across the cell; Case 2-2-2: Green grid points are separated; and the green isosurface breaks into two pieces inside this cell. On the other hand, red grid points are joined. Configuration 3-2: There are two possible cases: Case 3-2-1: Red grid points are joined inside the cell; and the red material spreads across the cell; Case 3-2-2: Red grid points are separated; and the green isosurface breaks into two pieces inside this cell. 4.2. Handling 2D ambiguity in triangulation After distinguishing topology ambiguities, now we mesh the 2D multiple-material domain. For those simple cases, the standard DC is sufficient to represent the correct topology. However, in order to handle the ambiguous cases, modifications to the standard DC is necessary. As discussed above, ambiguities are due to the alternative choice of topology for a quad cell containing two opposite corners with the same material ID. In fact, these ambiguities can be removed entirely by tessellating the ambiguous quad cell with triangles [28].
The disadvantage of this scheme is that it creates more triangles in the final mesh. However, considering the small portion of ambiguous cells in the quadtree, the element number increase is negligible. Following the above principle, for those ambiguous cells we choose a quad diagonal along which the material ID does not change, and use this diagonal to split the quad into two triangles. If no such diagonal exists (Case 3-2-2), we choose either one of the two diagonals. These triangles are treated as ‘‘triangular cells’’, and therefore the original quadtree has become a hybrid tree consisting of quad and triangular cells. Till now the tessellation is done. Now, we construct meshes with respect to the hybrid tree. Material-change edges and interior edges are analyzed to construct triangular meshes for each material region. Material-change edge: Each material-change edge is shared by two boundary cells in a uniform case, and one minimizer is computed for each of them. Note that these cells may be quads or triangles. These two minimizers and each end point of the material-change edge construct a triangle. If one material is the background, then no triangle is created from it. Therefore, a total of one or two triangles are obtained. In adaptive case, we always consider the minimal edges. With such restriction, each material-change edge is also shared by two boundary cells, but these two cells may have different sizes. From these two cells we can also get two minimizers to construct two triangles.
Fig. 3. 2D iso-surface topology possibilities. The cell has up to four materials. The iso-surface (red, green, blue and yellow) intersects the cell. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
170
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
Fig. 4. The meshing of 2D topology ambiguities. The dash lines form quad/triangular cells, the grey nodes are minimizers, while the colored nodes are grid points. In (a, b), the red nodes are connected in the middle while in (c, d) the red nodes are separated. In (a) and (c) the red, blue and green color represents the material ID of each triangle; in (b) and (d) the color denotes the type of edges around which triangles are built. The yellow triangles are generated by analyzing material-change edges; and the tan triangles are generated by interior edges (cell diagonal). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 5. An adaptive local region with ambiguity. The rendering color represents different materials. (a) The red material spreads the central domain; and (b) the red domain is separated in the central domain. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Interior edge: Each interior edge is shared by two quad or triangle cells. For a boundary cell, we determine the minimizer using the QEF. Otherwise the center of the cell is chosen as the minimizer. This interior edge and each of the minimizer point construct a triangle, therefore two triangles are obtained. Uniform and adaptive grids only differ in cell size, and the same meshing method can be applied. We demonstrate the resolving and meshing of topology ambiguities for Case 3-2-1 (a, b) and Case 3-2-2 (c, d) in Fig. 4. Note in (a) and (c) the rendering color represents the material ID while in (b) and (d) the color denotes the edge type. Fig. 5 shows the meshing example of an adaptive ambiguous region. To make the drawing clear, we do not draw all the triangle edges. 5. Resolving topology ambiguity in 3D Resolving 3D topology ambiguity is much more complicated. Any grid point can belong to eight possible materials, yielding a total of 88 cases. Even after considering symmetry and complementary, there are still more than 300 unique cases, making enumeration impracticable. Therefore, a general approach is needed to solve this problem automatically without considering too much tedious details. Here we first introduce the categorization of topology ambiguities, then present how to resolve them, and finally discuss how to integrate these topologies into mesh generation. Topology ambiguity is caused by tessellating the space using cubes. There can be multiple choices of topologies for a cube with any face or body diagonal containing the same material ID. Again we can use the same collapsing rule [5] to find out the topology ambiguities. We can repeatedly collapse edges of the cubic cell whose two end points have the same material ID to a single vertex. If the result of this reduction turns out to be an edge or a vertex, the cubic cell is free of ambiguity.
Topology ambiguity can be categorized into two groups: face ambiguity and interior ambiguity. Face ambiguity arises when two grid points with the same material ID occupy one cubic face diagonal, as shown in Fig. 6(a) and (b). The red grid points belong to Material 1, while the green ones belong to Material 2. The four edges of this face are all material-change edges, and we need to decide how to triangulate this local region. Meanwhile, there are additional interior ambiguities inside the cubic cell, as shown in Fig. 6(c) and (d). In these situations the iso-surface of Material 1 can either be two separated surfaces or a tunnel piece passing through the cubic cell. These interior structures decide what the topology looks like inside this cell. 5.1. Identifying 3D ambiguity In 3D cells, the material IDs are only given at the eight corner points. In order to resolve ambiguity, we need to know the material ID of interior points. Once again we use an indicator variable. In one cubic cell, suppose there are n different materials, and their corresponding material IDs are ID1 ; ID2 ; . . . ; IDn , where n ranges from 1 to 8. For each of the material i (i = 1 to n) and grid point j (j = 0 to 7) of the cubic cell, we assign the following indicator variable:
vi;j ¼
1 If IDi ¼ material ID of Grid Point j 0
Otherwise
To illustrate, we use the example in Fig. 7. In this cubic cell, there are three materials, red (ID1 ), green (ID2 ) and blue (ID3 ), so n = 3. For material ID1 , grid points P 0 ; P3 and P5 belong to it; therefore their indicator variables for ID1 are 1. On the other hand, the indicator variables of all the other grid points are 0, as shown in (b). Similarly, for material ID2 , grid points P2 and P 6 belong to it; therefore their indicator variables for ID2 are 1. On the other hand, the indicator variables of all the other grid points are 0, as shown in (c). Finally, for material ID3 , grid points P1 P4 and P7 belong to it; therefore
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
171
Fig. 6. Examples of 3D face and interior ambiguities. The cubic cell has two grid points (red) belonging to Material 1 and all the other grid points (green) belonging to Material 2. The iso-surface (red) intersects with this cell. (a, b) One face ambiguity; and (c, d) one interior ambiguity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 7. An example of indicator variables: (a) contains three materials and their material IDs are ID1 (red), ID2 (green) and ID3 (blue); (b–d) show the indicator variables of material material ID1 ; ID2 and ID3 , respectively, where a solid circle denotes vi;j ¼ 1, and an empty circle denotes vi;j ¼ 0. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
their indicator variables are 1. On the other hand, the indicator variables of all the other grid points are 0, as shown in (d). The indicator variable is tri-linearly interpolated in the interior domain. The tri-linear interpolation, Fðn; g; fÞ, is written in terms of the function values at the eight grid points, F ijk ði; j; k ¼ 0; 1Þ:
Fðn; g; fÞ ¼ F 000 ð1 nÞð1 gÞð1 fÞ þ F 001 ð1 nÞð1 gÞf þ F 010 ð1 nÞgð1 fÞ þ F 011 ð1 nÞgf þ F 100 nð1 gÞð1 fÞ þ F 101 nð1 gÞf þ F 110 ngð1 fÞ þ F 111 ngf:
ð2Þ
Given an interior point P, suppose there are n materials in the cell, the material ID is simply the material ID whose indicator variable is the maximum. For each material the cell contains, we study the topology classification according to the criteria used in the singlematerial domain [33]. The only difference is that instead of using the function value of the face/body saddle points, we now study the material IDs at the face and body centers of each cell. If any of the materials in the cube has topology ambiguity, the cell is considered ambiguous. With the topology classification known, we are able to generate the dual mesh and incorporate the specific topology into mesh generation. 5.2. Handling 3D ambiguity in tetrahedralization It is obvious that resolving 3D ambiguity is much more complicated. Since enumeration is impossible, we develop a general approach similar to 2D domains. This approach involves with splitting ambiguous cells and mesh generation. In fact, topology ambiguities can be removed entirely by tessellating the space using tetrahedra rather than cubes [28,27]. There are 24 ¼ 16 cases and after considering symmetry there are only five unique ones, see Fig. 8. The main disadvantage of this tessellating scheme is that more elements are generated in the final mesh.
To solve this problem, we only split the ambiguous cubic cells into tetrahedra. Because the number of ambiguous cells is small compared with the total cell number, this increase can be neglected. Based on this idea, we propose the following four steps to handle topology ambiguities. These four steps include: (1) find the ambiguous octree cells and split each of them into 12 tetrahedra, forming a hybrid tree consisting of cubes and tetrahedra; (2) calculate one minimizer for each cube and choose the center of each tetrahedron as the minimizer; (3) analyze each material-change edge and interior edge to build a volumetric mesh. Details of each step are described as follows. 5.2.1. Splitting ambiguous cells To eliminate the ambiguities, we first split each ambiguous cell into tetrahedra. Moreover, each of the six faces is split into two triangles using the same rule as in 2D ambiguous faces. The cell center and each of these triangles construct a tetrahedron. This scheme is able to resolve all the topology ambiguities. Face ambiguities are resolved automatically because the 2D scheme applied on cubic faces distinguishes 2D ambiguities. The body ambiguities are also resolved because the material ID of the body center decides either a ‘‘tunnel’’ or two separated pieces exist. To illustrate, we use the cube in Fig. 9 to show the splitting schemes. For better visualization, only the front face is split. If the front face belongs to 2D Case 3-2-1 (see Fig. 3), then we choose the diagonal shown in (a). Otherwise we choose the diagonal shown in (b). 5.2.2. Computing one minimizer for each hybrid cell With all topology ambiguities identified, we are able to compute one minimizer for each of the hybrid cells. There are two types of octree cells: boundary cells and interior cells. Boundary cells include cubic cells and tetrahedral cells. For cubic and tetrahedral cells, we follow the standard DC method to compute one
172
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
Fig. 8. Five possible configurations of one tetrahedral cell.
minimizers to form a volumetric mesh. For each edge type, we have different treatments as follows.
Fig. 9. Cube splitting schemes when red nodes are joined (a) or separated (b) on the front face.
minimizer by minimizing the predefined QEF. For other interior cells we also choose the center as the minimizer. Our approach assumes that only one minimizer is generated in each cell, and since all topology ambiguities are resolved beforehand, our approach is free of any topological mis-matches. 5.2.3. Building tetrahedral mesh by analyzing edges To construct tetrahedral meshes from the hybrid octree, we need to analyze each material-change edge as well as interior edges. All the material-change edges and interior edges belong to the cubic cells or the tetrahedral cells, and they can be cell edges, internal edges and face diagonals, as shown in Fig. 10. Cell edges (Fig. 10(a)) are surrounded by a combination of cubic and tetrahedral cells. There can be 4 to 8 surrounding minimizers, corresponding to the situations where the edge is shared by 4 cubic cells, and where the edge is shared by 8 tetrahedral cells, respectively. The internal edges (Fig. 10(b)) are those which connect the center to cell corners. These edges are surrounded by 3 to 6 tetrahedra. Face diagonals (Fig. 10(c)) can be shared by either 1 cubic cell and 2 tetrahedral cells, or 4 tetrahedral cells. Therefore there can be 3 or 4 surrounding minimizers. For adaptive octree, we follow the standard DC method and always use the minimal edges. We analyze each material-change edge and interior edge in the analysis domain, connect the QEF
Material-change edge: Each material-change edge is shared by several cubic or tetrahedral cells, and one minimizer is computed for each cell. We use these minimizers to form a convex polygon. Depending on the position of the edge end points, we now have the following two tetrahedralization schemes. (1) Material-change edge with one end point in the background material: For these edges, we first triangulate the convex polygon, and construct one tetrahedron by connecting each triangle with the edge end point not in the background, as shown in Fig. 11(a). The nodes V 1 V 3 V 2 V 0 form one tetrahedron, and V 1 V 4 V 3 V 0 form the other. (2) Other material-change edges: For these edges, we construct a polyhedron using the polygon and two edge end points, and then split it into tetrahedra. Note that the tetrahedra in both sides of the polygon belong to different materials. If we choose the tetrahedralization scheme in Fig. 11(a), we may experience inverted elements in adaptive octree. Therefore we split the polyhedron by inserting one edge center node, as shown in Fig. 11(b). We construct a tetrahedron using this center node and each triangle in the polyhedral surface. In Fig. 11(b), the node V 6 is the edge center node, and we connect it with every polyhedral surface triangle, for instance, triangle V 1 V 2 V 5 to construct a tetrahedron. Interior edges: For these edges, we construct a polyhedron using the polygon and two edge end points, and then split it into tetrahedra. Again in order to avoid inverted elements, we construct the tetrahedra using the two edge end points and one edge of the polygon. In Fig. 11(c), four tetrahedra are created, and they are V 5 V 0 V 2 V 1 , V 5 V 0 V 3 V 2 , V 5 V 0 V 4 V 3 and V 5 V 0 V 1 V 4 . Fig. 12 shows the resulting meshes of five local regions before and after resolving topology ambiguity. To show the connecting
Fig. 10. Three kinds of edges (red): cell edges (a), internal edges (b) and face diagonals (c). Only the adjacent tetrahedra are drawn. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
173
Fig. 11. The tetrahedralization of three types of edges: material-change edges with one end point lying in background (a), other material-change edges (b), and interior edges (c). Different meshing schemes are shown, and the material ID of each tetrahedron is rendered in color. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 12. Meshes of five regions before (first row) and after (second row) the re-connection. The red circles indicate the locations of the ambiguity, green lines or the blue dot indicates non-manifold situations, and red dashed lines denote the axes of a tunnel topology. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
schemes more clearly, we only show the mesh of one material. The green lines or the blue node denotes non-manifold situations, where one edge is shared by more than two triangles on the surface, or one vertex whose neighborhood does not form a disc. Note in (c), along the red dashed axes, a hollow tunnel topology is missing, and our re-connection algorithm recovers it. In these five figures, it is obvious that the topology ambiguities are resolved and appropriate meshes are built. 6. Quality improvement Till now the generated mesh has correct topology and the element quality is good especially at the core; however, it may have poor-quality elements around the material boundary that influence the convergence and stability of the numerical simulation. Therefore, we need to take an additional step to improve the mesh quality. Our quality improvement scheme includes vertex classification, face swapping, edge removal, and geometric flow-based smoothing. To begin with, we need metrics to measure the quality of the triangular and tetrahedral meshes. For triangular meshes, the metric we choose is the aspect ratio, defined as 2 r in =r out , where rin is the radius of the inscribed circle, and rout is the radius of the circumcircle. For tetrahedral meshes, we choose the edge-ratio, a minimum volume bound and the Joe-Liu parameter [13] to measure the mesh quality. The edge-ratio is the ratio of the longest
edge length over the shortest edge length of one tetrahedron; the minimum volume bound is the minimum volume of all the tetrahedra in the mesh; the Joe-Liu parameter is defined as follows:
24=3 3 jsj2=3 Fðs; dÞ ¼ P ; 2 06i
ð3Þ
where jsj denotes the volume of one tetrahedron, V i (i = 0 to 3) are the four vertices, and eij denotes the edge connecting V i to V j . With these measures, the mesh quality can be judged by observing the worst element quality, and the distribution of elements in terms of their quality values. For microstructures, the datasets are defined in a cubic representative volume element (RVE). We need to preserve the RVE boundary during quality improvement. According to the relative positions to the RVE and to each grain, all the vertices can be categorized into seven distinct groups, as shown in Fig. 13 and each group needs different improving techniques. Group 1: The eight corners of the RVE box, which are fixed during the improvement in order to keep the RVE boundary. Group 2: Vertices on the twelve edges of the RVE box, which only move along the edge. If they lie inside one grain (2a), they are smoothed only along the edge; if they are shared by two or more than two grains (2b), they are fixed during the improvement. Group 3: Vertices on the six faces of the RVE box, which are smoothed only on the plane. If they are shared by more than two grains (3a), they are planar non-manifold points and are fixed dur-
174
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
Fig. 13. Vertex classification. (a) Groups 1–3 on the RVE boundary; and (b) Groups 4–7 inside the RVE volume.
ing the improvement; If they are on the planar curve shared by two grains (3b), they are smoothed along the tangent direction of the planar curve during the improvement; finally if they are the vertices of one grain (3c), they are only moved on the plane. Group 4: Vertices inside one grain, which will be improved using the weighted-averaging method. These vertices do not exist in triangular meshes. Group 5: Vertices located on the grain boundary surface patches shared by two grains, which are smoothed on the tangent plane of the grain boundary during the improvement. Group 6: Vertices located on interior curves, which can only move along the tangent direction of the interior curve. Group 7: Interior non-manifold vertices shared by more than two interior curves, which are fixed during the improvement. During the improvement process, the vertex classification and the corresponding moving method need to be followed. For instance, when we relocate a vertex in Group 3, the movement should be projected to the RVE plane. For planar or interior curves, we use a curve fairing method based on the mean curvature flow to improve them. We fix the starting and ending non-manifold points (Groups 3a, 7), and move the interior vertices (Groups 3b, 6) along the tangent direction using a B-spline interpolation and re-sampling technique [21]. For vertices on grain boundary surface patches (Group 5), we improve them using a surface diffusion flow developed in [35]. Finally for the interior vertices (Group 4), we relocate them to the mass center of their neighboring elements. In many cases, elements with bad aspect ratios may still remain in the mesh after smoothing. Therefore, we consider modifying the mesh topology and use it together with smoothing to improve the quality. In adaptive meshes, there might be vertices embedded in a triangle with an element valence number of 3 (triangular mesh), or
a tetrahedron with an element valence of 4 (tetrahedral mesh). Such vertices are first removed. In addition, if one tetrahedron has two edges lying on the same boundary curve, these two edges form a big angle and the local quality is usually bad. To handle them, we change the local connectivities to make sure that each tetrahedron has at most one edge on the boundary curve. Furthermore, face swapping and edge removal are also used to improve the local mesh topology [11]. The face swapping operation swaps the faces of several neighboring tetrahedra. There are four basic templates including 2–3 flip, 3–2 flip, 2–2 flip and 4–4 flip, as shown in Fig. 14. One swapping operation is legal only when it does not generate inverted elements. The algorithm begins by checking the Joe-Liu parameter of all the tetrahedra. When the quality is less than a given threshold, we perform the face swapping technique. For each tetrahedron, we need to check all its 4 faces and 6 edges. The following procedures are taken: (1) Check 4 faces: We examine each face of the tetrahedron. Whenever a 2–3 flip is legal, we save its resulting local qualities. (2) Check 6 edges: We analyze each edge of the tetrahedron. Depending on the element valence number of that edge, we can use 3–2 flip (the valence number equals to 3), 2–2 flip (the valence number equals to 2 and the edge lies on the RVE boundary) or 4–4 flip (the valence number equals to 4). Once a flip is legal, we also save its resulting local qualities. (3) Carry out the flip: After checking all its faces and edges, we pick the one that leads to the best local quality and carry out the face swapping operation according to the templates.
Fig. 14. Face swapping operations: (a, b) 2–3 and 3–2 flip; (b) 2–2 flip; and (c) 4–4 flip.
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
175
Fig. 15. An edge removal operation. (a) The original polyhedron with the red center edge; (b) the center edge and the surrounding minimizers; and (c) the resulting mesh after edge removal. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 16. Zoom-in details of a 308-grain polycrystalline material mesh (multi-material domain) before (a) and after (b) the quality improvement.
Edge swapping is another efficient approach in eliminating sliver tetrahedral elements. We use an edge removal operation to swap interior edges. These edges are shared by multiple tetrahedra. In Fig. 15, we show the removal of one red edge. In this process the middle polygon needs to be triangulated in the most optimized way. Note that to preserve boundary features of the multiple-material domain, the face swapping and edge removal operations cannot be carried out across grain boundaries. Fig. 16 shows the effects of our quality improvement. 7. Results and discussion Our developed algorithm is automatic and robust in resolving topology ambiguity for both homogeneous and multiple-material domains. We have applied our algorithm to two multiple-material datasets (Figs. 1 and 18) as well as one homogeneous dataset (Fig. 17). The computations were carried out on a PC with Intel Q9600 CPU and 4 GB DDR-II memories. Tables 1 and 2 show statistics of the generated triangular and tetrahedral meshes. Fig. 19 shows the histogram of aspect ratios for triangular meshes and Joe-Liu parameters for tetrahedral meshes. The minimum– maximum octree level of the Foam, 92-grain and 308-grain models are 5–8, 4–7 and 4–7, respectively. The mesh size is increased after resolving the ambiguities, since one cubic leaf cell may be split into 12 tetrahedral cells. The mesh quality is improved greatly after quality improvement. Note that the computation time includes the time for the mesh generation, resolving topology ambiguities, together with quality improvement, which takes most of the time. Our technique combines the idea of DC and Marching Tetrahedra. The spatial grid is a hybrid tree structure consisting of cubic and tet-
rahedral cells. Since the percentage of ambiguous cells is small (usually 1 in 1000 cells), such combination solves the ambiguity problem without increasing the number of vertices too much. Furthermore, compared with the topology preservation technique developed for homogeneous materials [36], which categorizes all topologies into 31 cases, our new technique has two advantages. First, it is able to handle both homogeneous and multiple-material domains. In this algorithm, the homogeneous material is treated as one special case of multiple-material domains. Second, instead of based on enumeration, our new algorithm is more general. This makes it much easier to implement and extend to tetrahedral mesh generation. For our octree-based method, it comes with several inherent limitations. Bad quality elements are generated along the boundary although good quality elements are present at the core, and the obtained mesh quality may be inferior to Delaunay meshes. For complex models with small local features, the memory consumption can be very high. Therefore it is critical to choose an adaptive octree and generate adaptive meshes. During the octree subdivision, the maximum depth is determined by surface accuracy. The maximum number of materials is limited to four for 2D domains and eight for 3D domains. This can be overcome by refining the cells locally. 8. Conclusion In this study, we have developed an efficient and automatic algorithm to resolve topology ambiguity in generating triangular and tetrahedral meshes for multiple-material domains. We begin with an indicator variable to decide the topology inside each octree leaf cell. When one ambiguity is detected, we split the cell into 12
176
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
Fig. 17. (a) Adaptive triangular and tetrahedral meshes of a foam (homogeneous) material. (b) shows one local region of the tetrahedral mesh with some elements removed; (c) shows a zoom-in picture of one local boundary; (d) and (e) show zoom-in details of an interior local region with topology ambiguity (the red circle). The green line in (d) denotes one non-manifold edge, which is resolved in (e). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 18. (a) Adaptive triangular and tetrahedral meshes of a 308-grain polycrystalline material (multi-material domains). (b) shows the adaptive mesh between two interior grains; (c) shows a zoom-in picture of the planar surface. (d) and (e) show zoom-in details of an interior local region with topology ambiguity (the red circle). The green line in (d) denotes one non-manifold edge, which is resolved in (e). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
177
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178 Table 1 Statistics of the generated triangle meshes. Model
Mesh
Mesh size (Vert #, Elem #)
Aspect ratio (Worst, Best)
Ambiguity Cell #
Time (s)
Foam (Fig. 17)
Original Improved
(188,352, 378,755) (189,423, 380,412)
(0.01, 1.0) (0.10, 1.0)
214 214
135 161
92 Grain data (Fig. 1)
Original Improved
(15,789, 34,918) (16,000, 36,286)
(0.01, 1.0) (0.13, 1.0)
30 30
14 21
308 Grain data (Fig. 18)
Original Improved
(17,092, 40,715) (18,713, 42,092)
(0.04, 1.0) (0.14, 1.0)
14 14
23 30
Table 2 Statistics of the generated tetrahedral meshes. Model
Mesh
Mesh size (Vert #, Elem #)
Edge ratio (Worst, Best)
Volume (Min, Max)
Joe–Liu parameter (Worst, Best)
Time (s)
Foam (Fig. 17)
Original Improved
(553,612, 2,712,452) (554,706, 2,721,620)
(411.5, 1.03) (11.50, 1.01)
(1.0e5, 0.87) (6.2e2, 0.67)
(2.0e6, 1.0) (0.1, 1.0)
258 395
92 Grain data (Fig. 1)
Original Improved
(76,765, 407,525) (77,333, 408,664)
(312.32, 1.01) (13.21, 1.0)
(1.2e4, 1.52) (2.7e-3, 1.20)
(2.5e4, 0.98) (0.10, 1.0)
37 101
308 Grain data (Fig. 18)
Original Improved
(92,110, 494,455) (92,769, 495,353)
(631.44, 1.01) (13.12, 1.02)
(1.0e6, 1.30) (4.6e3, 0.84)
(1.0e6, 1.0) (0.10, 1.0)
65 123
Fig. 19. The histogram of aspect ratios for triangular meshes (a) and Joe-Liu Parameters for tetrahedral meshes (b).
tetrahedra. In the following we study each material-change edge, and use its surrounding minimizers to create triangles or boundary tetrahedra. In addition, we analyze each interior edge to create interior tetrahedra. Finally, we apply quality improvement techniques to attain a quality mesh. In this work we build up a hybrid octree and generate quality meshes with all topology ambiguities resolved. We have applied our algorithm to three complicated datasets and obtained good results. In the future we would like to test more multiple material datasets in various application fields. Acknowledgement The research was supported in part by an AFOSR grant FA955011-1-0346 and a NSF/DoD-MRSEC seed grant. References [1] E. Chernyaev, Marching cubes 33: construction of topologically correct isosurfaces, Technical Report CH/95-17, CERN, 1995. [2] P. Cignoni, F. Ganovelli, C. Montani, R. Scopigno, Reconstruction of topologically correct and adaptive trilinear isosurfaces, Computers and Graphics 24 (2000) 399–418.
[3] M. Garland, P. Heckbert, Simplifying surfaces with color and texture using quadric error metrics, in: IEEE Visualization’98, 1998, pp. 263–270. [4] M. Garland, E. Shaffer, A multiphase approach to efficient surface simplification, in: IEEE Visualization’02, 2002, pp. 117–124. [5] T. Gerstner, R. Paarola, Topology preserving and controlled topology simplifying multiresolution isosurface extraction, in: IEEE Visualization’00, 2000, pp. 259–266. [6] M.A. Groeber, B.K. Haley, M.D. Uchic, D.M. Dimiduk, S. Ghosh, 3D reconstruction and characterization of polycrystalline microstructures using a FIBCSEM system, Materials Characterization 57 (4–5) (2006) 399–418. [7] T. Ju, F. Losaasso, S. Schaefer, J. Warren, Dual contouring of Hermite data, ACM Transactions on Graphics 21 (2002) 339–346. [8] L. Kobbelt, M. Botsch, U. Schwanecke, Feature sensitive surface extraction from volume data, in: SIGGRAPH’01, 2001, pp. 57–66. [9] F. Labelle, J. Shewchuk, Isosurface stuffing: fast tetrahedral meshes with good dihedral angles, ACM Transactions on Graphics 26 (3) (2007) 57.1– 57.10. [10] S.A. Langer, E.R. Fuller Jr., W.C. Carter, OOF: an image-based finite-element analysis of material microstructures, Computing in Science Engineering 3 (3) (2001) 15–23. [11] J. Leng, Y. Zhang, G. Xu, A novel geometric flow-driven approach for quality improvement of segmented tetrahedral meshes, in: International Meshing Roundtable, 2011, pp. 347–364. [12] A.C. Lewis, A.B. Geltmacher, Image-based modeling of the response of experimental 3D microstructures to mechanical loading, Scripta Materialia 55 (1) (2006) 81–85. [13] A. Liu, B. Joe, Relationship between tetrahedron shape measures, BIT Numerical Mathematics 34 (2) (1994) 268–287.
178
Y. Zhang, J. Qian / Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 166–178
[14] A. Lopes, K. Brodlie, Improving the robustness and accuracy of the marching cubes algorithm for isosurfacing, IEEE Transactions on Visualization and Computer Graphics 9 (1) (2003) 16–29. [15] W. Lorensen, H. Cline, Marching cubes: a high resolution 3D surface construction algorithm, in: SIGGRAPH’87, vol. 21, 1987, pp. 163–169. [16] J. MacSleynea, M.D. Uchicb, J.P. Simmonsb, M. De Graefa, Three-dimensional analysis of secondary precipitates in René-88 dt and UMF-20 superalloys, Acta Materialia 57 (20) (2009) 6251–6267. [17] J. Madisona, J. Spowartb, D. Rowenhorstc, L.K. Aagesend, K. Thorntona, T.M. Pollocka, Modeling fluid flow in three-dimensional single crystal dendritic structures, Acta Materialia 58 (8) (2010) 2864–2875. [18] B. Natarajan, On generating topologically correct isosurfaces from uniform samples, The Visual Computer 11 (1994) 52–62. [19] G. Nielson, Dual marching cubes, in: IEEE Visualization’04, 2004, pp. 489– 496. [20] G. Nielson, B. Hamann, The asymptotic decider: resolving the ambiguity in marching cubes, in: IEEE Visualization’92, 1992, pp. 83–91. [21] J. Qian, Y. Zhang, W. Wang, A.C. Lewis, M.A.S. Qidwai, A.B. Geltmacher, Quality improvement of non-manifold hexahedral meshes for critical feature determination of microstructure materials, International Journal for Numerical Methods in Engineering 82 (11) (2010) 1406–1423. [22] D.J. Rowenhorst, A.C. Lewis, G. Spanos, Three-dimensional analysis of grain topology and interface curvature in a b-titanium alloy, Acta Materialia 58 (16) (2010) 5511–5519. [23] S. Schaefer, J. Warren, Dual marching cubes: primal contouring of dual grids, in: Computer Graphics and Applications’04, 2004, pp. 70–76. [24] S. Schaefer, T. Ju, J. Warren, Manifold dual contouring, IEEE Transactions on Visualization and Computer Graphics 13 (3) (2007) 610–619. [25] L. Schmitz, L.F. Scheidegger, D.K. Osmari, C.A. Dietrich, J.L.D. Comba, Efficient and quality contouring algorithms on the GPU, Computer Graphics Forum 29 (8) (2010) 2569–2578.
[26] Z. Shan, A.M. Gokhale, Micromechanics of complex three-dimensional microstructures, Acta Materialia 49 (11) (2001) 2001–2015. [27] B.S. Sohn, Topology preserving tetrahedral decomposition applied to trilinear interval volume tetrahedrization, KSII Transaction one Internet and Information Systems 3 (6) (2009) 667–681. [28] G.M. Treece, R.W. Prager, A.H. Gee, Regularised marching tetrahedra: improved iso-surface extraction, Computers and Graphics 23 (1998) 583–598. [29] J.R. Westermann, L. Kobbelt, T. Ertl, Real-time exploration of regular volume data by adaptive reconstruction of isosurfaces, The Visual Computer 15 (2) (1999) 100–111. [30] Z. Wu, J.M. Sullivan Jr., Multiple material marching cubes algorithm, International Journal for Numerical Methods in Engineering 58 (2) (2003) 189–207. [31] N. Zhang, W. Hong, A. Kaufman, Dual contouring with topology preserving simplification using enhanced cell representation, in: IEEE Visualization’04, 2004, pp. 505–512. [32] Y. Zhang, C. Bajaj, Adaptive and quality quadrilateral/hexahedral meshing from volumetric data, Computer Methods in Applied Mechanics and Engineering 195 (9–12) (2006) 942–960. [33] Y. Zhang, J. Qian, Dual contouring for domains with topology ambiguity, Computer Methods in Applied Mechanics and Engineering 217–220 (2012) 34–45. [34] Y. Zhang, C. Bajaj, B.-S. Sohn, 3D finite element meshing from imaging data, Computer Methods in Applied Mechanics and Engineering 194 (2005) 5083– 5106. [35] Y. Zhang, C. Bajaj, G. Xu, Surface smoothing and quality improvement of quadrilateral/hexahedral meshes with geometric flow, Communications in Numerical Methods in Engineering 25 (1) (2009) 1–18. [36] Y. Zhang, T. Hughes, C. Bajaj, An automatic 3D mesh generation method for domains with multiple materials, Computer Methods in Applied Mechanics and Engineering 199 (5–8) (2010) 405–415.