FEM单元合并方法

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng (2010) Published online in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/nme.2963 A novel energy-based approach for merging finite elements Or Yogev, Andrew A. Shapiro and Erik K. Antonsson∗,† California Institute of Technology, Pasadena, CA 91125, U.S.A. SUMMARY A novel approach for merging two intersecting finite elements is presented and demonstrated. The solution mimics concepts from biology and uses principles rooted in continuum mechanics. The problem of attaching (or merging) two coincident finite elements is common when using the plastering technique as part of the advancing front method. This problem is particularly challenging for 3-D meshes of non-convex shapes. Some automatic meshing methods require portions of the partially formed mesh to coincide and merge. This problem is generally solved with heuristic rules, which lack generality, and may have difficulties with unforeseen situations. The problem of merging two overlapping polyhedra may also appear in other applications such as computer graphics and CAD software. A new approach to address the problem of merging is presented here. This solution does not utilize heuristic rules, but rather uses an approach based on minimization of strain energy. A fully automatic merging routine has been created that can address, in an optimum way, any situation of two nearby or overlapping elements that are to be merged. This approach, with minor adjustments, is suitable for most types of 3-D elements. Copyright  2010 John Wiley & Sons, Ltd. Received 9 June 2009; Revised 18 May 2010; Accepted 25 May 2010 KEY WORDS: meshing; element merging; morphogen; strain energy; energy minimization; 3-D 1. INTRODUCTION The problem of attaching (or merging) two coincident finite elements is common when using the plastering technique as part of the advancing front method [1, 2]. Generation of a mesh for non- convex shapes utilizing the advancing front method usually starts at multiple locations along the boundary and proceeds towards the center. At some point during the process, each of the individual portions of the mesh must coincide and merge. This problem, from a geometric viewpoint, is complex because the adjacent element orientations are generally not able to be predicted, and it is difficult to determine a priori the best configuration for merging adjacent elements. As a result, this type of problem is generally solved with heuristic rules or a lookup table approach [3]. These ∗Correspondence to: Erik K. Antonsson, California Institute of Technology, Mail Stop: 104-44, 1200 East California Blvd., Pasadena, CA 91125, U.S.A. †E-mail: erik.antonsson@caltech.edu Copyright  2010 John Wiley & Sons, Ltd. O. YOGEV, A. A. SHAPIRO AND E. K. ANTONSSON rules handle different situations in different ways and are primarily determined by the skill and experience of the programmer. The disadvantages of using heuristics include lack of generality, and difficulties with unforeseen situations. A new approach to address the problem of merging is introduced here. This solution does not utilize any heuristic rules. Instead it uses a natural approach based on minimization of strain energy that can automatically merge two cells into an optimum merged configuration. Other processes that are usually integrated with the advancing front method to generate good quality meshes (such as healing techniques [4] and mesh quality control [5, 6]) are integrated into the new energy minimization merging approach presented here. Meshed regions can consist of a large number of elements, so a robust and computationally efficient method for merging elements is critical. The method introduced here produces valid meshes for all cases, and utilizes an analytical expression for the conjugate gradient method to minimize computational cost. 1.1. Prior related work Despite considerable work, the problem of generalizing meshes for arbitrary non-convex three- dimensional bodies remains an open area of research. Excluding mesh-free methods, 3-D meshes generally comprise either tetrahedral or hexahedral elements. In the tetrahedral case, three popular methods have been developed for automated mesh generation for three-dimensional bodies: octree, Delaunay triangulation, and advancing front. The octree method was developed in the 1980’s by Mark Shephard [7, 8]. The idea behind the octree method is that a given geometrical domain is subdivided into cubes that are successively subdivided until a desired resolution has been achieved. New methods have been recently developed to accelerate an octree-based mesh generation method using parallel algorithms and clusters of many computer processors [9, 10]. The main shortcoming of the octree method is that a pre-defined surface or volume is weakly matched to its representative mesh: the surface of the mesh is not conformal. Instead, more elements are formed near the boundary in order to refine the mesh to match the boundary [11]. On the contrary, the Delaunay triangulation method matches the body surface more accurately [12–14]. This method is based on the Delaunay criterion, which tests a triangulated set of vertices such that any circumsphere that predefines the mesh does not contain any other vertices [15]. This criterion is used in generating tetrahedral meshes for 3-D shapes [16–19]. Although the Delaunay triangulation method can generate a mesh for any given convex body defined as a convex hull, it usually fails to generate meshes for non-convex domains [20], although some interesting attempts to deal with this impediment have demonstrated partially successful results [17, 21]. The advancing front (AF) method deals with the issues of non-convexity successfully and today is becoming a common method for generating meshes for both convex and non-convex bodies. In the AF method, the tetrahedral or hexahedral elements are generated toward the interior part of a domain from its pre-meshed surface [2, 22]. This method is analogous to a wave transmitted in the body. In the case of tetrahedral elements, the algorithm first generates a triangular base on each facet of the front, and then optimizes the location of a fourth node. The fourth node might be a new or an existing node of an adjacent element. The algorithm keeps running a collision detection sub-algorithm to ensure that no elements are overlapping. The advancing front method has also been used with meshes of hexahedral elements [1, 3, 23, 24]. Copyright  2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme A NOVEL ENERGY-BASED APPROACH In a similar way to the tetrahedral case, the plastering method begins with predefined meshed boundaries that advance the fronts inwardly into the domain [25–28]. As the fronts collide, multiple geometric operations are performed to smooth and seam the elements in order to minimize the remaining void of the un-meshed area. By the end of this process a new element is inserted into the void, which completes the meshing procedure. This algorithm suffers from many deficiencies when the fronts collide because of the complexity of dealing with the remaining void. One of the solutions to this problem is to use a tetrahedral mesh for the remaining void and then to transform the tetrahedral elements into hexahedral ones [29, 30]. The drawback of this method is the poor quality of the hexahedral elements generated during the transformation step. Acknowledging the complexity of hexahedral elements, the Q-Morph method starts with a triangulated surface mesh, which is transformed into hexahedral elements as the front progressively moves in the solid body [31]. The final shape of the quadrilateral mesh will be on the same scale as the initial surface mesh. The Q-Morph method also avoids expensive intersection calculations which are common to the advancing front method. Similarly, the H-Morph method starts with a tetrahedral mesh, which is transformed into hexa- hedral elements as the front (initially a set of prescribed quadrilateral surface facets) moves progressively into the volume [3]. The H-Morph method is able to generate a mesh without the need to identify or process special element geometries. The whisker weaving method is a purely topological mesh approach, defined as the geometrical dual of a mesh as a collection of surfaces, called the Spatial Twist Continuum (STC) [32, 33]. This method starts with a pre-quadrilateral surface mesh and then constructs a hexahedral connectivity mesh by advancing into the solid. Though improved, the whisker weaving approach has only been successfully used to generate meshes of relatively simple shapes [34]. Utilizing the simplicity and reliability of the plastering method, while avoiding the need to mesh complex voids, the unconstrained plastering method shows promising results in generating 3-D hexahedral meshes [24]. The unconstrained plastering begins with un-meshed surfaces where new hexahedral elements are formed by the intersection of the projected fronts of the surfaces, as they inwardly progress into the body. The unconstrained plastering is not a finished concept yet, but rather a work in progress [35]. 1.2. New method This paper provides a novel method to deal with the challenge of the plastering method to create hexahedral element connectivity when several fronts intersect. Our approach is motivated by the challenge of attaching (or merging) two coincident finite elements in an optimum way as part of the plastering approach. The problem of merging two overlapping polyhedra may also appear in other applications such as computer graphics and CAD software. The common solution for addressing merging of finite elements is to create a set of heuristic rules [3]. These rules test and categorize states of elements according to their relative position and orientation. Each state corresponds to a process that merges elements in a predefined manner according to a lookup table. The use of heuristic rules creates many problems. First, it requires a long set of conditional statements that introduce complexity into the code. Additionally, the rules must address and capture all situations regarding the relative position and orientation of the elements, and decide which kind of merging process to execute. Since the number of possible orientations is large, it is difficult to ensure that all possible situations are addressed. This is the motivation for developing a new method that automates the merging process between two finite Copyright  2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme O. YOGEV, A. A. SHAPIRO AND E. K. ANTONSSON elements regardless of their orientations. In this paper a novel approach is presented that mimics merging phenomena observed in nature, such as the merging of biological cells and bubbles. By using the basic principles used in nature, a fully automatic merging routine has been created that can address any situation of two nearby or overlapping elements that are to be merged. 1.3. Overview In Section 2 the idea of formulating each finite element as an artificial cell containing some particular properties is introduced. The principle of minimum potential energy, as observed in bubble merging processes, is also established and explained. Section 3 describes the governing equations for applying this method to an eight-node hexahedral linear brick finite element. Section 4 presents the derivation of the conjugate gradient method used here to minimize strain energy in the configuration of the merged elements. An analytical expression for the strain energy gradient has been derived, which significantly increases the speed of convergence. Section 5 discusses the computational cost of merging elements using this approach and compares it with other methods. Section 6 shows the results of the merging of two elements using this method. Several examples illustrate the capabilities of this approach in automating the merging process without using heuristic rules in Section 7. Conclusions are presented in Section 8. 2. APPROACH Consider two finite elements (such as those shown in Figure 1(a)) that under certain conditions need to go through a merging process. Our approach is independent of the pre-condition setting, which can be the distance between center–center, face–face, node–face or intersection [36] of the paired elements. In this work we use an intersection condition to illustrate the technique as we believe it is the most challenging condition. The merging process takes one face from each element and merges their corresponding nodes. The method introduced here comprises three steps: first, determine the corresponding faces to be merged; second, determine the pairs of nodes (one from each element) to be merged; and third, optimize the positions of the merged nodes. Cells in nature, prior to merging, communicate by way of a chemical signaling mechanism [37] which activates the operation of other genes, regulated by the level of the signal exceeding a threshold [38]. This process is mimicked here by considering each finite element to be the analog of a cell. Each cell generates a communication signal at its center, and this signal propagates through the faces of the cell and into the adjacent cells, schematically illustrated in Figure 1(a). Each signal is characterized by an exponential function (1) that decays with distance from the center of the cell. Equation (1) approximates the behavior of a steady-state diffusion equation (as occurs in cell signaling). Although it is not a rigorous solution, solving the diffusion equation more accurately will not have any effect on the performance of this method (but will slow it down). Equation (1) is used to describe this process: S =exp(−xc) (1) where  is a constant, xc is the vector from the coordinates of the center of the cell, and S is the gradient vector of the concentration (or strength) of the signal. Copyright  2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme A NOVEL ENERGY-BASED APPROACH 1 2 3 0 4 5 6 7 1 2 3 0 4 5 6 7 Intersection points(a) (b) (c) Figure 1. (a) Two cells prior to merging with radiating signals from the center of each cell; (b) the faces of the two cells to be merged; and (c) the cells after merging. The face of each cell that contains the maximum intensity of the signal from the adjacent cell is chosen for merging, as illustrated schematically in Figure 1(b). Inverse isoparametric mapping [39] is then applied, which maps each cell to its local configuration. By inverse mapping the signal as well, the detection of the point with the maximum signal gradient on the adjacent face becomes straightforward. Once the faces to be merged have been identified, the best four pairs of nodes to be merged must be selected. Each pair of nodes will be merged into a single node initially positioned at the mean of the coordinates of both original nodes. At the end of the process, the two cells merge to each other without any overlap, as illustrated in Figure 1(c). For two hexahedral elements, there are a total of 24 different configurations of merged node pairs from which to choose. The selection of the optimum configuration mimics the process of bubble merging [40]. When bubbles merge, they adopt a shape that reduces their total merging surface energy. This principle minimizes the distortion of each bubble from a spherical shape. Using this energy minimization principle for all 24 possible combinations of nodes to merge, the pairing of nodes that results in the minimum distortion to both cells, and thus introduces the minimum change in strain energy (when each pair of nodes is merged into a single node positioned at the mean of the original coordinates of the pair) is selected. Any contribution to the change in strain energy will originate only from the distortion of the cells. This idea is applied rigorously here by considering the constitutive properties of the materials of each cell to be quasi neo-Hookean, but only for the purpose of merging the cells. Copyright  2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme O. YOGEV, A. A. SHAPIRO AND E. K. ANTONSSON Prior to the merging process, both cells are in a reference configuration. While the merging process is being applied, a finite displacement constrains each pair of nodes. The strain energy density relationship for neo-Hookean materials is given by w in (2). Note that Equation (2) is not the general neo-Hookean relation but rather a representation of incompressible material, nevertheless this relation is sufficient for the technique. w=C1(I1 −3) (2) where C1 is a constant, and: I1 = C ′ · I = ( F ′T F ′)· I (3) F ′ = det( F)−1/3 F (4) F = x  X = ( X +u)  X = I + u  X (5) where u is the displacement, x is the new coordinate, X is the initial coordinates of the cell, C is the Cauchy stress tensor, F ′ is the distortional part of the deformation gradient, and F is defined in (5) such that det( F ′)=1 [41]. This is the pure measure of distortion. The deformation gradient F is composed of two parts, the identity tensor, I , and the derivatives of the displacements with respect to the reference configuration. The displacements of the nodes are known, and the derivatives can be computed using the isoparametric mapping [23, 42]. Once the strain energy density w is determined, the total strain energy can be computed by integrating w over the volume of the cell as shown in (6): U = ∫  wdV (6) Out of the 24 possible combinations, the one that produces the minimum strain energy U is chosen. At the end of this process, the two cells are merged with no overlap, and the strain energy required for the operation is the minimum of all 24 cases. Out of all 24 combinations, some may result in non-convex cells. This represents an ill- conditioned case where the one-to-one isoparametric mapping breaks down. A way to overcome this problem is to condition the determinant of the Jacobian matrix, evaluated at the Gauss quadra- ture points, as part of the energy calculation process. An energy estimator will be assigned into particular state if, and only if, the condition presented in Equation (7) is satisfied. In case the latter condition fails, the resultant energy will receive the value MAX ENERGY and thus will prevent this state from being selected. det( Ji )>0 ∀ i =1,2, . . . ,8 (7) Once a state is been selected it can still be improved, in terms of energy minimization. The initial assumption of merging a pair of nodes by taking the average of their corresponding coordinates will not generally result in a global energy minimization for the two merged cells. In order to Copyright  2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme A NOVEL ENERGY-BASED APPROACH improve the initial assumption, the conjugate gradient method is applied to the displacements to find the global minimum of the strain energy, as shown in (8), where u1k are the displacements of the nodes corresponding to the selected face of the first cell, and u2k is defined in the same way for the selected face of the second cell. The iteration of umk j , j0 in the conjugate gradient method satisfies the recursion shown in (9) [43, 44]. min(U (u1k +u2k) : u ∈ R3) k =1, . . . ,4 (8) umk j+1 = umk j + j d j (9) The directions d j are generated by the rule in (10): d j+1 =−g j+1+ j d j (10) where d0 = −g0 g j = ∇U (umk j ) (11) The step size  j is determined by a global line search developed by Fletcher–Reeves [45], shown in (12):  j =− gTj · d j dTj · g′j (12) An additional constraint was developed by Wolf [46], shown in (14), with the inequality (13): 0<<<1 (13) where  and  are two scalars, and  is determined according to Fletcher–Reeves algorithm (15): U (umk j )−U (umk j + j d j )− j (gTj ·