polyhedral mesh generation
POLYHEDRAL MESH GENERATION W. Oaks1, S. Paoletti2 adapco Ltd, 60 Broadhollow Road, Melville 11747 New York, USA 1
2
[email protected], [email protected]
ABSTRACT A new methodology to generate a hex-dominant mesh is presented. From a closed surface and an initial all-hex mesh that contains the closed surface, the proposed algorithm generates, by intersection, a new mostly-hex mesh that includes polyhedra located at the boundary of the geometrical domain. The polyhedra may be used as cells if the field simulation solver supports them or be decomposed into hexahedra and pyramids using a generalized mid-point-subdivision technique. This methodology is currently used to provide hex-dominant automatic mesh generation in the preprocessor pro*am of the CFD code STAR-CD. Keywords: mesh generation, polyhedra, hexahedra, mid-point-subdivision, computational geometry.
1. INTRODUCTION In the past decade automatic mesh generation has received significant attention from researchers both in academia and in the industry with the result of considerable advances in understanding the underlying topological and geometrical properties of meshes. Most of the practical methodologies that have been developed rely on using simplical cells such as triangles in 2D and tetrahedra in 3D, due to the availability of relatively robust algorithms for such shapes. Nevertheless, in the industry, users of finite volume and finite element technologies in 3D tend to prefer hexahedral cells because in many cases hexahedra can better capture the characteristics of the solution and because a larger number of tetrahedral (arguably 3 to 5 times) is needed to achieve an equivalent accuracy in the solution. As a result, truly automatic unstructured hexahedral mesh generation is one of the most desired, and not completely fulfilled, features in a mesh generator. Several methodologies have been proposed to automatically create an unstructured hexahedral mesh. These methodologies can be ideally divided into two families: the first one starts the generation from a given closed surface tessellated with quadrilaterals and attempts
to create a mesh that satisfies the boundary constraint, while the second family does not require a quadrilateral surface mesh. Its creation becomes part of the generation process. The first family of methodologies such as spatial twist continuum STC and related whisker weaving [1], hexahedral advancing-front [2], rely on a proof of existence [3],[4] for combinatorial hexahedral meshes whose only requirement is an even number of quadrilaterals on the boundary surface. Although the methods of the first family have provided an enormous improvement in understanding the topological properties of hexahedral meshes, their application to reallife geometries has been limited by the mainly topological nature of the algorithms and by problems in ’closing’ the mesh i.e. in generating the last cells. In fact for several simple closed configurations of quadrilaterals, which should admit a hexahedral mesh, there is no known valid mesh. The second family of methodologies includes those based on medial surface [5],[6] and the so-called grid-based methods [7],[8],[9],[10]. The medial surface technique creates a two-dimensional skeleton of a three-dimensional volume as the locus of all
inscribed spheres of maximal diameter. The medial surface is then used to decompose the volume into parts meshable by existing techniques.
accommodate the solution of the PDE's. In fact all current FV methods use polyhedra. It happens to be that these polyhedra are hexahedra, prisms, pyramids or tetrahedra.
The so-called grid methods start covering the domain with a Cartesian hexahedral grid and removing from it the cells that either are outside the domain or are closer than a specified distance to the boundary. The boundary of this ‘initial mesh’ is made by quadrilaterals that can be extruded up to the surface, generating an isomorphic hexahedral mesh on the boundary.
Here we propose the use of the polyhedra that result from the intersection of an initial all-hex mesh that cover the whole domain of interest with a triangulated boundary surface. The result of the intersection or 'trimming' is the set of hexahedra of the initial mesh which lie totally inside the boundary surface plus a set of polyhedra generated by the intersection of the 'surface-related' hexahedra with the surface itself.
One of the problems of the medial-surface-based algorithms is the abundance of degeneracies in the medial surface leading to considerable topological complexity. The grid-based methods are relatively simple and robust. One limitation lies in the representation of the discontinuities (edges and corners) of the boundary surface. In fact it would be difficult for a mesh generated with such method to match a corner in the surface with an arbitrary number of sharp edges attached to it. A carefully handmade initial mesh could overcome the limitation, but this would also undermine the automatic character of the mesh generation.
1.1 Motivations The all-hex methods just mentioned, although very promising, may present some limitations when applied to the most extreme geometrical cases that can be encountered in the industry-related Computational Fluid Dynamics (CFD) and Stress Analysis (SA) simulations. Hexdominant algorithms have been proposed in an attempt to overcome the problem. The hex-dominant algorithms use one (or more) all-hex method but, when the all-hex method fails or in regions that are likely too complex for such method, then create more tractable shapes (typically tetrahedra, prisms and pyramids). In this paper we present a different approach to generate a hex-dominant mesh. If the constraint of using only hexahedra is released, it becomes acceptable to generate a hex-dominant mesh using any kind of shape or more generally a polyhedron, provided that the overwhelming majority of the mesh is still composed of hexahedra.
1.2 Polyhedra in the framework of Finite Volume and Finite Element methods Meshing may be defined as the process of breaking up a geometrical domain into smaller and geometrically simpler sub-domains (the cells) which should allow the solution of the discretized partial differential equations (PDE) that describe the behaviors of the fields involved in the physics of the phenomenon. In CFD one of the most popular ways of solving the Euler or Navier-Stokes equations is the Finite Volume (FV) method, which is based on the concepts of cell faces and flux through those faces. The method itself does not depend on the shape of the cells if the cell faces can be correctly defined. In this case a (convex) polyhedron may be used to
The method has a certain resemblance with the grid-method proposed in [7]. In fact it would create the same internal hexahedral structure, while the hexahedral buffer, which in Schneiders’ method is created using an isomorphic extrusion of the boundary faces of the interior mesh, here is replaced by the polyhedra created by the intersection. Since polyhedra with an arbitrary number of faces and arbitrary node valence (here valence is defined as the number of edges sharing a node) may represent arbitrary sharp edges and corners on the surface, the problem of a correct representation of the surface features which may affects Schneiders’ method is eliminated. On the one hand the direct usage of polyhedra as cells requires a solver able to handle them so that STAR-CD list of acceptable cell shapes has been extended to include a certain number of simple polyhedra. On the other hand, if an extension of the mid-pointsubdivision [11] is applied to the whole new mesh (including the polyhedra), as explained later, the result is a set of hexahedra and a much smaller set of pyramids. As a result, the proposed methodology has the potential to generate meshes for all the solvers that accept hexahedra and pyramids, i.e. the majority of Finite Elements (FE) and FV programs.
2. THE INITIAL MESH The initial mesh is a hexahedral mesh that must cover the whole domain of interest. The initial mesh may be structured or unstructured. Trivial ways to create the initial mesh include: •
Cartesian mesh generation
•
Cylindrical mesh generation (if the domain to be meshed presents an axis of symmetry).
•
Trans Finite Interpolation [12].
Moreover any combination of such techniques may be used to mesh different parts of the domain. We found that unmodified initial grids may lead to the generation of very small cells when the surface is close to one or more nodes of the cells. To overcome this difficulty we perform an adjustment procedure on the nodes of the initial grid.
We propose two different simple methods of accomplishing the task: 1.
Projecting the nodes of the initial mesh onto the surface when they are close to the surface. Here close means that the distance is smaller than a specified small fraction of the cell size.
2.
Moving the nodes of the initial mesh away from the surface when they are initially close to the surface.
To move a node far from the surface, first we find the vector from the node to its projection onto the surface, and then we move the node in the opposite direction up to a specified distance from the surface Another reason to move the mesh nodes close to the surface onto the surface or away from the surface is that the adjustment modifies the topology of the final mesh. In a polyhedron a trivalent node is a node that has exactly three edges and three faces attached, while for a nontrivalent node the number of edges and faces is not equal to three (typically it is greater). If all the nodes of a polyhedron are trivalent the polyhedron is trivalent. An interesting property of trivalent polyhedra is that when a b-rep solid represented by a trivalent polyhedron is intersected by another b-rep solid, the result is another trivalent polyhedron if the faces of the two solid intersect only at the edges and not at the corners, their sharp edges have no common point and there is no corner with more than three sharp edges inside the trivalent polyhedron. The proposed methods generate quite different final polyhedral cells. Method 1 tends to generate a larger number of polyhedra with non-trivalent nodes than the second one. In fact, the initial mesh after the adjustment has several nodes on the surface, and the intersection of a cell with a node on the surface with the surface itself may create a non-trivalent polyhedron. On the other hand it has the advantage of generating simpler polyhedra with a smaller number of faces and this is a desirable property because not all the polyhedra can be accepted as valid cell by STAR-CD (see 3.1.) If method 2 is applied, no node of the initial mesh is on the surface. The intersection creates a non-trivalent polyhedron only if the cell edges and the sharp edges of the surface intersect or if a surface corner with more that three sharp edges attached is inside the cell. This characteristic is important to maximize the number of the resulting hexahedra, when using the generalized mid-pointsubdivision (MPS), as we will explain later. The initial mesh does not need to have the property of conformity (for a definition of conformity see[8]) and an initial fully connected mesh may be locally and recursively refined according to several criteria including: •
Distance from the boundary
•
Curvature of the closest boundary
•
Interactive user input
A cell can be marked by one of the above criteria and split into 8 cells. All the information related to the cell neighbors is automatically updated. Other interesting improvements are under evaluation. For examples it would be useful to automatically create a roughly body-fitted initial mesh.
2.1 Basic Methodology We may consider the core of the polyhedral mesh generation as an intersection operation between two b-rep solids: the triangulated boundary surface S and the initial hexahedral cell C, whose faces, without loss of generality, can be triangulated in a unique manner choosing a point on each face (the face centroid) and connecting it to the four edges, creating four triangles per face. Once the triangulation is performed we have two b-rep solids whose boundaries are triangular faces. There are fast ways to evaluate approximated intersections such as the Marching Cube Algorithm[13], but it tends not to resolve the surface discontinuity leading to a poor representation of corners and sharp edges. Since it may be important to represent such details in the final mesh, we chose to perform an exact intersection. Given two polyhedral solids S and C, the intersection is a Boolean operation and the result P is yet another polyhedral solid: P = S ∩C (1) that is formed by: •
All components of S that lie inside C and all components of C that lie inside S.
There are several methods available in literature to perform Boolean operations on b-rep solids (a small bibliography can be found in [14]). However, the computational cost of several thousands or even hundreds of thousands of intersections between the cells of the initial mesh and the surface can be extremely high. An efficient implementation of the intersection operation is important for the practical applicability of this mesh generation method. As an extension of the boundary operator defined in [15], let us define b() as the operator that extracts the boundary of an entity, providing as result a set of one-dimensionlower entities. For completeness we may also define the identity operator as I(X)=X. Applying the b() to eqn.(1) yields the following description of the boundary of P (its faces): b ( P ) = b ( S ∩ C ) = b ( S ) ∩ C + S ∩ b (C )
(2) If we go down another level of dimensionality, we find the boundaries of the faces, or edges. By applying the
boundary operator b() to eqn.(2) and discarding duplicated entities in the final result, we obtain:
We take advantage of an octree-based search and of Cartesian box intersection to limit the number of operations.
L = b( S ) ∩ b(C ) + b(b( S )) ∩ C + S ∩ b(b(C ))
Once the set of line segments L is found we perform a merge operation on the nodes of the segments using a small fraction of the initial cell size as tolerance. This step would not be necessary if we had used exact arithmetic [16] in the intersection operation. On the one hand exact arithmetic allows basic geometric algorithms to be built without local or global tolerances, whose usage may lead to catastrophic errors when the comparison of floating-point numbers is used at a logical branch point. On the other hand it has the potential to greatly slow-down the performance of the algorithm.
(3) where L is the set of edge segments bounding the faces of P. The intersection to be calculated consists of three parts: 1. 2. 3.
Intersection of the faces of the two solids S and C. The part of the edges of C that lies inside S. The part of the edges of S that lies inside C.
Part 2 and 3 of the intersection rely on the ability to determine if a point is inside or outside a solid. Several techniques are available to answer point classification queries for solids. A brief overview can be found in [10] (chap.22 pp. 9-10). For our purposes we have implemented a three-dimensional version of the so-called ray-casting algorithm. Basically the algorithm shoots a ray in a coordinate direction and counts the number of intersections with the triangulated surfaces. If the number of intersection is odd the point is inside, otherwise it is outside. In order to make the technique relatively insensible to small gaps and overlaps in the surface all three coordinates direction are used and a majority rule is applied to the answers. We have developed two algorithms to intersect the cells of the original mesh and the surface keeping in mind the need of a small computational cost. The first algorithm performs an exact intersection with a minimum of inside-outside calculations, which are typically the most expensive part of the intersection, while the second accepts minor approximations in the intersection as a prize to achieve further reductions in the computational cost. The availability of two algorithms with slightly different behaviors may be useful in case of failure of one of them. In the following sections brief descriptions of both are presented.
2.2 Cell-based algorithm A direct implementation of eqn.(3) leads to a cell-based algorithm, in which each cell of the initial mesh is processed separately and in the same way. Part 2 and 3 of the intersection do not present particular difficulties. They are based on the intersection of a line segment with a triangle, and an inside-outside calculation based on the ray-shooting technique. Part 1 of the intersection involves nc triangles of the boundary of C and ns triangles of the boundary of S. A naïve implementation of this intersection would generate an algorithm with an operation count proportional to ns x nc, where nc=24 and nc is of the order of at least several thousands for any geometrically interesting object. Clearly only a fraction of b(S) actually intersect b(C).
Once the set of all the line segments L=b(b(P)) has been found, a graph G is built with them. The edges of the graph represent the edges of the faces of the polyhedron P. The task here is to define a reconstruction operator t(L) b( P ) = t ( L ) = t (b(b( P )))
(4) that reconstructs the faces from their boundaries. It appears from eqn.(4) that the operator t() may be considered the inverse of b(): t (b( X )) = I ( X ) = X
(5) The target polyhedron, in order to be consistent, must have two, and only two, faces sharing an edge and each of its nodes must be at least trivalent. Having these requirements in mind, the edges of the graph are merged until all the nodes are at least trivalent. If we define the faces of a polyhedron as those closed loops in the graph with the shortest length, we can apply the socalled “greedy” algorithm from graph theory [17],[18] to find them. The “greedy” algorithm is able to find the path with the shortest distance between two nodes of a graph, where distance (or weight) actually represents any quantity that can be associated with an edge. A question that arises is: what definition of distance leads the greedy algorithm to find all and only the faces of the polyhedron? Even though we do not claim to have the general answer we found that the following definition is robust enough to find the faces of the polyhedral cells. The weight wij to be associated with an edge that connects node i and node j and belong to a path p is: wij = 1 + f (l j ) (6) where the function f(lj) is a measure of how much the closed loop lj is out-of-plane. The loop lj consists of the part of path p that goes from node 1 to node j plus the edge from the node j to node 1.
If we consider the unit vector uk connecting each node k with the node k+1 of loop lj and vk=uk x uk+1 then we can define:
Filling in the missing edges generates closed loops (faces). Additionally, if a face piercing is inside the closed loop, a generated edge is further subdivided to include it.
∑v
Volume Processing (3D): in this step, all the faces that define a volume inside the surface are found. The volume bounded by the mesh cell faces is divided into regions that are either outside or inside the surface. The outside regions are ignored. The inside regions are completed by constructing the faces required to close the volume. They are generated from edges that are defined in only one face. These edges form closed loops. Additionally, if surface corners are inside the cell, the closed loops are further subdivided to include them using a process similar to the face processing.
f (l j ) =
k
⋅ vk +1
k
(7) as the out-of-plane penalty for the loop lj.
2.3 Mesh-based algorithm A reduction of the computational effort required by the intersection between the initial mesh and the surface that is the boundary of the volume to be meshed (hereafter known just as the surface) can be achieved in two different ways: •
Take advantage of common components in adjacent cells.
•
Simplify some parts of the intersection.
The mesh-based algorithm starts by building a hierarchy of information about the intersections between the entire initial mesh and the surface. Information is captured and stored starting with the lowest dimensionality of the mesh (0D or corner points of each cell) and finally looking at the volume of the hexahedral cell. This approach serves two purposes. Since mesh details at lower dimensions are shared by many cells (in a structured mesh, a corner is shared by 8 cells, an edge by 4 cells and a face by 2 cells), it reduces the computational time by processing each component only once instead of once for each sharing cell as in the cell-based algorithm. Additionally, by using a common definition for each entity, mesh conformity is automatically preserved. The algorithm executes the following steps to evaluate the sidedness of the hierarchy of components: Corner Processing (0D): for each corner point of the mesh cells, determine its sidedness using the Inside-Outside rayshooting based procedure. Edge Processing (1D): the mesh cell edges are divided into segments that are inside the surface, outside the surface or on the surface. This is done using the information found in the previous step for the end points of the edge and calculating the intersections between the edge and the surface. The algorithm considers two different types of intersections: a surface crossing that indicates a transition of sidedness from inside to outside or vice versa and a segment that is entirely contained on the surface. For the latter case, the sidedness of the segment is set to boundary. Face Processing (2D): now, we take the faces of the mesh cells and divide them into areas that are inside the surface, outside the surface or on the surface. This process is done using the edge sidedness from the previous step and the information about the face interaction with the surface. In this case, the surface information needed is the location of the surface edges piercing the face. A set of edges is formed by retaining the inside part of the face edges.
Looking at these results, we see that information from each level of dimensionality of S and C is required to perform the intersection operation. The cells in the mesh can be looked upon as a probe of the surface S extracting information from each level of dimensionality. The following table summarizes the flow of information extracted from the surface. Mesh Entity
Previous Information
Mesh Entity Dimension
Operation With Surface
Surface Entity Dimension
Corner
None
0
Inside / Outside
3
Edge
Corner
1
Intersects Surface
2
Face
Edge
2
Pierced by Surface Edge
1
Volume
Face
3
Contains Surface Corner
0
“Previous Information” is the information extracted from the lower level operation that is needed during the current operation. “Operation With Surface” is the operation performed between the mesh entity and the surface and is needed to extract information at this level of mesh dimensionality. “Mesh Entity Dimension” and “Surface Entity Dimension” show, respectively, the dimensions of the information used and extracted in the current operation. For what concerns the simplification of the intersection we note that eqn.(3) can be formally modified by applying the two operator b() and t() one after the other to the three parts of the intersection: L = b( S ) ∩ b(C ) + b(b( S )) ∩ C + S ∩ b(b(C )) = = t[b(b( S ) ∩ b(C ) + b( S ) ∩ b(b(C ))] + + t[b(b(b( S ))) ∩ C + b(b( S )) ∩ b(C )] + + t[b(b(b(C ))) ∩ S + b(b(C )) ∩ b( S )] (8)
The modified eqn.(8) shows that, with appropriate reconstruction operators, the line segment set L may be derived from the following node information: •
The nodes resulting from intersecting the edges of S and the faces of C.
•
The nodes resulting from intersecting the edges of C and the faces of S.
•
The sidedness of the nodes of C respect to S.
•
The sidedness of the nodes of S respect to C.
As an example of reconstruction, let us concentrate on the second part of eqn. (3) b (b ( S )) ∩ C = t[b (b (b ( S ))) ∩ C + b (b ( S )) ∩ b (C )]
(9) that represents the part of the edges of S that are inside C. Let us suppose that the situation is the one described in Figure 1 From the R.H.S of eqn.(9) we have a corner p of S that is inside C and the nodes i1, i2 and i3 that are the intersections between the edges of S and the faces of C. Counting the number of edges attached to the internal corner p and noting that they are equal to the number of intersection between the edges of S and the faces of C it is possible to reconstruct the edges e1, e2 and e3 by simply connecting node p with nodes i1, i2 and i3. The simplification of the intersection has an interesting useful side effect. In fact if any of the edges e1, e2, or e3 is not continuous inside the cell C, i.e. there is a hole or an overlap in the surface S, still the algorithm would be able to reconstruct them. This means that the ‘mesh algorithm’ is relatively insensible to holes and overlaps smaller than the cell size.
Figure 1. Reconstruction of internal edges from node information
3. CELL GENERATION 3.1 Polyhedral cells In this section the actual cell generation is described. Given a set of faces that form a closed volume, the program finds the equivalent three-dimensional cell shape that the faces define by combinatorially comparing them with the faces of a list of acceptable cells. Many of these cells will be ‘traditional’ cells found in most FV and FE codes: tetrahedra, prisms, pyramids and hexahedra. Others will be more complex shapes that FV or FE programs may or may not handle. STAR-CD has been modified to recognize 6 basic polyhedral shapes (see Figure 2). However, it must be noted that STAR-CD accepts cells with collapsed (pinched) edges. From each of the six basic polyhedra shown in Figure 2, many other shapes can be derived by collapsing one or more edges. The edge collapsing operation is acceptable when it does not collapse any face. With this extension the number of acceptable polyhedral topologies is much larger. They can cover the overwhelming majority of the intersections. In the rare case of an unacceptable polyhedron, an automatic split in acceptable shapes is performed.
3.2 Generalized mid-point-subdivision Once all the faces of the polyhedral solid P have been found we may apply a generalization of the mid-pointsubdivision (MPS). Since MPS, in its standard formulation [11], can be used to split a convex polyhedron with all trivalent nodes, it is necessary first to decompose each concave polyhedron into two or more convex polyhedra. Existing general procedures for convex decomposition of polyhedra [19], when applied in the mesh generation framework, do not maintain the conformity of the mesh. On the contrary, they tend to create a large number of small convex pieces. Attempts to develop methods more related to mesh generation have been made [20], but we are not
Figure 2. Basic polyhedral cells accepted by STAR-CD
aware of techniques that guarantee the conformity of the resulting mesh when applied to an arbitrary non-convex set of polyhedra In the mesh generation framework here proposed, the task is simplified because, in practice, the polyhedra that can be generated by intersecting the initial hexahedral cell with the triangulated surface appear to be a very restricted class. Moreover, decreasing the size of the cells of the initial mesh tends to reduce the geometrical and topological complexity of the resulting polyhedra. An algorithm that performs a convex decomposition has been implemented with the following basic ideas: • • •
Split all concave faces by adding new edges passing through the concave points. Apply the “greedy” algorithm again to find new closed loops that can be built using the new edges. Use the new faces as splitting faces to cut the concave polyhedron and generate two or more new convex polyhedra.
Figure 3. Wireframe of a concave solid generated by intersecting a hexahedral cell and the surface
This algorithm preserves the conformity of the mesh because a concave face either is a boundary or is shared by two polyhedra. In both cases the splitting is unique and the conformity of the mesh is maintained. Let us now focus our attention on MPS as applied to a regular trivalent node. The node is attached to three other nodes through three edges and there are three faces sharing that node. Three other faces can be constructed by connecting the three mid-face nodes to the center of the polyhedron. These two sets of three faces form a hexahedron. An example of splitting a concave volume and a subsequent MPS is shown in Figure 3 and Figure 4. If the initial node is n-valent and we apply the same procedure the result is two sets of n faces, that we could call a 2n-hedron. For example, if n=4 we have an octahedron with quadrilateral faces, if n=5, a decahedron with quadrilateral faces, and so on. For these polyhedra there is no known valid hexahedral mesh so we propose the alternative of splitting them into pyramids maintaining a connected mesh by placing an additional point in the center of the 2n-hedron and connecting it with the quadrilateral faces. The result is a set of pyramids. The result is a generalization of MPS that is insensible to the valence of the nodes of the polyhedron. Since this method still divides each face of each polyhedron into quadrilaterals, a fully connected (conformity property) initial hexahedral mesh results in a fully connected hex-dominant mesh with some pyramids. If NP is the set of non-trivalent polyhedra created by the intersection, NN the set of non-trivalent nodes of the i-th polyhedron and vij is the valence of the j-th node of the i-th
Figure 4. Mesh in a concave solid after convex decomposition and MPS.
polyhedron, then the number of pyramids in the mesh is: n pyr =
∑ ∑ 2v
ij
i∈NP j∈NN
(10) The number of hexahedra is simply equal to the sum, for each polyhedron, of the number of trivalent nodes.
4. TEST CASES In this section we present three different test cases that show an increasing geometric complexity and have been successfully meshed. The first, a pyramid, is included because is an example of the generalized MPS. The second one is the geometry for an internal flow simulation of a valve, while the third is an internal combustion engine head that represent a class of geometries common in FE analysis. The meshes of the test
cases have not been smoothed, optimized, or projected onto the surface in order to show the result of the basic methodology.
Test Case Pyramid Valve Engine Head
Number of cells of initial mesh 27 35910 39688
Number of surface triangles
MPS
Number of cells
Total Time (sec.)
6 16864
Yes Yes
100 102178
10 220
31571
Yes
132133
325
Table 1 Figure 5. Mesh inside a pyramid. Table 1 summarizes the results for these test cases showing the number of cells in the initial meshes, the number of triangles in the surface, the number of final cells as well as the total elapsed time spent in the generation. An IBM IntelliStation with a Pentium II Xeon 450 Mhz and 256 MB Ram running Linux SUSE 6.0 has been used.
4.1 Pyramid A simple pyramid with quadrilateral base is presented. The edges of the base and the height have length 1. In Figure 5 the mesh for a pyramid is shown. Here the size of the initial mesh has been chosen to be to 0.4 and a generalized MPS has been applied to all the resulting polyhedra, including the internal hexahedra. The top of the pyramid has a quadvalent node. As a result there must be an octahedron split into 8 pyramids as shown in Figure 6.
Figure 6. The wireframe of the top octahedron plus its center.
4.2 Valve This test case represents the domain around a valve for an internal flow simulation. It has been chosen because, even if it is relatively simple from a geometric viewpoint, it presents highly concave regions and narrow passages between the valve and the boundary of the domain. See Figure 7, Figure 8 and Figure 9.
4.3 Engine head This test case has been chosen because it represents the class of meshing problems of structural analysis. The geometry presents holes, local concavities, convexities and some relatively thin parts. See Figure 10, Figure 11, Figure 12 and Figure 13.
Figure 7. Mesh of the Valve
Figure 8. Cross section of the mesh normal to the axis of the valve.
Figure 9. Detail of mesh close to the valve
5. CONCLUSIONS A method to generate hex-dominant meshes by intersecting an initial hexahedral mesh and a surface has been presented. The method, in principle, can be used to generate meshes for an arbitrarily complex threedimensional domain, preserving the significant surface details. It is also interesting to note that the method, in principle, may be applied in two dimensions and generate a truly allquad mesh. More basic research is needed to reliably perform the convex decomposition of concave solids. Such improvement will be useful both when using polyhedral cells and when splitting them with MPS. The meshes that can be created using the proposed methodology may have low quality cells and, if MPS is applied to eliminate the polyhedra, may not follow the surface curvature at the boundary as closely as possible. The first problem may be addressed by improving the mesh with an optimization-based smoother, while the second can be solved with a projection of the boundary nodes of the mesh onto the original surface and a subsequent surface mesh smoothing. Recent advances in both surface and volume mesh optimization [21],[22],[23],[24], let us believe that the quality of both a polyhedral mesh and a MPS mesh can be effectively improved. Other areas of improvements include local refinement preserving the conformity of the mesh, as shown in [8] and a more sophisticated adjustment procedure that could avoid the intersection of the edges of the initial cell with the sharp edges of the surface, further reducing the number of nontrivalent polyhedra and consequently of pyramids, if MSP is applied.
Figure 10. View of the surface mesh generated for the Engine Head test case.
Figure 11. Cross section of the mesh for the engine head
[6]
[7]
[8]
Figure 12. Detail of the mesh close to a hole in the engine head.
[9]
[10]
[11]
[12]
[13] Figure 13. Detail of the mesh of the engine head in a narrow passage.
REFERENCES [1] Tautges, T.J. and Mitchell, S., Progress report on the whisker weaving all-hexahedral meshing algorithm, Proceedings of the 5th Int. Conf. On Numerical Grid Generation in Computational Field Simulation, pp 659-670, 1996. [2] Blacker, T.D. and Meyers, R.J., Seams and wedges in plastering: a 3D hexahedral mesh generation algorithm, Engineering with computers, 9, pp 83-93,1993. [3] Mitchell, S.A., A characterization of the quadrilateral meshes of a surface which admit a compatible hexahedral mesh of the enclosed volume, Proceeding STACS ’96, Grenoble, 1996. [4] Thurston, B., Hexahedral decomposition of polyhedra, in Geometry in action: http://www.ics.uci.edu/~eppstein/gina/Thurstonhexahedra.html. [5] Li T.S., Armstrong C.G. and McKeag R.M., Automatic Partitioning of Analysis Models using
[14]
[15]
[16]
[17]
[18]
[19]
the Medial Axis Transform, Advances in Parallel and Vector Processing for Structural Mech, Sept 1994 pp 165-171. Sheffer, A., Etzion M. Rappoport A, Bercovier M., Hexahedral Mesh Generation using the Embedded Voronoi Graph, Proceedings of the 7th International Meshing Roundtable, Sandia National Laboratories, pp 347-364, 1998. Schneider, R., A grid-based algorithm for the generation of hexahedral element meshes, Engineering with computer, V.12, No.3-4, pp 168-177, 1996. Weiler F., Shindler R., Schneiders R., Automatic Geometry-Adaptive Generation of quadrilateral and Hexahedral Element Meshes for the FEM, Proceedings of the 5th Int. Conf. On Numerical Grid Generation in Computational Field Simulation, pp 689-697, 1996. Tchon, K., Hirtsch C. Schneiders R., OctreeBases Hexahedral Mesh Generation for Viscous Flow Simulations, Proceedings of the 13th AIAA Computational Fluid Dynamics Conference, Snowmass, CO, 1997. Thompson J.F., Soni B.K., Weatherill N.P. editors, Handbook of Grid Generation, chap.21 pp 7-18, CRC PRESS, ISBN 0-8493-2687-7, 1999. Li T.S., McKeag R.M. and Armstrong C.G., Hexahedral meshing using midpoint subdivision and integer programming, Computer Meth Appl Mech Engrg, 124, 171-193, 1995. Thompson J.F., Warsi Z.U.A., Mastin C.W., Numerical Grid Generation - Foundation and Application, North-Holland, ISBN 0-444-00985X, pp 310-315, 1985. Lorensen W.E. and Cline H.E. Marching cubes: A High Resolution 3D Surface Reconstruction Algorithm, ACM Computer Graphics, 21(4):163169,July 1987. Krishnan S. Narkhede A. Manocha D., BOOLE: A System to Compute Boolean Compinations of Sculptured Solids, Proceedings of the ACM Symposium on Computational Geomtery, 1995. Vanecek G Jr., Nau D.S. and Karinthi R.R., a simple Approach to performing Set Operation on Polyhedra, University of Maryland, Technical Research Report, TR 91-34. Edelsbrunner H. and Ernst P., Simulation of simplicity: A Technique to Cope with Degenerate cases in Geometric Algorithms, ACM Transaction on Graphics, 9(1):66-104, 1990. Kruse R.L., Bruce P., Tondo C.L., Data Structure and Program Design in C, Prentice-Hall, ISBN 0-13-725649-3, 1991. Dijkstra E.W., A note on two problems in connection with graphs, Numerische Matematik I, pp.269-271, 1959. Chazelle B. and Palios L., Convex partition of polyhedra: A lower bound and worst case optimal algorithm, SIAM J. Comp., 13(3):488-507,1984.
[20] Lu Y., Gadh R., Tautges J.T. Volume Decomposition and Feature Recognition for Hexahedral Mesh Generation, Proceedings of the 8th International Meshing Roundtable, South Lake Tahoe, 1999. [21] Castillo, J., A discrete variational method. SIAM, J. Sci. Stat. Comput., Vol. 12, No. 2 pp 454468,1991. [22] Zavattieri P., Optimization Strategies in Unstructured Mesh Generation, Intl. J. Num. Meth. Engr., Vol. 39, pp2055-2071,1996. [23] Ivanenko S.A., Harmonic Mappings, in Handbook of Grid Generation, pp 8-1,8-43, CRC Press 1999. [24] Knupp P., Matrix Norms & The Condition Number: A General Framework to Improve Mesh Quality Via Node-Movement, Proc.of the 8th Int. Meshing Roundtable, South Lake Tahoe, 1999.
2
[email protected], [email protected]
ABSTRACT A new methodology to generate a hex-dominant mesh is presented. From a closed surface and an initial all-hex mesh that contains the closed surface, the proposed algorithm generates, by intersection, a new mostly-hex mesh that includes polyhedra located at the boundary of the geometrical domain. The polyhedra may be used as cells if the field simulation solver supports them or be decomposed into hexahedra and pyramids using a generalized mid-point-subdivision technique. This methodology is currently used to provide hex-dominant automatic mesh generation in the preprocessor pro*am of the CFD code STAR-CD. Keywords: mesh generation, polyhedra, hexahedra, mid-point-subdivision, computational geometry.
1. INTRODUCTION In the past decade automatic mesh generation has received significant attention from researchers both in academia and in the industry with the result of considerable advances in understanding the underlying topological and geometrical properties of meshes. Most of the practical methodologies that have been developed rely on using simplical cells such as triangles in 2D and tetrahedra in 3D, due to the availability of relatively robust algorithms for such shapes. Nevertheless, in the industry, users of finite volume and finite element technologies in 3D tend to prefer hexahedral cells because in many cases hexahedra can better capture the characteristics of the solution and because a larger number of tetrahedral (arguably 3 to 5 times) is needed to achieve an equivalent accuracy in the solution. As a result, truly automatic unstructured hexahedral mesh generation is one of the most desired, and not completely fulfilled, features in a mesh generator. Several methodologies have been proposed to automatically create an unstructured hexahedral mesh. These methodologies can be ideally divided into two families: the first one starts the generation from a given closed surface tessellated with quadrilaterals and attempts
to create a mesh that satisfies the boundary constraint, while the second family does not require a quadrilateral surface mesh. Its creation becomes part of the generation process. The first family of methodologies such as spatial twist continuum STC and related whisker weaving [1], hexahedral advancing-front [2], rely on a proof of existence [3],[4] for combinatorial hexahedral meshes whose only requirement is an even number of quadrilaterals on the boundary surface. Although the methods of the first family have provided an enormous improvement in understanding the topological properties of hexahedral meshes, their application to reallife geometries has been limited by the mainly topological nature of the algorithms and by problems in ’closing’ the mesh i.e. in generating the last cells. In fact for several simple closed configurations of quadrilaterals, which should admit a hexahedral mesh, there is no known valid mesh. The second family of methodologies includes those based on medial surface [5],[6] and the so-called grid-based methods [7],[8],[9],[10]. The medial surface technique creates a two-dimensional skeleton of a three-dimensional volume as the locus of all
inscribed spheres of maximal diameter. The medial surface is then used to decompose the volume into parts meshable by existing techniques.
accommodate the solution of the PDE's. In fact all current FV methods use polyhedra. It happens to be that these polyhedra are hexahedra, prisms, pyramids or tetrahedra.
The so-called grid methods start covering the domain with a Cartesian hexahedral grid and removing from it the cells that either are outside the domain or are closer than a specified distance to the boundary. The boundary of this ‘initial mesh’ is made by quadrilaterals that can be extruded up to the surface, generating an isomorphic hexahedral mesh on the boundary.
Here we propose the use of the polyhedra that result from the intersection of an initial all-hex mesh that cover the whole domain of interest with a triangulated boundary surface. The result of the intersection or 'trimming' is the set of hexahedra of the initial mesh which lie totally inside the boundary surface plus a set of polyhedra generated by the intersection of the 'surface-related' hexahedra with the surface itself.
One of the problems of the medial-surface-based algorithms is the abundance of degeneracies in the medial surface leading to considerable topological complexity. The grid-based methods are relatively simple and robust. One limitation lies in the representation of the discontinuities (edges and corners) of the boundary surface. In fact it would be difficult for a mesh generated with such method to match a corner in the surface with an arbitrary number of sharp edges attached to it. A carefully handmade initial mesh could overcome the limitation, but this would also undermine the automatic character of the mesh generation.
1.1 Motivations The all-hex methods just mentioned, although very promising, may present some limitations when applied to the most extreme geometrical cases that can be encountered in the industry-related Computational Fluid Dynamics (CFD) and Stress Analysis (SA) simulations. Hexdominant algorithms have been proposed in an attempt to overcome the problem. The hex-dominant algorithms use one (or more) all-hex method but, when the all-hex method fails or in regions that are likely too complex for such method, then create more tractable shapes (typically tetrahedra, prisms and pyramids). In this paper we present a different approach to generate a hex-dominant mesh. If the constraint of using only hexahedra is released, it becomes acceptable to generate a hex-dominant mesh using any kind of shape or more generally a polyhedron, provided that the overwhelming majority of the mesh is still composed of hexahedra.
1.2 Polyhedra in the framework of Finite Volume and Finite Element methods Meshing may be defined as the process of breaking up a geometrical domain into smaller and geometrically simpler sub-domains (the cells) which should allow the solution of the discretized partial differential equations (PDE) that describe the behaviors of the fields involved in the physics of the phenomenon. In CFD one of the most popular ways of solving the Euler or Navier-Stokes equations is the Finite Volume (FV) method, which is based on the concepts of cell faces and flux through those faces. The method itself does not depend on the shape of the cells if the cell faces can be correctly defined. In this case a (convex) polyhedron may be used to
The method has a certain resemblance with the grid-method proposed in [7]. In fact it would create the same internal hexahedral structure, while the hexahedral buffer, which in Schneiders’ method is created using an isomorphic extrusion of the boundary faces of the interior mesh, here is replaced by the polyhedra created by the intersection. Since polyhedra with an arbitrary number of faces and arbitrary node valence (here valence is defined as the number of edges sharing a node) may represent arbitrary sharp edges and corners on the surface, the problem of a correct representation of the surface features which may affects Schneiders’ method is eliminated. On the one hand the direct usage of polyhedra as cells requires a solver able to handle them so that STAR-CD list of acceptable cell shapes has been extended to include a certain number of simple polyhedra. On the other hand, if an extension of the mid-pointsubdivision [11] is applied to the whole new mesh (including the polyhedra), as explained later, the result is a set of hexahedra and a much smaller set of pyramids. As a result, the proposed methodology has the potential to generate meshes for all the solvers that accept hexahedra and pyramids, i.e. the majority of Finite Elements (FE) and FV programs.
2. THE INITIAL MESH The initial mesh is a hexahedral mesh that must cover the whole domain of interest. The initial mesh may be structured or unstructured. Trivial ways to create the initial mesh include: •
Cartesian mesh generation
•
Cylindrical mesh generation (if the domain to be meshed presents an axis of symmetry).
•
Trans Finite Interpolation [12].
Moreover any combination of such techniques may be used to mesh different parts of the domain. We found that unmodified initial grids may lead to the generation of very small cells when the surface is close to one or more nodes of the cells. To overcome this difficulty we perform an adjustment procedure on the nodes of the initial grid.
We propose two different simple methods of accomplishing the task: 1.
Projecting the nodes of the initial mesh onto the surface when they are close to the surface. Here close means that the distance is smaller than a specified small fraction of the cell size.
2.
Moving the nodes of the initial mesh away from the surface when they are initially close to the surface.
To move a node far from the surface, first we find the vector from the node to its projection onto the surface, and then we move the node in the opposite direction up to a specified distance from the surface Another reason to move the mesh nodes close to the surface onto the surface or away from the surface is that the adjustment modifies the topology of the final mesh. In a polyhedron a trivalent node is a node that has exactly three edges and three faces attached, while for a nontrivalent node the number of edges and faces is not equal to three (typically it is greater). If all the nodes of a polyhedron are trivalent the polyhedron is trivalent. An interesting property of trivalent polyhedra is that when a b-rep solid represented by a trivalent polyhedron is intersected by another b-rep solid, the result is another trivalent polyhedron if the faces of the two solid intersect only at the edges and not at the corners, their sharp edges have no common point and there is no corner with more than three sharp edges inside the trivalent polyhedron. The proposed methods generate quite different final polyhedral cells. Method 1 tends to generate a larger number of polyhedra with non-trivalent nodes than the second one. In fact, the initial mesh after the adjustment has several nodes on the surface, and the intersection of a cell with a node on the surface with the surface itself may create a non-trivalent polyhedron. On the other hand it has the advantage of generating simpler polyhedra with a smaller number of faces and this is a desirable property because not all the polyhedra can be accepted as valid cell by STAR-CD (see 3.1.) If method 2 is applied, no node of the initial mesh is on the surface. The intersection creates a non-trivalent polyhedron only if the cell edges and the sharp edges of the surface intersect or if a surface corner with more that three sharp edges attached is inside the cell. This characteristic is important to maximize the number of the resulting hexahedra, when using the generalized mid-pointsubdivision (MPS), as we will explain later. The initial mesh does not need to have the property of conformity (for a definition of conformity see[8]) and an initial fully connected mesh may be locally and recursively refined according to several criteria including: •
Distance from the boundary
•
Curvature of the closest boundary
•
Interactive user input
A cell can be marked by one of the above criteria and split into 8 cells. All the information related to the cell neighbors is automatically updated. Other interesting improvements are under evaluation. For examples it would be useful to automatically create a roughly body-fitted initial mesh.
2.1 Basic Methodology We may consider the core of the polyhedral mesh generation as an intersection operation between two b-rep solids: the triangulated boundary surface S and the initial hexahedral cell C, whose faces, without loss of generality, can be triangulated in a unique manner choosing a point on each face (the face centroid) and connecting it to the four edges, creating four triangles per face. Once the triangulation is performed we have two b-rep solids whose boundaries are triangular faces. There are fast ways to evaluate approximated intersections such as the Marching Cube Algorithm[13], but it tends not to resolve the surface discontinuity leading to a poor representation of corners and sharp edges. Since it may be important to represent such details in the final mesh, we chose to perform an exact intersection. Given two polyhedral solids S and C, the intersection is a Boolean operation and the result P is yet another polyhedral solid: P = S ∩C (1) that is formed by: •
All components of S that lie inside C and all components of C that lie inside S.
There are several methods available in literature to perform Boolean operations on b-rep solids (a small bibliography can be found in [14]). However, the computational cost of several thousands or even hundreds of thousands of intersections between the cells of the initial mesh and the surface can be extremely high. An efficient implementation of the intersection operation is important for the practical applicability of this mesh generation method. As an extension of the boundary operator defined in [15], let us define b() as the operator that extracts the boundary of an entity, providing as result a set of one-dimensionlower entities. For completeness we may also define the identity operator as I(X)=X. Applying the b() to eqn.(1) yields the following description of the boundary of P (its faces): b ( P ) = b ( S ∩ C ) = b ( S ) ∩ C + S ∩ b (C )
(2) If we go down another level of dimensionality, we find the boundaries of the faces, or edges. By applying the
boundary operator b() to eqn.(2) and discarding duplicated entities in the final result, we obtain:
We take advantage of an octree-based search and of Cartesian box intersection to limit the number of operations.
L = b( S ) ∩ b(C ) + b(b( S )) ∩ C + S ∩ b(b(C ))
Once the set of line segments L is found we perform a merge operation on the nodes of the segments using a small fraction of the initial cell size as tolerance. This step would not be necessary if we had used exact arithmetic [16] in the intersection operation. On the one hand exact arithmetic allows basic geometric algorithms to be built without local or global tolerances, whose usage may lead to catastrophic errors when the comparison of floating-point numbers is used at a logical branch point. On the other hand it has the potential to greatly slow-down the performance of the algorithm.
(3) where L is the set of edge segments bounding the faces of P. The intersection to be calculated consists of three parts: 1. 2. 3.
Intersection of the faces of the two solids S and C. The part of the edges of C that lies inside S. The part of the edges of S that lies inside C.
Part 2 and 3 of the intersection rely on the ability to determine if a point is inside or outside a solid. Several techniques are available to answer point classification queries for solids. A brief overview can be found in [10] (chap.22 pp. 9-10). For our purposes we have implemented a three-dimensional version of the so-called ray-casting algorithm. Basically the algorithm shoots a ray in a coordinate direction and counts the number of intersections with the triangulated surfaces. If the number of intersection is odd the point is inside, otherwise it is outside. In order to make the technique relatively insensible to small gaps and overlaps in the surface all three coordinates direction are used and a majority rule is applied to the answers. We have developed two algorithms to intersect the cells of the original mesh and the surface keeping in mind the need of a small computational cost. The first algorithm performs an exact intersection with a minimum of inside-outside calculations, which are typically the most expensive part of the intersection, while the second accepts minor approximations in the intersection as a prize to achieve further reductions in the computational cost. The availability of two algorithms with slightly different behaviors may be useful in case of failure of one of them. In the following sections brief descriptions of both are presented.
2.2 Cell-based algorithm A direct implementation of eqn.(3) leads to a cell-based algorithm, in which each cell of the initial mesh is processed separately and in the same way. Part 2 and 3 of the intersection do not present particular difficulties. They are based on the intersection of a line segment with a triangle, and an inside-outside calculation based on the ray-shooting technique. Part 1 of the intersection involves nc triangles of the boundary of C and ns triangles of the boundary of S. A naïve implementation of this intersection would generate an algorithm with an operation count proportional to ns x nc, where nc=24 and nc is of the order of at least several thousands for any geometrically interesting object. Clearly only a fraction of b(S) actually intersect b(C).
Once the set of all the line segments L=b(b(P)) has been found, a graph G is built with them. The edges of the graph represent the edges of the faces of the polyhedron P. The task here is to define a reconstruction operator t(L) b( P ) = t ( L ) = t (b(b( P )))
(4) that reconstructs the faces from their boundaries. It appears from eqn.(4) that the operator t() may be considered the inverse of b(): t (b( X )) = I ( X ) = X
(5) The target polyhedron, in order to be consistent, must have two, and only two, faces sharing an edge and each of its nodes must be at least trivalent. Having these requirements in mind, the edges of the graph are merged until all the nodes are at least trivalent. If we define the faces of a polyhedron as those closed loops in the graph with the shortest length, we can apply the socalled “greedy” algorithm from graph theory [17],[18] to find them. The “greedy” algorithm is able to find the path with the shortest distance between two nodes of a graph, where distance (or weight) actually represents any quantity that can be associated with an edge. A question that arises is: what definition of distance leads the greedy algorithm to find all and only the faces of the polyhedron? Even though we do not claim to have the general answer we found that the following definition is robust enough to find the faces of the polyhedral cells. The weight wij to be associated with an edge that connects node i and node j and belong to a path p is: wij = 1 + f (l j ) (6) where the function f(lj) is a measure of how much the closed loop lj is out-of-plane. The loop lj consists of the part of path p that goes from node 1 to node j plus the edge from the node j to node 1.
If we consider the unit vector uk connecting each node k with the node k+1 of loop lj and vk=uk x uk+1 then we can define:
Filling in the missing edges generates closed loops (faces). Additionally, if a face piercing is inside the closed loop, a generated edge is further subdivided to include it.
∑v
Volume Processing (3D): in this step, all the faces that define a volume inside the surface are found. The volume bounded by the mesh cell faces is divided into regions that are either outside or inside the surface. The outside regions are ignored. The inside regions are completed by constructing the faces required to close the volume. They are generated from edges that are defined in only one face. These edges form closed loops. Additionally, if surface corners are inside the cell, the closed loops are further subdivided to include them using a process similar to the face processing.
f (l j ) =
k
⋅ vk +1
k
(7) as the out-of-plane penalty for the loop lj.
2.3 Mesh-based algorithm A reduction of the computational effort required by the intersection between the initial mesh and the surface that is the boundary of the volume to be meshed (hereafter known just as the surface) can be achieved in two different ways: •
Take advantage of common components in adjacent cells.
•
Simplify some parts of the intersection.
The mesh-based algorithm starts by building a hierarchy of information about the intersections between the entire initial mesh and the surface. Information is captured and stored starting with the lowest dimensionality of the mesh (0D or corner points of each cell) and finally looking at the volume of the hexahedral cell. This approach serves two purposes. Since mesh details at lower dimensions are shared by many cells (in a structured mesh, a corner is shared by 8 cells, an edge by 4 cells and a face by 2 cells), it reduces the computational time by processing each component only once instead of once for each sharing cell as in the cell-based algorithm. Additionally, by using a common definition for each entity, mesh conformity is automatically preserved. The algorithm executes the following steps to evaluate the sidedness of the hierarchy of components: Corner Processing (0D): for each corner point of the mesh cells, determine its sidedness using the Inside-Outside rayshooting based procedure. Edge Processing (1D): the mesh cell edges are divided into segments that are inside the surface, outside the surface or on the surface. This is done using the information found in the previous step for the end points of the edge and calculating the intersections between the edge and the surface. The algorithm considers two different types of intersections: a surface crossing that indicates a transition of sidedness from inside to outside or vice versa and a segment that is entirely contained on the surface. For the latter case, the sidedness of the segment is set to boundary. Face Processing (2D): now, we take the faces of the mesh cells and divide them into areas that are inside the surface, outside the surface or on the surface. This process is done using the edge sidedness from the previous step and the information about the face interaction with the surface. In this case, the surface information needed is the location of the surface edges piercing the face. A set of edges is formed by retaining the inside part of the face edges.
Looking at these results, we see that information from each level of dimensionality of S and C is required to perform the intersection operation. The cells in the mesh can be looked upon as a probe of the surface S extracting information from each level of dimensionality. The following table summarizes the flow of information extracted from the surface. Mesh Entity
Previous Information
Mesh Entity Dimension
Operation With Surface
Surface Entity Dimension
Corner
None
0
Inside / Outside
3
Edge
Corner
1
Intersects Surface
2
Face
Edge
2
Pierced by Surface Edge
1
Volume
Face
3
Contains Surface Corner
0
“Previous Information” is the information extracted from the lower level operation that is needed during the current operation. “Operation With Surface” is the operation performed between the mesh entity and the surface and is needed to extract information at this level of mesh dimensionality. “Mesh Entity Dimension” and “Surface Entity Dimension” show, respectively, the dimensions of the information used and extracted in the current operation. For what concerns the simplification of the intersection we note that eqn.(3) can be formally modified by applying the two operator b() and t() one after the other to the three parts of the intersection: L = b( S ) ∩ b(C ) + b(b( S )) ∩ C + S ∩ b(b(C )) = = t[b(b( S ) ∩ b(C ) + b( S ) ∩ b(b(C ))] + + t[b(b(b( S ))) ∩ C + b(b( S )) ∩ b(C )] + + t[b(b(b(C ))) ∩ S + b(b(C )) ∩ b( S )] (8)
The modified eqn.(8) shows that, with appropriate reconstruction operators, the line segment set L may be derived from the following node information: •
The nodes resulting from intersecting the edges of S and the faces of C.
•
The nodes resulting from intersecting the edges of C and the faces of S.
•
The sidedness of the nodes of C respect to S.
•
The sidedness of the nodes of S respect to C.
As an example of reconstruction, let us concentrate on the second part of eqn. (3) b (b ( S )) ∩ C = t[b (b (b ( S ))) ∩ C + b (b ( S )) ∩ b (C )]
(9) that represents the part of the edges of S that are inside C. Let us suppose that the situation is the one described in Figure 1 From the R.H.S of eqn.(9) we have a corner p of S that is inside C and the nodes i1, i2 and i3 that are the intersections between the edges of S and the faces of C. Counting the number of edges attached to the internal corner p and noting that they are equal to the number of intersection between the edges of S and the faces of C it is possible to reconstruct the edges e1, e2 and e3 by simply connecting node p with nodes i1, i2 and i3. The simplification of the intersection has an interesting useful side effect. In fact if any of the edges e1, e2, or e3 is not continuous inside the cell C, i.e. there is a hole or an overlap in the surface S, still the algorithm would be able to reconstruct them. This means that the ‘mesh algorithm’ is relatively insensible to holes and overlaps smaller than the cell size.
Figure 1. Reconstruction of internal edges from node information
3. CELL GENERATION 3.1 Polyhedral cells In this section the actual cell generation is described. Given a set of faces that form a closed volume, the program finds the equivalent three-dimensional cell shape that the faces define by combinatorially comparing them with the faces of a list of acceptable cells. Many of these cells will be ‘traditional’ cells found in most FV and FE codes: tetrahedra, prisms, pyramids and hexahedra. Others will be more complex shapes that FV or FE programs may or may not handle. STAR-CD has been modified to recognize 6 basic polyhedral shapes (see Figure 2). However, it must be noted that STAR-CD accepts cells with collapsed (pinched) edges. From each of the six basic polyhedra shown in Figure 2, many other shapes can be derived by collapsing one or more edges. The edge collapsing operation is acceptable when it does not collapse any face. With this extension the number of acceptable polyhedral topologies is much larger. They can cover the overwhelming majority of the intersections. In the rare case of an unacceptable polyhedron, an automatic split in acceptable shapes is performed.
3.2 Generalized mid-point-subdivision Once all the faces of the polyhedral solid P have been found we may apply a generalization of the mid-pointsubdivision (MPS). Since MPS, in its standard formulation [11], can be used to split a convex polyhedron with all trivalent nodes, it is necessary first to decompose each concave polyhedron into two or more convex polyhedra. Existing general procedures for convex decomposition of polyhedra [19], when applied in the mesh generation framework, do not maintain the conformity of the mesh. On the contrary, they tend to create a large number of small convex pieces. Attempts to develop methods more related to mesh generation have been made [20], but we are not
Figure 2. Basic polyhedral cells accepted by STAR-CD
aware of techniques that guarantee the conformity of the resulting mesh when applied to an arbitrary non-convex set of polyhedra In the mesh generation framework here proposed, the task is simplified because, in practice, the polyhedra that can be generated by intersecting the initial hexahedral cell with the triangulated surface appear to be a very restricted class. Moreover, decreasing the size of the cells of the initial mesh tends to reduce the geometrical and topological complexity of the resulting polyhedra. An algorithm that performs a convex decomposition has been implemented with the following basic ideas: • • •
Split all concave faces by adding new edges passing through the concave points. Apply the “greedy” algorithm again to find new closed loops that can be built using the new edges. Use the new faces as splitting faces to cut the concave polyhedron and generate two or more new convex polyhedra.
Figure 3. Wireframe of a concave solid generated by intersecting a hexahedral cell and the surface
This algorithm preserves the conformity of the mesh because a concave face either is a boundary or is shared by two polyhedra. In both cases the splitting is unique and the conformity of the mesh is maintained. Let us now focus our attention on MPS as applied to a regular trivalent node. The node is attached to three other nodes through three edges and there are three faces sharing that node. Three other faces can be constructed by connecting the three mid-face nodes to the center of the polyhedron. These two sets of three faces form a hexahedron. An example of splitting a concave volume and a subsequent MPS is shown in Figure 3 and Figure 4. If the initial node is n-valent and we apply the same procedure the result is two sets of n faces, that we could call a 2n-hedron. For example, if n=4 we have an octahedron with quadrilateral faces, if n=5, a decahedron with quadrilateral faces, and so on. For these polyhedra there is no known valid hexahedral mesh so we propose the alternative of splitting them into pyramids maintaining a connected mesh by placing an additional point in the center of the 2n-hedron and connecting it with the quadrilateral faces. The result is a set of pyramids. The result is a generalization of MPS that is insensible to the valence of the nodes of the polyhedron. Since this method still divides each face of each polyhedron into quadrilaterals, a fully connected (conformity property) initial hexahedral mesh results in a fully connected hex-dominant mesh with some pyramids. If NP is the set of non-trivalent polyhedra created by the intersection, NN the set of non-trivalent nodes of the i-th polyhedron and vij is the valence of the j-th node of the i-th
Figure 4. Mesh in a concave solid after convex decomposition and MPS.
polyhedron, then the number of pyramids in the mesh is: n pyr =
∑ ∑ 2v
ij
i∈NP j∈NN
(10) The number of hexahedra is simply equal to the sum, for each polyhedron, of the number of trivalent nodes.
4. TEST CASES In this section we present three different test cases that show an increasing geometric complexity and have been successfully meshed. The first, a pyramid, is included because is an example of the generalized MPS. The second one is the geometry for an internal flow simulation of a valve, while the third is an internal combustion engine head that represent a class of geometries common in FE analysis. The meshes of the test
cases have not been smoothed, optimized, or projected onto the surface in order to show the result of the basic methodology.
Test Case Pyramid Valve Engine Head
Number of cells of initial mesh 27 35910 39688
Number of surface triangles
MPS
Number of cells
Total Time (sec.)
6 16864
Yes Yes
100 102178
10 220
31571
Yes
132133
325
Table 1 Figure 5. Mesh inside a pyramid. Table 1 summarizes the results for these test cases showing the number of cells in the initial meshes, the number of triangles in the surface, the number of final cells as well as the total elapsed time spent in the generation. An IBM IntelliStation with a Pentium II Xeon 450 Mhz and 256 MB Ram running Linux SUSE 6.0 has been used.
4.1 Pyramid A simple pyramid with quadrilateral base is presented. The edges of the base and the height have length 1. In Figure 5 the mesh for a pyramid is shown. Here the size of the initial mesh has been chosen to be to 0.4 and a generalized MPS has been applied to all the resulting polyhedra, including the internal hexahedra. The top of the pyramid has a quadvalent node. As a result there must be an octahedron split into 8 pyramids as shown in Figure 6.
Figure 6. The wireframe of the top octahedron plus its center.
4.2 Valve This test case represents the domain around a valve for an internal flow simulation. It has been chosen because, even if it is relatively simple from a geometric viewpoint, it presents highly concave regions and narrow passages between the valve and the boundary of the domain. See Figure 7, Figure 8 and Figure 9.
4.3 Engine head This test case has been chosen because it represents the class of meshing problems of structural analysis. The geometry presents holes, local concavities, convexities and some relatively thin parts. See Figure 10, Figure 11, Figure 12 and Figure 13.
Figure 7. Mesh of the Valve
Figure 8. Cross section of the mesh normal to the axis of the valve.
Figure 9. Detail of mesh close to the valve
5. CONCLUSIONS A method to generate hex-dominant meshes by intersecting an initial hexahedral mesh and a surface has been presented. The method, in principle, can be used to generate meshes for an arbitrarily complex threedimensional domain, preserving the significant surface details. It is also interesting to note that the method, in principle, may be applied in two dimensions and generate a truly allquad mesh. More basic research is needed to reliably perform the convex decomposition of concave solids. Such improvement will be useful both when using polyhedral cells and when splitting them with MPS. The meshes that can be created using the proposed methodology may have low quality cells and, if MPS is applied to eliminate the polyhedra, may not follow the surface curvature at the boundary as closely as possible. The first problem may be addressed by improving the mesh with an optimization-based smoother, while the second can be solved with a projection of the boundary nodes of the mesh onto the original surface and a subsequent surface mesh smoothing. Recent advances in both surface and volume mesh optimization [21],[22],[23],[24], let us believe that the quality of both a polyhedral mesh and a MPS mesh can be effectively improved. Other areas of improvements include local refinement preserving the conformity of the mesh, as shown in [8] and a more sophisticated adjustment procedure that could avoid the intersection of the edges of the initial cell with the sharp edges of the surface, further reducing the number of nontrivalent polyhedra and consequently of pyramids, if MSP is applied.
Figure 10. View of the surface mesh generated for the Engine Head test case.
Figure 11. Cross section of the mesh for the engine head
[6]
[7]
[8]
Figure 12. Detail of the mesh close to a hole in the engine head.
[9]
[10]
[11]
[12]
[13] Figure 13. Detail of the mesh of the engine head in a narrow passage.
REFERENCES [1] Tautges, T.J. and Mitchell, S., Progress report on the whisker weaving all-hexahedral meshing algorithm, Proceedings of the 5th Int. Conf. On Numerical Grid Generation in Computational Field Simulation, pp 659-670, 1996. [2] Blacker, T.D. and Meyers, R.J., Seams and wedges in plastering: a 3D hexahedral mesh generation algorithm, Engineering with computers, 9, pp 83-93,1993. [3] Mitchell, S.A., A characterization of the quadrilateral meshes of a surface which admit a compatible hexahedral mesh of the enclosed volume, Proceeding STACS ’96, Grenoble, 1996. [4] Thurston, B., Hexahedral decomposition of polyhedra, in Geometry in action: http://www.ics.uci.edu/~eppstein/gina/Thurstonhexahedra.html. [5] Li T.S., Armstrong C.G. and McKeag R.M., Automatic Partitioning of Analysis Models using
[14]
[15]
[16]
[17]
[18]
[19]
the Medial Axis Transform, Advances in Parallel and Vector Processing for Structural Mech, Sept 1994 pp 165-171. Sheffer, A., Etzion M. Rappoport A, Bercovier M., Hexahedral Mesh Generation using the Embedded Voronoi Graph, Proceedings of the 7th International Meshing Roundtable, Sandia National Laboratories, pp 347-364, 1998. Schneider, R., A grid-based algorithm for the generation of hexahedral element meshes, Engineering with computer, V.12, No.3-4, pp 168-177, 1996. Weiler F., Shindler R., Schneiders R., Automatic Geometry-Adaptive Generation of quadrilateral and Hexahedral Element Meshes for the FEM, Proceedings of the 5th Int. Conf. On Numerical Grid Generation in Computational Field Simulation, pp 689-697, 1996. Tchon, K., Hirtsch C. Schneiders R., OctreeBases Hexahedral Mesh Generation for Viscous Flow Simulations, Proceedings of the 13th AIAA Computational Fluid Dynamics Conference, Snowmass, CO, 1997. Thompson J.F., Soni B.K., Weatherill N.P. editors, Handbook of Grid Generation, chap.21 pp 7-18, CRC PRESS, ISBN 0-8493-2687-7, 1999. Li T.S., McKeag R.M. and Armstrong C.G., Hexahedral meshing using midpoint subdivision and integer programming, Computer Meth Appl Mech Engrg, 124, 171-193, 1995. Thompson J.F., Warsi Z.U.A., Mastin C.W., Numerical Grid Generation - Foundation and Application, North-Holland, ISBN 0-444-00985X, pp 310-315, 1985. Lorensen W.E. and Cline H.E. Marching cubes: A High Resolution 3D Surface Reconstruction Algorithm, ACM Computer Graphics, 21(4):163169,July 1987. Krishnan S. Narkhede A. Manocha D., BOOLE: A System to Compute Boolean Compinations of Sculptured Solids, Proceedings of the ACM Symposium on Computational Geomtery, 1995. Vanecek G Jr., Nau D.S. and Karinthi R.R., a simple Approach to performing Set Operation on Polyhedra, University of Maryland, Technical Research Report, TR 91-34. Edelsbrunner H. and Ernst P., Simulation of simplicity: A Technique to Cope with Degenerate cases in Geometric Algorithms, ACM Transaction on Graphics, 9(1):66-104, 1990. Kruse R.L., Bruce P., Tondo C.L., Data Structure and Program Design in C, Prentice-Hall, ISBN 0-13-725649-3, 1991. Dijkstra E.W., A note on two problems in connection with graphs, Numerische Matematik I, pp.269-271, 1959. Chazelle B. and Palios L., Convex partition of polyhedra: A lower bound and worst case optimal algorithm, SIAM J. Comp., 13(3):488-507,1984.
[20] Lu Y., Gadh R., Tautges J.T. Volume Decomposition and Feature Recognition for Hexahedral Mesh Generation, Proceedings of the 8th International Meshing Roundtable, South Lake Tahoe, 1999. [21] Castillo, J., A discrete variational method. SIAM, J. Sci. Stat. Comput., Vol. 12, No. 2 pp 454468,1991. [22] Zavattieri P., Optimization Strategies in Unstructured Mesh Generation, Intl. J. Num. Meth. Engr., Vol. 39, pp2055-2071,1996. [23] Ivanenko S.A., Harmonic Mappings, in Handbook of Grid Generation, pp 8-1,8-43, CRC Press 1999. [24] Knupp P., Matrix Norms & The Condition Number: A General Framework to Improve Mesh Quality Via Node-Movement, Proc.of the 8th Int. Meshing Roundtable, South Lake Tahoe, 1999.
