Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976 www.elsevier.com/locate/cma

3D tetrahedral, unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric Cyril Gruau *, Thierry Coupez Material Processing Center (CEMEF), E´cole des Mines de Paris, BP 207, 06904 Sophia-Antipolis, France Received 28 November 2003; received in revised form 25 June 2004; accepted 26 November 2004

Abstract In this paper we present a 3D tetrahedral, unstructured and anisotropic mesh generator that is not based on the Delaunay, frontal or octree method. Instead, it proceeds by local optimizations and uses an anisotropic shape criterion to ﬁt a metric ﬁeld. Then, we introduce a new 3D metric ﬁeld that tightens the mesh around interfaces when the calculation domain is divided in several subdomains, and a 3D metric ﬁeld that places enough elements through each subdomain thickness, without introducing too many nodes in the other directions. Finally, we show some applications for material forming geometries. 2005 Elsevier B.V. All rights reserved. Keywords: Topology and shape optimization; Elliptic interpolation; Thickness detection and curvature treatment; Interface reﬁnement

1. Introduction Recent advances in numerical methods together with high performance computing made it possible to simulate 3D applications involving several objects (i.e. several subdomains). In material processing, these multidomain simulations may concern injection moulding process [1], foam expansion [2], forging process [3] or composite forming [4]. In this context, numerical software has to treat increasingly complex three-dimensional geometries, especially with industrial thermoplastic parts [5] and with the casting process [6]. Today, the great challenge is to

*

Corresponding author. E-mail addresses: [email protected] (C. Gruau), [email protected] (T. Coupez).

0045-7825/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.11.020

4952

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

generate suitable meshes for computation, without much human intervention. In our case, a suitable mesh means: • a mesh with several element layers everywhere through the thickness (to observe non-Newtonian behavior and thermal eﬀects, for example), even for thin and curved geometries, • a mesh that captures precisely the interfaces between subdomains (between ﬂuid and structure for example), • a mesh that does not contains too many nodes (because of memory capacity and computational time limitations). Therefore, we tend to generate anisotropic meshes, because isotropic meshes would contain too many nodes for an implicit three-dimensional ﬁnite element solver. So, we need a mesh generator that can take a metric ﬁeld into account. Furthermore, we have to deﬁne a metric ﬁeld that meets the common needs detailed above. The originality of this paper is to propose an alternative to the existing 3D, unstructured and anisotropic mesh generators [7–10], and to present two kinds of 3D metric ﬁelds that are neither produced by a posteriori information [11], nor designed for boundary layer improvement [12], i.e.: • the natural metric, that places enough element layers through the thickness while saving maximum nodes in the other directions (see a natural mesh, Fig. 1(a)), • and the multidomain metric, that gathers various subdomains in a single mesh (see a multidomain mesh, Fig. 1(b) and (c)).

Fig. 1. A natural mesh of an optical lens quarter (2644 nodes), a bidomain mesh of this lens in its mould (22,893 nodes) and a sectional view.

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4953

This paper is organized as follows: In Section 2, we introduce notations and deﬁnitions to explain our anisotropic mesh generation process. In Section 3, we deﬁne the natural metric, we indicate a way to compute it and we show some results. In Section 4, we do the same with the multidomain metric. Finally, concluding remarks are gathered in Section 5.

2. Mesh generation and anisotropic adaptation by topological optimization In this section, an adaptive method for generating anisotropic and unstructured meshes in any dimension is described. Firstly, useful properties about meshes and mesh topologies are established (Sections 2.1 and 2.2). Secondly, Section 2.3 details a mesh generation by topological optimization, which was ﬁrst introduced in [13] and reused by [14]. Thirdly, we expose how this adaptive method can be used to conform a mesh to a nonEuclidian metric ﬁeld (Section 2.4). From now, our calculation domain is a polytope X 2 Rd (with d = 2, 3, . . .), which is bounded, orientable and with non-empty interior. 2.1. Basic deﬁnitions and properties For the sake of completeness, the purpose of this paragraph is to deﬁne the classical notions that are widely used in this paper. Deﬁnition 1. A metric is a d-rowed and d-columned, real, symmetrical, positive deﬁnite matrix. In fact, according to the following proposition, a metric M consists in d orthogonal directions (the columns of R) and one size hi in each direction. Proposition 2. A metric M can be diagonalized in 3 2 1 0 7 6 h2 7 6 1 7 > 6 .. 7R ; M ¼ R6 7 6 . 7 6 4 15 0 h2d

ð1Þ

where R is a rotation matrix whose columns are eigenvectors of M and ðh2 i Þ16i6d the corresponding eigenvalues. Proof. Since M is symmetrical, M can be diagonalized in the orthogonal group (at this stage, R = R>). Since M is positive deﬁnite, the eigenvalues k of M are strictly positive (so, $ h 5 0 such that k = 1/h2). For the columns of R, we can choose unit eigenvectors of M (then, jdet(R)j = 1). If det(R) = 1, we can choose the opposite of one eigenvector so as to get det(R) = 1 (hence, R is a rotation matrix). h Before deﬁning the requirement for a mesh, we need the deﬁnition of a simplex in any dimension and with any number of vertices. Deﬁnition 3. A k-simplex of Rd (with 0 6 k 6 d) is the convexhull of k + 1 points in Rd (called vertices) and a k-face of a simplex T is a k-simplex whose vertices are vertices of T. It is usual to say that a (k 1)-face of a k-simplex is shortly a face and a (1)-face is the empty set. Now, in many steps the calculus of the volume of a simplex is required.

4954

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

Proposition 4. Let T be k-simplex of Rd with vertices S0, . . . , Sk. We denote by A the rectangular matrix (S0S1, . . . , S0Sk). The volume of T is pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ detðA> AÞ jT j ¼ . ð2Þ k! So, we can compute the volume of a simplex for each dimension d and each simplex type k. Formula (2) can be rewritten as • when k = 1

jT j ¼ kS 0 S 1 k, j detðAÞj , • when k = d jT j ¼ d! kS 0 S 1 ^ ^ S 0 S d1 k • when k = d 1 jT j ¼ , ðd 1Þ! where ^ ^ stands for the cross product in Rd . Proof. Let Tb be the k-simplex of Rk whose vertices are (dij)16i,j6k, one can easily prove that j Tb j ¼ k!1 . We denote by f the linear function such that f ð Tb Þ ¼ T and $f = A. From Riemannian geometry [15,16], we know that Rd (provided with its canonical Euclidian metric) induces on T (with its map f) the metric pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ R pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ R pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ M = A>A and that jT j ¼ T^ detðMÞ dx. Since M is constant we have T^ detðMÞ dx ¼ j Tb j detðMÞ. Hence, we obtain Formula (2). h A simplex T is said to be non-degenerated if and only if jTj > 0. Now, we are able to deﬁne a mesh. Deﬁnition 5. Let N be a ﬁnite node set in X, let T be a set of non-degenerated d-simplices whose vertices belong to N and let F be the faces of those simplices. ðN; TÞ is a mesh of X if and only if: (i) S each face of F does not share more than two simplices of T, (ii) T 2T T ¼ X, (iii) for all simplices T 1 6¼ T 2 2 T, the intersection T1 \ T2 is a k-face of both simplices with k 2 {1, 0 , . . . , d 1}. The (topological) boundary of a mesh ðN; TÞ is oT, the faces of F that belong to only one simplex in T. We have the same deﬁnition for a mesh ðN; UÞ of oX where U is a set of (d 1)-simplices (which is used in Deﬁnition 7). It is obvious that if ðN; TÞ is a mesh of X, then ðN; oTÞ may not be a mesh of oX, because condition (i) may not be satisﬁed. However, the volume of X can be computed with oT alone. Proposition 6. Let ðN; TÞ be a mesh of X. The volume of this mesh is X jXj ¼ jT j;

ð3Þ

T 2T

Z 1 ¼ d S

x ~ n dx F 2oT

ð4Þ

F

with ~ n the unit normal vector pointing outward. Owing to Formula (4) the volume of a mesh can be computed using only its boundary. R R Proof. We have jXj ¼ X dx ¼ S T dx. Since simplices of T do not overlap themselves, we have T 2T RS R P P dx ¼ dx ¼ jT j. T 2T

T

T 2T T

T 2T

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4955

Besides, according theorem [17], as the divergence $ Æ x =S d and since X is an orientable R to Stokes R R S polytope, we have X dx ¼ d1 X r x dx ¼ d1 oX x ~ n dx. Finally, since X ¼ T 2T T , we have oX ¼ F 2oT F , which leads to Formula (4). h 2.2. The mesh topology deﬁnition and the minimal volume criterion Since a mesh is a very diﬃcult object to build directly, especially when d P 3, we introduce an easier object (a mesh topology) in the following deﬁnition. The purpose of this section is to rise a suﬃcient condition to obtain a mesh from a mesh topology. Deﬁnition 7. Let N be a ﬁnite set of nodes in X, let T be a set of d-simplices whose vertices belong to N and let F be the faces of those simplices. T is a mesh topology if and only if: (i) each face of F does not share more than two simplices of T, (ii) ðN; oTÞ is a mesh of oX. Again, in this deﬁnition, the boundary oT of a mesh topology T is the faces of F that belong to only one simplex in T. We have the same deﬁnition for a mesh topology ðN; UÞ of oX where U is a set of (d 1)-simplices (which is used in Lemma 10). R Thanks to condition (ii), we still have jXj ¼ d1 S F x ~ n dx. However, all mesh topologies are not F 2oT

meshes because their simplices can be degenerated, overlap themselves or overﬂow X. Hence, with this P definition we do not have jXj ¼ T 2T jT j but only the following result. Lemma 8. Let T be a mesh topology, then X

S

T 2T T .

That is, thanks to condition (ii) of Deﬁnition 7, each point of the calculation domain X lies inside an element of T. Proof. We consider the following advancing front: • C0 ¼ oT and T0 ¼ ; • 8n 2 N – if Cn = ; then Cn+1 = Cn – if Cn 5 ; then Cn+1 = Cn 4 o{Tn} (with the operator A 4 B = (A [ B)n(A \ B)) and Tnþ1 ¼ Tn [ fT n g where T n 62 Tn is such that Cn \ o{Tn} 5 ;. Since T is a ﬁnite set, 9 N 2 N such that CN1 5 ; and "n P N Cn = ;. Now, it can be proved that 80 6 n < N Cn ¼ oðT n Tn Þ: • With n = 0, we have C0 ¼ oT ¼ oðT n T0 Þ because T0 ¼ ;. • Now with 0 6 n < N, if Cn ¼ oðT n Tn Þ then let F 2 Cn+1. – If F 2 Cn, then F 62 o{Tn} (because Cn+1 = Cn 4 o{Tn}) so, 9 !T 6¼ T n 2 T n Tn such that F 2 o{T}. Hence, F 2 oðT n Tn n fT n gÞ ¼ oðT n Tnþ1 Þ. – If F 62 Cn, then F 2 o{Tn} and F 62 oT so, condition (ii) of the Deﬁnition 7 gives 9 !T 6¼ T n 2 T such that F 2 o{T}. We have T 62 Tn , otherwise F 2 oðT n Tn Þ ¼ Cn . Hence, F 2 oðT n Tn n fT n gÞ ¼ oðT n Tnþ1 Þ. So, Cnþ1 oðT n Tnþ1 Þ. Now, let F 2 oðT n Tnþ1 Þ ¼ oðT n Tn n fT n gÞ. We have 9 !T 2 T n Tn n fT n g such that F 2 o{T}.

4956

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

– If F 2 o{Tn}, then F 62 oðT n Tn Þ because T 2 T n Tn and T n 2 T n Tn . So, F 2 oðT n Tn Þ4 fT n g ¼ Cn 4 fT n g ¼ Cnþ1 . – If F 62 o{Tn}, then T is unique in T n Tn , so F 2 oðT n Tn Þ. Hence, F 2 oðT n Tn Þ 4 fT n g ¼ Cn 4 fT n g ¼ Cnþ1 . It comes oðT n Tnþ1 Þ Cnþ1 , so Cnþ1 ¼ oðT n Tnþ1 Þ. So, the advancing front ðCn Þn2N is well deﬁned (because if Cn 5 ;, then oðT n Tn Þ 6¼ ;, so S 9 T n 62 Tn such that Cn \ o{Tn} 5 ;). Now, it can be proved that 80 6 n < N oðX 4 T 0 4 4 T n1 Þ F 2Cn F : S S • With n = 0, since ðN; oTÞ is a mesh of oX, we have oX ¼ F 2oT F ¼ SF 2C0 F . • Now with 0 6 n < N, we suppose S that oðX 4 T 0 4 4 T n1 Þ F 2Cn F . We ﬁrst prove that ðoðX 4 T 0 4 4 T n1 Þ 4 oT n Þ F 2Cnþ1 F : – if x 2 o(X 4 T0 4 4 Tn1) and x 62 oTn, then $ F 2 Cn such that x 2 F and "F 0 2 o{Tn} we have x 62 F 0 . So, F 2 Cn+1. – if x 62 o(X 4 T0 4 4 Tn1) and x 2 oTn, then $ F 2 o{Tn} such that x 2 F and "F 0 2 Cn we have x 62S F. So, F 2 Cn+1. S Since F 2Cnþ1 F is a closed set, we also have oðX 4 T 0 4 4 T n1 Þ 4 oT n F 2Cnþ1 F . Finally, we use S oðA 4 BÞ ¼ oA 4 oB and it comes oðX 4 T 0 4 4 T n Þ F 2Cnþ1 F . S So, oðX 4 T 0 4 4 T N Þ F 2CN F ¼ ; and, since (X 4 T0 4 4 TN) is bounded, we obtain (X 4 T0 4 4 TN) = ;. Finally, since (AnB) (A 4 B), it is easy to prove that (Xn(T0 [ [ TN)) (X 4 T0 4 4 TN). Hence, Xn(T0 [ [ TN) = ;. h With this lemma, we avoid the hypothesis of orientability assumed in [18] when proving the following theorem. Theorem 9. Let T be a mesh topology on a finite node set N in X. Then ðN; TÞ is a mesh of X if and only if the simplices of T are non-degenerated and X jT j ¼ jXj. ð5Þ T 2T

That is, among all mesh topologies, valid meshes are those that satisﬁes the minimal volume criterion (5). P Proof. We just have to prove that if T 2T jT j ¼ jXj then ðN; TÞ is a mesh of X. Since condition (i) is already satisﬁed, we focusSon conditions (ii) and (iii) S of the Deﬁnition S 5. Firstly, suppose that T ¼ 6 X, since X T it comes j T 2T T 2T T 2T T j > jXj. But we always have S P T 2T jT j P j T 2T T j, which raises a contradiction. Hence, condition (ii) is true. P Secondly, Ssuppose that 9 T 1 6¼ T 2 2 T such that T1 \ T2 6¼ ;, then jT1j + jT2j > jT1 [ T2j. It comes T 2T jT j > j T 2T T j, which is the same contradiction. Thirdly, suppose that 9 T 1 6¼ T 2 2 T such that oT1 \ oT2 is not a k-face of T1 and T2 with k 2 {1, 0, . . . , d 1}. Let F1 be a face of T1 and let F2 be a face of T2 such that F1 \ F2 is not a k-face of F1 and F2 with k 2 {1, 0 , . . . , d 2}. If F 1 62 ofTg then exists T 3 2 T n fT 1 ; T 2 g whose F1 is a face. Since T3 is not degenerated, T3 overlap T2, it gives the same contradiction. So F 1 2 oT and similarly F 2 2 oT. But in this case, since F1 \ F2 is not a k-face of F1 and F2 with k 2 {1, 0 , . . . , d 2}, ðN; oTÞ cannot be a mesh of oX. h Finally, the only primitive used to construct a mesh topology is the star operation which consists in connecting one node S to a set of faces F T ðS; FÞ ¼ fconvexhull of fSg [ F where F 2 F and S 62 F g. ð6Þ

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4957

Lemma 10. Let S be a node and F a set of faces. If F is a mesh topology without boundary (i.e. oF ¼ ;), then T ðS; FÞ is a mesh topology and its boundary oðT ðS; FÞÞ ¼ F: This lemma validates the star operation as a mesh topology generator. Proof. Let F0 be a face of T ðS; FÞ that shares three distinct elements: T, T 0 and T 00 respectively convexhull of {S} union F, F 0 and F 00 . Then, we have F \ F0 = F 0 \ F0 = F 00 \ F0 = G. So, G shares three distinct elements of F, which is impossible because F is a mesh topology. Besides, by construction we have F oT ðS; FÞ. Let F be a face of T ðS; FÞ such that F 62 F. F is a face of the convexhull of {S} [ F 0 where F 0 2 F. We denote G = F \ F 0 , since F is without boundary, 9 F 00 6¼ F 0 2 F such that G shares F 0 and F00 . But in this case, F shares the convexhull of {S} [ F 0 and {S} [ F 00 , so F 62 oT ðS; FÞ. h Now, with Theorem 9 and the star operation (6), it is possible to generate a mesh of X. 2.3. Mesh generation The purpose of this section is to construct a mesh of X starting from a mesh of oX (given by some CAD tools for example). The generation technique proposed here does not relie on classical octree, frontal of Delaunay methods. Instead it proceeds by local cut and paste operations. We suppose that X has a connected boundary oX (this condition can be removed with Section 4). The topological mesh generator proceeds thus: Algorithm 1. Starting from ðN; oTÞ a mesh of oX, – build a simple mesh topology T ¼ T ðS 0 ; oTÞ of X where S0 is a node of the boundary (at this step, if X is not convex, ðN; TÞ is not a valid mesh, see Fig. 2 at the top) while the topology T has been changed do for each node or each edge of the mesh topology do – cut the local topology Tc around this node or this edge P – then, try to paste a better local topology Te so as to minimize T 2T jT j among all candidates obtained by the star operator: T ðS 00 ; oTc Þ done done Empirically, this algorithm always leads to the minimal volume jXj (computed with Formula (4)), and then, according to Theorem 9, ðN; TÞ is a valid mesh of X as in Fig. 2 at the end. However, the convergence of this algorithm has not been proved yet. The local topology Tc around a node S is composed by the elements of T whose vertices belong to S ¼ fSg [ NðSÞ where NðSÞ are the neighbors of S (Fig. 3(a)). The candidates for Te are T ðS 00 ; oTc Þ where S00 is S or one vertices of oTc or C the isobarycenter of the vertices of oTc . The local topology Tc around an edge [S, S 0 ] is composed by the elements of T whose vertices belong to S ¼ fS; S 0 g [ NðS; S 0 Þ where NðS; S 0 Þ ¼ NðSÞ \ NðS 0 Þ. In this case, the candidates for Te are T ðS 00 ; oTc Þ where S00 is S, S 0 or one vertices of oTc or C the middle of [S, S 0 ]. The cut and paste operation is licit in any dimension, because of the following result, partially proved in [18]. Proposition 11. With the notations of the Algorithm 1, we denote by Tf the elements of T whose all vertices, except exactly one, belong to S. If:

4958

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

Fig. 2. Mesh generation process by local optimization.

S5 S4

S6

S4 S3

S5

S1

S

S4

S6 S1

S3

S2

S5

S1

S3

S2

(a)

S6

C

(b)

S2 (c)

Fig. 3. The local topology Tc of S and two candidates Te .

(i) oTc is a mesh topology, (ii) the faces of oTc are also the faces of Tf whose all vertices belong to S, then T n Tc [ Te is also a mesh topology of X. That is, after the cut and paste operation, we still have a mesh topology (except when condition (ii) is not fulﬁlled).

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4959

Proof. Let F be a face of T n Tc [ Te that shares strictly more than two simplices. Firstly, since oTc is a mesh topology (condition (i)), T n Tc and Te are mesh topologies. So F is necessary a face of T n Tc and a face of Te . Secondly, let T be an element of T n Tc whose F is a face. Since T 62 Tc , all the vertices of T cannot belong to S. Since F is a face of Te , the vertices of F belong to S. Since F is a face of T, it comes T 2 Tf . Then, condition (ii) says that F belongs to oTc . Thirdly, since F 2 oTc , we have F 2 oT or F 2 oðT n Tc Þ, but in any case F cannot share more than one simplex in T n Tc . Besides, Lemma 10 says that oTe ¼ oTc , so F cannot share more than one simplex in Te , which raises a contradiction. Finally, it remains the proof of oðT n Tc [ Te Þ ¼ oT. • Let F 2 oT, we have 9 !T 2 T such that F 2 o{T}. – If T 62 Tc then F 62 oTc , so 8T 0 2 Te F 62 ofT 0 g. Hence, T is unique in T n Tc [ Te . – If T 2 Tc then F 2 oTc (because T is unique in Tc ) so, 9 !T 0 2 Te such that F 2 o{T 0 }. Since 8T 00 2 T n Tc F 62 ofT 00 g, T 0 is unique in T n Tc [ Te . At this point we have oT oðT n Tc [ Te Þ. • Let F 2 oðT n Tc [ Te Þ, we have 9 !T 2 T n Tc [ Te such that F 2 o{T}. – If T 62 Te then T 2 T n Tc and F 62 oTe . We have 8T 0 2 Tc F 62 fT 0 g, otherwise, T 0 is unique in Tc (because of the presence of T in T) so F 2 oTc , which is impossible because oTe ¼ oTc (Lemma 10). – If T 2 Te then F 2 oTe ¼ oTc (Lemma 10). So, 9 !T 0 2 Tc such that F 2 o{T 0 }. By construction, we have Te \ ðT n Tc Þ ¼ ;. Since T is unique in T n Tc [ Te , 8T 00 2 T n Tc F 62 ofT 00 g. Hence, T 0 is unique in T. It comes oðT n Tc [ Te Þ oT. h With the star operation, edges and faces can be swapped, the node S can be removed (only if S 62 oTc ) and the centroid C can be introduced. But, for a boundary node S 2 oX, the solution to get S 62 oTc is to connect all boundary faces F 2 oX to a ﬁctitious node (say node 0). With these (virtual) elements, the topology is without boundary (i.e. oT ¼ ;, Fig. 4). In this case, the faces of oX are the node 0 opposite faces. The introduction of node 0 allows the treatment of boundary nodes as internal ones. Thus, they do not have a particular treatment, which simpliﬁes the implementation of our dimension independent mesh generator. In particular, boundary nodes can be eliminated, moved or introduced if it does not induce volume changing (the coordinates of node 0 are locally the coordinates of node S).

0

Fig. 4. A topology without boundary and its virtual elements connected to node 0.

4960

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

Finally, this algorithm is more time consuming than a Delaunay mesh generator because it tries many more topologies. However, to improve a mesh by local optimizations turns out to be a good strategy in the following section. 2.4. Anisotropic adaptation Now that the minimal volume has been reached by topological optimization, we usually want to adapt the mesh to an isotropic or anisotropic metric ﬁeld, may it be Euclidian (i.e. constant over X) or not. Let ðN; TÞ be a mesh of X and let M be a continuous P1 metric ﬁeld on this mesh. We only use continuous P1 ﬁelds during mesh generation process, because they are easier to transport (than discontinuous P0 ﬁelds for example). To adapt the mesh to its metric ﬁeld, we proceed by the same local cut and paste technique. But in this case, we have to keep the minimal volume and to improve the quality of the elements. The anisotropic quality criterion [19] we use for an element T whose vertices are S0, . . . , Sd is c0

jT jM T hdM T

ð7Þ

;

where we compute the average metric

the Euclidian volume 2

the L average size

d 1 X MðS i Þ; d þ 1 i¼0 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ jT jM T ¼ jT j detðM T Þ;

MT ¼

hM T ¼

X 2 2 kS j S i kM T dðd þ 1Þ 06i Mx. This anisotropic quality criterion is Euclidian because a constant metric MT is considered over the element T. Note that c0 is such that (7) = 1 when T is equilateral (in the metric MT). Furthermore, in order to introduce or to remove enough nodes before optimizing the shape of the elements, a size criterion is also used d 1 min hM T ; . ð8Þ hM T The ﬁnal criterion c(T) is the worst between the quality criterion (7) and the size criterion (8). At this point, we can introduce an ordering relationship to compare local topologies: let T1 ¼ ðT i1 Þ16i6I 1 and T2 ¼ ðT i2 Þ16i6I 1 (such that 81 6 i1 < i2 6 I j cðT ij1 Þ 6 cðT ij2 Þ) be two local topologies of the same cavity (Fig. 3(b) and (c), for example) 8i < i0 cðT i1 Þ ¼ cðT i2 Þ; T1 < T2 () 9 i0 such that ð9Þ and cðT i10 Þ < cðT i20 Þ. Following this ordering relationship, the next algorithm ensure that the worst element of a maximal topology will be as good as possible. Instead, the use of a cost function [20], only provides an average result that enable degenerated or very poor elements. Moreover, compared to [20,21], which uses the same optimization technique, the only primitive needed here is the star operator.

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4961

Algorithm 2. The algorithm to make the mesh ðN; TÞ ﬁt the metric ﬁeld M is: while the topology T has been changed do for each node and each edge of the mesh do – cut the local topology around this node or this edge – among all minimal volume candidates (obtained by the star operator) paste one maximal topology in the ordering relationship between local topologies – transport the metric ﬁeld, if necessary done done The metric ﬁeld needs to be transported only if a new node C has beenP introduced. In this case, when C is n the isobarycenter of {S1, . . . , Sn}, the vertices of oTc , we take MðCÞ ¼ 1n i¼1 MðS i Þ and when C is the mid0 0 dle of [S, S ], we take M(C) = (M(S) + M(S ))/2. This choice is made for the sake of simplicity, but it remains inaccurate. The mesh we obtain is not directly perfect. By computing again the metric ﬁeld on this mesh and restarting the Algorithm 2, a good mesh can be iteratively produced by this adaptive process (as in Fig. 5 where the metric ﬁeld is deﬁned by [20]).

2.5. Conclusion The introduction of a metric M only changes volume and length calculation. This anisotropic mesh generation method is very robust and its implementation is less technical than for an hybrid Delaunay/frontal [22]. Furthermore, since the algorithm only uses local operations, this mesh generation and adaptation has been easily parallelized [23].

Fig. 5. Anisotropic adaptive process. (a) Initial mesh; ﬁrst step: (b) 1943 nodes in 73 s, (c) 25,066 nodes in 1301 s, and (d) 23,694 nodes in 352 s.

4962

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

The remainder of this paper focuses on the construction of a metric ﬁeld M in the a priori context of mesh generation.

3. Natural metric ﬁeld In this section, we want to address a diﬃcult problem in mesh generation: the meshing of thin and curved geometries (Fig. 6). For that purpose, we propose to take advantage of the anisotropic adaptivity driven by a metric ﬁeld M that places enough elements where it is necessary (i.e. through the thickness), while saving maximum nodes everywhere it is possible (i.e. in other directions), without much human intervention and only from geometrical informations. With several notations, the local anisotropy of a given geometry can be deﬁned (Section 3.1). Then, an approximation of the natural metric can be computed (Sections 3.2–3.4). Finally, the natural metric can be used in placing several element layers through the thickness (Sections 3.5 and 3.6). 3.1. Natural anisotropy deﬁnition Firstly, we need to introduce an ordering relationship between metrics: let M1 and M2 be two metrics ( 81 6 i < i0 k1i ¼ k2i ; ð10Þ M 1 < M 2 () 9 1 6 i0 6 d such that and k1i0 > k2i0 ; where kj1 6 6 kjd are the ordered eigenvalues of Mj. Secondly, the unit ball of a metric M is an ellipsoid EðMÞ ¼ fy 2 Rd such that kykM 6 1g. Thirdly, we denote by o12 X the median surface of X [7].

(a)

(b) Fig. 6. Initial mesh (a) and natural mesh (b) with 4 element layers (193 nodes) in 2D.

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4963

Fig. 7. Natural metric deﬁnition notations.

Finally, we deﬁne the natural metric on x 2 o12 X by the largest metric M such that x þ EðMÞ X (Fig. 7). We deﬁne the natural metric on x 2 X o12 X by the natural metric of x0 2 o12 X such that kx x0k2 is the distance of x to o12 X. If x0 is not unique, we take the averaged natural metric. This deﬁnition cannot be directly implemented because it would involve a inﬁnite number of operations. The proposed method to compute approximately the natural metric is only performed on the nodes S of the mesh. It can be divided in three steps: • detection of the local anisotropy by aggregation of elements around S (Section 3.2), • elliptic interpolation of this element patch (Section 3.3), • and interpolation error estimation to decide when to stop the aggregation (Section 3.4).

3.2. Successive neighborhoods To aggregate a collection of elements around a node S we proceed by successive neighborhoods. Therefore, we deﬁne Nk ðSÞ the kth order neighborhood around S (with k 2 N) by 8 < N0 ðSÞ ¼ fSg; S S ð11Þ Nl ðSÞ [ ðNl ðSÞ \ oXÞ : Nk ðSÞ ¼ NðNk1 ðSÞÞ n 06l kmax the elliptic interpolation error is too large, i.e. when 2 max kS 0 C k kMðSÞ 1 > 0.75. ð18Þ S0 2 N e k ðSÞ f k ðSÞ becomes really non-elliptic (75% is an In other words, we stop the neighborhood growth when N empirically good value for discriminating elliptic and non-elliptic neighborhoods). Then, the appropriate order is K = kmax (or K = kmin when kmin > kmax). Hence, the computed natural metric is M natural ðSÞ ¼ cK LðC K Þ1 with LðC K Þ ¼

Z

ð19Þ

ðy C K Þ ðy C K Þ dy

ð20Þ

CK ðSÞ

and

P cK ¼ P

S0 2 N e K ðSÞ S0 2 N e K ðSÞ

kS 0 C K k2LðCK Þ1 4

kS 0 C K kLðCK Þ1

.

ð21Þ

3.5. Several element layers through the thickness At this point, we know the local anisotropy of the geometry, especially the size and the direction of the thickness. Remember that we want to place enough elements where it is necessary (i.e. through the thickness) while saving maximum nodes everywhere it is possible (i.e. in other directions). 2 We denote by h2 1 P P hd the classiﬁed eigenvalues of Mnatural(S) and by R the rotation matrix such that 3 2 1 0 7 6 h2 7 6 1 7 > 6 .. 7R . M natural ðSÞ ¼ R6 ð22Þ 7 6 . 7 6 4 15 0 h2d

4966

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

In order to get n element layers through the thickness (n being provided by the end user), we have to divide the shortest size h1 by n. But, when the local geometry is isotropic, i.e. when the metric Mnatural(S) is isotropic, all the sizes h1, . . . , hd should be divided by n (in order to generate isotropic elements). In practice, the sizes h1, . . . , hp considered as the thickness are those between h1 and (1 + a)h1 where a > 0 (empirically, we use a = 1). Then, h1, . . . , hp are divided by n, while hp+1, . . . , hd remain unchanged. Finally, we obtain the partially divided natural metric 3 2 2 n 07 6 2 7 6 h1 7 6 7 6 . .. 7 6 7 6 7 6 2 n 7 6 7 6 2 7 > 6 hp 7R . M nnatural ðSÞ ¼ R6 ð23Þ 7 6 1 7 6 7 6 7 6 h2pþ1 7 6 7 6 .. 7 6 7 6 . 7 6 15 4 0 h2d 3.6. Numerical results The previous metric can perfectly treat three-dimensional curvatures (Fig. 10 with n = 6).

Fig. 10. Thin and curved geometry: initial mesh, natural mesh (5441 nodes) and a sectional view.

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4967

Fig. 11. Natural mesh (174,759 nodes) of a complex electrical device (courtesy of Schneider Electric) and a partial view.

Furthermore, the algorithm shows a great robustness, which enables the natural mesh generation of very complex 3D geometries (Fig. 11). In what concerns the computational cost of natural metric construction, the complexity of the algorithm is linear and the computation time on a 220 Mﬂops workstation is about twenty minutes for 105 nodes. Note that the anisotropy of these meshes does not induce any loss of accuracy in the solution given by our ﬁnite element solver [25]. Indeed, it respects the anisotropy of the solution which is often driven by the anisotropy of the geometry. 3.7. Conclusion Our method does not use skeletonization to extract the local thickness [26], and it does not need to extrude a mesh from its boundary or to determine opposite vertices [27,28], so as to get several element layers through the thickness. Instead, we automatically perform a metric ﬁeld that detects the local anisotropy of the geometry. The only a priori information needed to build this natural metric ﬁeld is: an initial mesh ðN; TÞ that deﬁnes the geometry of the part and n, the number of element layers the user wants through the thickness. To compute this natural metric, we gather an element collection around each node, which seems to be a good way to determine the local anisotropy. Then, among all the techniques to determine an elliptic interpolation of this element patch, the length tensor integration and inversion turn out to be particularly robust.

4968

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

Nevertheless, we proceed by successive neighborhoods to collect these element patches, which is certainly not the only possible method. However, we obtain satisfactory results, even on complex parts. 4. Multidomain metric ﬁeld At this point, we are able to generate a natural mesh of a single domain X. Now, let X be composed of several subdomains X = x1 [ [ xp. We want to replace the usual method, which consists in generating one mesh for each subdomain [3,29–31] by the generation of a single mesh which captures the interfaces between subdomains. Our main goal is then to represent internal surfaces in a 3D mesh, up to a desired precision and only from geometrical informations concerning the subdomains xi. Of course inside each subdomain, a natural mesh should be considered (Section 4.1). The ﬁrst idea to represent such internal surfaces, consists in imposing exactly these surfaces in the mesh, either by sub-mesh merging or by well known operations in Delaunay mesh generation technique [22]. However, keeping exactly this internal surfaces in the mesh represents a great constraint on all remeshing operations after. The alternative we develop here relies on a metric ﬁeld that tighten the mesh around internal surfaces (cf. Fig. 12). For that purpose, an interface tensor is added to the natural metric (Section 4.2). This tensor is computed from presence functions which is discussed in Section 4.3. 4.1. Subdomain natural metric The natural mesh of whole X is not very interesting. Instead, the natural metric must be computed on each subdomain: if S 2 xi, then the metric we compute on S is the natural metric of xi. Note that, if a node S belongs to several subdomains, then the natural metric in S does not really matter, because of interface reﬁnement in the following paragraph.

Fig. 12. A 2D multidomain mould mesh (1771 nodes) with a cavity, three prints and two regulation channels.

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4969

To extend the algorithm of the Section 3 onto xi X, it is suﬃcient to ignore the neighbors outside xi when we deﬁne Nk ðSÞ and also to ignore elements whose vertices are not all in xi when oTk ðSÞ is deﬁned. Then, we easily have Mnatural, a metric ﬁeld on X that is natural for each subdomain. Again, we can divide this natural metric on xi by ni, the number of element layers through xi. Note that, in the notation M nnatural , n may be diﬀerent on each subdomain. Hence, the mesh of X will be adapted diﬀerently on each subdomain x1, . . . , xp. But it is not enough to capture the interfaces between them. An additional term is needed and described in the following paragraph. 4.2. Interface tensor deﬁnition Let x X be a subdomain with non-empty interior. We denote by d(x, ox) the signed distance between x 2 Rd and ox (that is, d(x, ox) 6 0 when x 2 x and d(x, ox) > 0 when x 62 x). A strategy to tighten the mesh around the border of x X can be • to take a function g such that g = 1 inside x, g = 0 outside and 0 < g < 1 in the neighborhood of ox. For example 1 gex ðxÞ ¼ ð24Þ 1 þ eadðx;oxÞ with a 1 (Fig. 13) • and to triangulate its graph G ¼ x 2 Rdþ1 such that ðx1 ; . . . ; xd Þ 2 X and xdþ1 ¼ gðx1 ; . . . ; xd Þ ;

ð25Þ

so as to get a mesh that captures ox by orthogonal projection on X, • the triangulation of G being driven in Rdþ1 by the metric M nnatural ðx1 ; . . . ; xd Þ Mðx1 ; . . . ; xdþ1 Þ ¼ ; m2

ð26Þ

where m is the number of element layers the user wants around the interface ox.

1 0.8 0.6 0.4 0.2 0 2 2

1 1

0

0

–1

–1 –2 –2

Fig. 13. Function gex with dðx; oxÞ ¼ kx ck2 r.

4970

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

In fact, this strategy is not directly applied. Proposition 12. The induced metric by M on G (and thus on X by orthogonal projection) is M multidomain ðx1 ; . . . ; xd Þ ¼ M nnatural ðx1 ; . . . ; xd Þ þ m2 rgðx1 ; . . . ; xd Þ rgðx1 ; . . . ; xd Þ.

ð27Þ

That is, the previous strategy in Rdþ1 is equivalent to the use of metric (27) in X Rd . Proof. We denote by f the orthogonal projection on Rd : f(x1, . . . , xd+1) = (x1, . . . , xd). According to Riemannian geometry [15,16], the induced metric by M on X (with map f ) have components 1 1 M ijmultidomain ¼ ðofoxj Þ> M ofoxi . Since on G we have f1(x1, . . . , xd) = (x1, . . . , xd,g(x1, . . . , xd)) we get

of 1 oxi

og > ¼ ðd1i ; . . . ; ddi ; ox Þ . Hence, we obtain i > ! " og og og og ¼ d1j ; . . . ; ddj M nnatural ðd1i ; . . . ; ddi Þ> þ m2 M ijmultidomain ¼ d1j ; . . . ; ddj ; M d1i ; . . . ; ddi ; oxj oxi oxj oxi ! n "ij ij 2 ¼ M natural þ m ðrg rgÞ ;

which gives exactly the metric (27). h Note that the interface tensor $g $g is not a metric, because 0 is always an eigenvalue of x x with x 2 Rd when d P 2. But in addition with a metric, the result is also a metric. At this stage, the multidomain metric involves the computation of the costly function gex. In the next paragraph, gex is substituted by a cheapest function g such that g = 1 inside x, g = 0 outside and 0 < g < 1 in the neighborhood of ox. 4.3. The presence function and its gradient A subdomain x X can be represented by its presence function x ðxÞ

x

: X ! f0; 1g such that

¼ 1 () x 2 x.

ð28Þ 0

In practice, to compute eﬃciently the interface tensor, the function g taken into account is the P discontinuous interpolation of x . So, g: X ! [0, 1] is deﬁned (almost everywhere) by 8T 2 T

8x 2 T

gðxÞ ¼ gðT Þ ¼

jT \ xj . jT j

Proposition 13. The interpolation errors are !1=2 X kg x kL2 ðXÞ ¼ gðT Þð1 gðT ÞÞjT j ;

ð29Þ

ð30Þ

T 2T

kg

x kL1 ðXÞ

¼ 0.

So, considering g instead of small.

ð31Þ x

is conservative and precise when all elements T such that 0 < g(T) < 1 are

R R P Proof. We have kg x kL2 ðXÞ ¼ ð X ðg x Þ2 dxÞ1=2 ¼ ð T 2T T ðg x Þ2 dxÞ1=2 . Since g is constant on T R R R R we get T ðg x Þ2 dx ¼ gðT Þ2 jT j 2gðT Þ T x dx þ T 2x dx. As T x dx ¼ jT \ xj and 2x ¼ x , we R 2 2 obtain T ðg x Þ dx ¼ gðT Þ jT j 2gðT ÞjT \ xj þ jT \ xj. And using jT \ xj = g(T)jTj, it comes R 2 which proves Formula (30). Besides, we have T ðg x Þ dx ¼ gðT ÞjT jðgðT Þ 2gðT Þ þ 1Þ

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4971

T F n(F)

T(F)

T

ω

(a)

(b) 0

Fig. 14. Pixelization and P gradient notations.

kg

x kL1 ðXÞ

mula (31).

¼

R

X ðg

x Þ dx

¼

P

R

T 2T T ðg

x Þ dx

¼

P

jT \xj T 2T ð jT j jT j

jT \ xjÞ ¼ 0 which proves For-

h

To calculate g we pixelize or voxelize X up to a precision (the choice of this precision gives the thickness of the reﬁnement). Then, we light on the pixels or the voxels that belong to x (Fig. 14(a)) and we take 8T 2 T

jT \ xj lit voxels in T ¼ . jT j all voxels in T

ð32Þ

Hence, the only information we need about x is a function that says if a point is inside or not. The numerical gradient of g is a P0 discontinuous vector ﬁeld rg: X ! Rd deﬁned (almost everywhere) by 8T 2 T

8x 2 T rgðxÞ ¼ rgðT Þ ¼

X gðT ðF ÞÞ gðT Þ jF jnðF Þ jT ðF Þj þ jT j F 2oT

ð33Þ

with oT the set of T faces, T(F) the opposite element of T through F, and n(F) the outgoing normal of T on F (Fig. 14(b)). This P0 gradient needs boundary conditions. When F 2 oT, we take jT(F)j = 1, which corresponds to adiabatic boundary conditions. Since there are several presence functions g1, . . . , gp corresponding to diﬀerent subdomains, for each component 1 6 i 6 d we take the maximum

Fig. 15. The mesh that deﬁnes x, its voxelization and the initial mesh of X.

4972

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

! " rgðT Þi ¼ max rg1 ðT Þi ; . . . ; rgp ðT Þi ;

ð34Þ

which leads to a good treatment of triple points. Finally, we obtain a continuous P1 multidomain metric ﬁeld, deﬁned on each node S of the mesh by M multidomain ðSÞ ¼ M nnatural ðSÞ þ m2 rgðSÞ rgðSÞ with rgðSÞ ¼

X T 2TðSÞ

!, jT jrgðT Þ

X

ð35Þ

! jT j ;

T 2TðSÞ

where TðSÞ is the elements of T whose S is a vertex.

Fig. 16. Iterations between the metric generator and the mesh generator.

ð36Þ

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

Fig. 17. The ﬁnal mesh of X with internal surface ox captured (150,683 nodes), and a sectional view.

Fig. 18. A partial view of the previous sectional view.

4973

4974

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4.4. Numerical results We consider half of a symmetrical part x (Fig. 15(a)) in a rectangular domain X (Fig. 15(c)). This part is a biomedical application (courtesy of Sophysa), which is a benchmark for the multidomain ﬁnite element solver [1]. The voxelization of x (Fig. 15(b)) is accurate enough to keep all geometrical informations and is used to calculate eﬃciently the presence function of x. Starting from a very coarse mesh (Fig. 15(c)), we generate the multidomain metric and we adapt the mesh to this metric. This process is iterated several times (Fig. 16).

Fig. 19. The extracted meshes where g = 1, 0 < g < 1 and g = 0.

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4975

At iteration 6, the boundary seems to be good (Fig. 17(a)) and we can see that the mesh is well adapted everywhere in X, not only on oX (Fig. 17(b)). The reﬁnement around the interface between x and X nx is really anisotropic (Fig. 18). Finally, we can extract the sub-meshes of x and Xnx and see the natural metric eﬀect (Fig. 19(a) and (c)). The elements T where 0 < g(T) < 1 are very thin (Fig. 19(b)). Hence, the interpolation error kg x kL2 ðXÞ is very low. 4.5. Conclusion To enrich a mesh around the interfaces between several subdomains, we use presence functions g to locate interfaces and their gradients $g to compute the interface tensor $g $g. Furthermore, this tensor gives a multidomain metric that works properly when it is used in addition to the natural metric we deﬁned previously. The only a priori information we need to build the multidomain metric ﬁeld is: for each subdomain, its natural metric and a function that returns if a point is inside or not. Like [32,33] we progressively detect and reﬁne the elements around interfaces. But, we can deal with complex geometries and the reﬁnement is anisotropic. Although this technique is still little node consuming, we obtain good results.

5. Conclusion The dimension independent mesh generator we introduce in this paper shows a great robustness in adapting an unstructured mesh to any anisotropic metric ﬁeld. Besides, we deﬁne a multidomain metric ﬁeld that captures internal surfaces between subdomains. This metric is computed only from geometry information and uses a natural metric ﬁeld inside each subdomain. Since this natural metric can treat thinness, curvature and changing thickness, we are now able to treat complex parts, whose meshes cannot be generated manually. The ﬁnite element solvers we use with the meshes discussed below, do not show obvious sensitivity to nearly degenerated elements (when they are stretched along the ﬂow direction). The metric ﬁeld method we propose here may be used with other anisotropic mesh generators. Although these algorithms were developed for material processing simulations, they are applicable to other numerical ﬁelds.

Acknowledgement The authors thank gratefully to the referee for helpful suggestions, and to the Rem3D Consortium partners for their support.

References [1] S. Batkam, T. Coupez, A new multi-domain approach for the thermal problem solution in 3D mold ﬁlling, in: ESAFORM Conference on Material Forming, Liege, 2001, pp. 31–34. [2] J. Bruchon, T. Coupez, A study of the 3D polymer foam formation based on the simulation of anisothermal bubble expansion, Mecanique Indust. 4 (2003) 331–338. [3] J. Barboza, L. Fourment, J.-L. Chenot, Contact algorithm for 3D multi-bodies problems: Application to forming of multimaterial parts and tool deﬂection, in: 5th World Congress on Computational Mechanics, Vienna, 2002.

4976

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

[4] P.O. Bouchard, F. Bay, Y. Chastel, I. Tovena, Crack propagation modelling using an advanced remeshing technique, Comput. Meth. Appl. Mech. Engrg. 189 (2000) 723–742. [5] C. Gruau, T. Coupez, Anisotropic and multidomain mesh automatic generation for viscous ﬂow ﬁnite element method, in: N.-E. Wiberg, P. Dı`ez, Adaptive Modelling and Simulation, Go¨teborg, 2003. ISBN: 84-95999-30-7. [6] M. Bellet, C. Aliaga, O. Jaouen, Finite elements for a thermomechanical of analysis of solidiﬁcation processes, in: 9th International Conference on Modeling of Casting, Welding and Advanced Solidiﬁcation Processes, Aachen, 2000, pp. 10–17. [7] P.J. Frey, P.-L. George, Mesh generation—Application to ﬁnite elements, Hermes, 2000. [8] J. Dompierre, M.-G. Vallet, Y. Bourgault, M. Fortin, W.G. Habashi, Anisotropic mesh adaptation: Towards user-independent, mesh-independent and solver-Independent CFD. Part III: Unstructured meshes, Int. J. Numer. Meth. Fluids 39 (2002) 675–702. [9] J. Peraire, K. Morgan, Unstructured mesh generation including directional reﬁnement for aerodynamic ﬂow simulation, Finite Elements Anal. Des. 25 (1997) 343–355. [10] R. Lo¨hner, C. Yang, E. On˜ate, S. Idelssohn, An unstructured grid-based, parallel free surface solver, Appl. Numer. Math. 31 (1999) 271–293. [11] F. Alauzet, P.J. Frey, B. Mohammadi, Anisotropic mesh adaptation. Application to transient CFD problem, in: Second MIT Conference on Computational Fluid and Solid Mechanics, Boston, 2003, pp. 1851–1854. [12] M.J. Castro-Dı´az, F. Hecht, B. Mohammadi, O. Pironneau, Anisotropic unstructured mesh adaptation for ﬂow simulations, Int. J. Numer. Meth. Fluids 25 (1997) 475–491. [13] T. Coupez, Grandes transformations et remaillage automatique, PhD Thesis, France, 1991, pp. 123–210. [14] P.D. Zavattieri, E.A. Dari, G.C. Buscaglia, Optimization strategies in unstructured mesh generation, Int. J. Numer. Meth. Eng. 39 (1996) 2055–2071. [15] M.P. doCarmo, Riemannian Geometry, Birkha¨user, 1988. [16] J. Jost, Riemannian Geometry and Geometric Analysis, Springer-Verlag, 1998. [17] P.M. Morse, H. Feshbach, Stokes theorem, in: Methods of Theoretical Physics, Part I, New York, 1953, p. 43. [18] T. Coupez, Generation de maillage et adaptation de maillage par optimisation locale, Revue europeenne des elements ﬁnis 9 (2000) 403–423. [19] J. Dompierre, M.-G. Vallet, P. Labbe´, F. Guibault, On simplex shape measures with extension for anisotropic meshes, in: Workshop on Mesh Quality and Dynamic Meshing, Livermore, 2003, pp. 46–71. [20] C. Bottasso, On the choice of the optimization strategy and of the objective function for metric-driven mesh adaptation, in: International Conference on Adaptative Modelling and Simulation, Go¨teborg, 2003. [21] C.C. Pain, A.P. Umpleby, C.R.E. de Oliveira, A.J.H. Goddard, Tetrahedral mesh optimisation and adaptativity for steady-state and transient ﬁnite element calculations, Comput. Methods Appl. Mech. Engrg. 190 (2001) 3771–3796. [22] P.-L. George, H. Borouchaki, Delaunay triangulation and meshing. Application to ﬁnite elements, Hermes, 1998. [23] T. Coupez, H. Digonnet, R. Ducloux, Parallel meshing and remeshing, Appl. Math. Modell. 25 (2000) 153–175. [24] S.G. Advani, C.L. Tucker, The use of tensors to describe and predict ﬁber orientation in short ﬁber composites, J. Rheol. 31 (1987) 751–784. [25] T. Coupez, E. Bigot, 3D anisotropic mesh generation and adaptation with applications, in: European Congress on Computational Methods in Applied Sciences and Engineering, Barcelona, 2000. [26] K.-F. Tchon, M. Khachan, F. Guibault, R. Camarero, Constructing anisotropic geometric metrics using octrees and skeletons, in: 12th International Meshing Roundtable, Santa Fe, 2003, pp. 293–304. [27] R.V. Garimella, M.S. Shephard, Boundary layer meshing for viscous ﬂows in complex domains, in: 7th International Meshing Roundtable, Dearborn, 1998, pp. 107–118. [28] R.V. Garimella, Anisotropic tetrahedral mesh generation, PhD thesis, Rensselaer Polytechnic Institute, 1998. [29] S. Chakravarthy, A uniﬁed-grid ﬁnite volume formulation for computational ﬂuid dynamics, Int. J. Numer. Meth. Fluids 31 (1999) 309–323. [30] T. Tezduyar, Y. Osawa, The multi-domain method for computation of the aerodynamics of a parachute crossing the far wake of an aircraft, Comput. Methods Appl. Mech. Engrg. 191 (2001) 705–716. [31] M. Lesoinne, C. Farhat, Geometric conservation law for ﬂow problems with moving boundaries and deformable meshes, and their impact on aeroelastic computations, Comput. Methods Appl. Mech. Engrg. 134 (1996) 71–90. [32] T. Stalicky´, H.-G. Roos, Anisotropic mesh reﬁnement for problems with internal and boundary layers, Int. J. Numer. Meth. Engrg. 46 (1999) 1933–1953. [33] T. Tezduyar, S. Aliabadi, M. Behr, Enhanced-discretisation interface-capturing technique (EDICT) for computation of unsteady ﬂows with interfaces, Comput. Methods Appl. Mech. Engrg. 155 (1998) 235–248.

3D tetrahedral, unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric Cyril Gruau *, Thierry Coupez Material Processing Center (CEMEF), E´cole des Mines de Paris, BP 207, 06904 Sophia-Antipolis, France Received 28 November 2003; received in revised form 25 June 2004; accepted 26 November 2004

Abstract In this paper we present a 3D tetrahedral, unstructured and anisotropic mesh generator that is not based on the Delaunay, frontal or octree method. Instead, it proceeds by local optimizations and uses an anisotropic shape criterion to ﬁt a metric ﬁeld. Then, we introduce a new 3D metric ﬁeld that tightens the mesh around interfaces when the calculation domain is divided in several subdomains, and a 3D metric ﬁeld that places enough elements through each subdomain thickness, without introducing too many nodes in the other directions. Finally, we show some applications for material forming geometries. 2005 Elsevier B.V. All rights reserved. Keywords: Topology and shape optimization; Elliptic interpolation; Thickness detection and curvature treatment; Interface reﬁnement

1. Introduction Recent advances in numerical methods together with high performance computing made it possible to simulate 3D applications involving several objects (i.e. several subdomains). In material processing, these multidomain simulations may concern injection moulding process [1], foam expansion [2], forging process [3] or composite forming [4]. In this context, numerical software has to treat increasingly complex three-dimensional geometries, especially with industrial thermoplastic parts [5] and with the casting process [6]. Today, the great challenge is to

*

Corresponding author. E-mail addresses: [email protected] (C. Gruau), [email protected] (T. Coupez).

0045-7825/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.11.020

4952

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

generate suitable meshes for computation, without much human intervention. In our case, a suitable mesh means: • a mesh with several element layers everywhere through the thickness (to observe non-Newtonian behavior and thermal eﬀects, for example), even for thin and curved geometries, • a mesh that captures precisely the interfaces between subdomains (between ﬂuid and structure for example), • a mesh that does not contains too many nodes (because of memory capacity and computational time limitations). Therefore, we tend to generate anisotropic meshes, because isotropic meshes would contain too many nodes for an implicit three-dimensional ﬁnite element solver. So, we need a mesh generator that can take a metric ﬁeld into account. Furthermore, we have to deﬁne a metric ﬁeld that meets the common needs detailed above. The originality of this paper is to propose an alternative to the existing 3D, unstructured and anisotropic mesh generators [7–10], and to present two kinds of 3D metric ﬁelds that are neither produced by a posteriori information [11], nor designed for boundary layer improvement [12], i.e.: • the natural metric, that places enough element layers through the thickness while saving maximum nodes in the other directions (see a natural mesh, Fig. 1(a)), • and the multidomain metric, that gathers various subdomains in a single mesh (see a multidomain mesh, Fig. 1(b) and (c)).

Fig. 1. A natural mesh of an optical lens quarter (2644 nodes), a bidomain mesh of this lens in its mould (22,893 nodes) and a sectional view.

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4953

This paper is organized as follows: In Section 2, we introduce notations and deﬁnitions to explain our anisotropic mesh generation process. In Section 3, we deﬁne the natural metric, we indicate a way to compute it and we show some results. In Section 4, we do the same with the multidomain metric. Finally, concluding remarks are gathered in Section 5.

2. Mesh generation and anisotropic adaptation by topological optimization In this section, an adaptive method for generating anisotropic and unstructured meshes in any dimension is described. Firstly, useful properties about meshes and mesh topologies are established (Sections 2.1 and 2.2). Secondly, Section 2.3 details a mesh generation by topological optimization, which was ﬁrst introduced in [13] and reused by [14]. Thirdly, we expose how this adaptive method can be used to conform a mesh to a nonEuclidian metric ﬁeld (Section 2.4). From now, our calculation domain is a polytope X 2 Rd (with d = 2, 3, . . .), which is bounded, orientable and with non-empty interior. 2.1. Basic deﬁnitions and properties For the sake of completeness, the purpose of this paragraph is to deﬁne the classical notions that are widely used in this paper. Deﬁnition 1. A metric is a d-rowed and d-columned, real, symmetrical, positive deﬁnite matrix. In fact, according to the following proposition, a metric M consists in d orthogonal directions (the columns of R) and one size hi in each direction. Proposition 2. A metric M can be diagonalized in 3 2 1 0 7 6 h2 7 6 1 7 > 6 .. 7R ; M ¼ R6 7 6 . 7 6 4 15 0 h2d

ð1Þ

where R is a rotation matrix whose columns are eigenvectors of M and ðh2 i Þ16i6d the corresponding eigenvalues. Proof. Since M is symmetrical, M can be diagonalized in the orthogonal group (at this stage, R = R>). Since M is positive deﬁnite, the eigenvalues k of M are strictly positive (so, $ h 5 0 such that k = 1/h2). For the columns of R, we can choose unit eigenvectors of M (then, jdet(R)j = 1). If det(R) = 1, we can choose the opposite of one eigenvector so as to get det(R) = 1 (hence, R is a rotation matrix). h Before deﬁning the requirement for a mesh, we need the deﬁnition of a simplex in any dimension and with any number of vertices. Deﬁnition 3. A k-simplex of Rd (with 0 6 k 6 d) is the convexhull of k + 1 points in Rd (called vertices) and a k-face of a simplex T is a k-simplex whose vertices are vertices of T. It is usual to say that a (k 1)-face of a k-simplex is shortly a face and a (1)-face is the empty set. Now, in many steps the calculus of the volume of a simplex is required.

4954

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

Proposition 4. Let T be k-simplex of Rd with vertices S0, . . . , Sk. We denote by A the rectangular matrix (S0S1, . . . , S0Sk). The volume of T is pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ detðA> AÞ jT j ¼ . ð2Þ k! So, we can compute the volume of a simplex for each dimension d and each simplex type k. Formula (2) can be rewritten as • when k = 1

jT j ¼ kS 0 S 1 k, j detðAÞj , • when k = d jT j ¼ d! kS 0 S 1 ^ ^ S 0 S d1 k • when k = d 1 jT j ¼ , ðd 1Þ! where ^ ^ stands for the cross product in Rd . Proof. Let Tb be the k-simplex of Rk whose vertices are (dij)16i,j6k, one can easily prove that j Tb j ¼ k!1 . We denote by f the linear function such that f ð Tb Þ ¼ T and $f = A. From Riemannian geometry [15,16], we know that Rd (provided with its canonical Euclidian metric) induces on T (with its map f) the metric pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ R pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ R pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ M = A>A and that jT j ¼ T^ detðMÞ dx. Since M is constant we have T^ detðMÞ dx ¼ j Tb j detðMÞ. Hence, we obtain Formula (2). h A simplex T is said to be non-degenerated if and only if jTj > 0. Now, we are able to deﬁne a mesh. Deﬁnition 5. Let N be a ﬁnite node set in X, let T be a set of non-degenerated d-simplices whose vertices belong to N and let F be the faces of those simplices. ðN; TÞ is a mesh of X if and only if: (i) S each face of F does not share more than two simplices of T, (ii) T 2T T ¼ X, (iii) for all simplices T 1 6¼ T 2 2 T, the intersection T1 \ T2 is a k-face of both simplices with k 2 {1, 0 , . . . , d 1}. The (topological) boundary of a mesh ðN; TÞ is oT, the faces of F that belong to only one simplex in T. We have the same deﬁnition for a mesh ðN; UÞ of oX where U is a set of (d 1)-simplices (which is used in Deﬁnition 7). It is obvious that if ðN; TÞ is a mesh of X, then ðN; oTÞ may not be a mesh of oX, because condition (i) may not be satisﬁed. However, the volume of X can be computed with oT alone. Proposition 6. Let ðN; TÞ be a mesh of X. The volume of this mesh is X jXj ¼ jT j;

ð3Þ

T 2T

Z 1 ¼ d S

x ~ n dx F 2oT

ð4Þ

F

with ~ n the unit normal vector pointing outward. Owing to Formula (4) the volume of a mesh can be computed using only its boundary. R R Proof. We have jXj ¼ X dx ¼ S T dx. Since simplices of T do not overlap themselves, we have T 2T RS R P P dx ¼ dx ¼ jT j. T 2T

T

T 2T T

T 2T

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4955

Besides, according theorem [17], as the divergence $ Æ x =S d and since X is an orientable R to Stokes R R S polytope, we have X dx ¼ d1 X r x dx ¼ d1 oX x ~ n dx. Finally, since X ¼ T 2T T , we have oX ¼ F 2oT F , which leads to Formula (4). h 2.2. The mesh topology deﬁnition and the minimal volume criterion Since a mesh is a very diﬃcult object to build directly, especially when d P 3, we introduce an easier object (a mesh topology) in the following deﬁnition. The purpose of this section is to rise a suﬃcient condition to obtain a mesh from a mesh topology. Deﬁnition 7. Let N be a ﬁnite set of nodes in X, let T be a set of d-simplices whose vertices belong to N and let F be the faces of those simplices. T is a mesh topology if and only if: (i) each face of F does not share more than two simplices of T, (ii) ðN; oTÞ is a mesh of oX. Again, in this deﬁnition, the boundary oT of a mesh topology T is the faces of F that belong to only one simplex in T. We have the same deﬁnition for a mesh topology ðN; UÞ of oX where U is a set of (d 1)-simplices (which is used in Lemma 10). R Thanks to condition (ii), we still have jXj ¼ d1 S F x ~ n dx. However, all mesh topologies are not F 2oT

meshes because their simplices can be degenerated, overlap themselves or overﬂow X. Hence, with this P definition we do not have jXj ¼ T 2T jT j but only the following result. Lemma 8. Let T be a mesh topology, then X

S

T 2T T .

That is, thanks to condition (ii) of Deﬁnition 7, each point of the calculation domain X lies inside an element of T. Proof. We consider the following advancing front: • C0 ¼ oT and T0 ¼ ; • 8n 2 N – if Cn = ; then Cn+1 = Cn – if Cn 5 ; then Cn+1 = Cn 4 o{Tn} (with the operator A 4 B = (A [ B)n(A \ B)) and Tnþ1 ¼ Tn [ fT n g where T n 62 Tn is such that Cn \ o{Tn} 5 ;. Since T is a ﬁnite set, 9 N 2 N such that CN1 5 ; and "n P N Cn = ;. Now, it can be proved that 80 6 n < N Cn ¼ oðT n Tn Þ: • With n = 0, we have C0 ¼ oT ¼ oðT n T0 Þ because T0 ¼ ;. • Now with 0 6 n < N, if Cn ¼ oðT n Tn Þ then let F 2 Cn+1. – If F 2 Cn, then F 62 o{Tn} (because Cn+1 = Cn 4 o{Tn}) so, 9 !T 6¼ T n 2 T n Tn such that F 2 o{T}. Hence, F 2 oðT n Tn n fT n gÞ ¼ oðT n Tnþ1 Þ. – If F 62 Cn, then F 2 o{Tn} and F 62 oT so, condition (ii) of the Deﬁnition 7 gives 9 !T 6¼ T n 2 T such that F 2 o{T}. We have T 62 Tn , otherwise F 2 oðT n Tn Þ ¼ Cn . Hence, F 2 oðT n Tn n fT n gÞ ¼ oðT n Tnþ1 Þ. So, Cnþ1 oðT n Tnþ1 Þ. Now, let F 2 oðT n Tnþ1 Þ ¼ oðT n Tn n fT n gÞ. We have 9 !T 2 T n Tn n fT n g such that F 2 o{T}.

4956

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

– If F 2 o{Tn}, then F 62 oðT n Tn Þ because T 2 T n Tn and T n 2 T n Tn . So, F 2 oðT n Tn Þ4 fT n g ¼ Cn 4 fT n g ¼ Cnþ1 . – If F 62 o{Tn}, then T is unique in T n Tn , so F 2 oðT n Tn Þ. Hence, F 2 oðT n Tn Þ 4 fT n g ¼ Cn 4 fT n g ¼ Cnþ1 . It comes oðT n Tnþ1 Þ Cnþ1 , so Cnþ1 ¼ oðT n Tnþ1 Þ. So, the advancing front ðCn Þn2N is well deﬁned (because if Cn 5 ;, then oðT n Tn Þ 6¼ ;, so S 9 T n 62 Tn such that Cn \ o{Tn} 5 ;). Now, it can be proved that 80 6 n < N oðX 4 T 0 4 4 T n1 Þ F 2Cn F : S S • With n = 0, since ðN; oTÞ is a mesh of oX, we have oX ¼ F 2oT F ¼ SF 2C0 F . • Now with 0 6 n < N, we suppose S that oðX 4 T 0 4 4 T n1 Þ F 2Cn F . We ﬁrst prove that ðoðX 4 T 0 4 4 T n1 Þ 4 oT n Þ F 2Cnþ1 F : – if x 2 o(X 4 T0 4 4 Tn1) and x 62 oTn, then $ F 2 Cn such that x 2 F and "F 0 2 o{Tn} we have x 62 F 0 . So, F 2 Cn+1. – if x 62 o(X 4 T0 4 4 Tn1) and x 2 oTn, then $ F 2 o{Tn} such that x 2 F and "F 0 2 Cn we have x 62S F. So, F 2 Cn+1. S Since F 2Cnþ1 F is a closed set, we also have oðX 4 T 0 4 4 T n1 Þ 4 oT n F 2Cnþ1 F . Finally, we use S oðA 4 BÞ ¼ oA 4 oB and it comes oðX 4 T 0 4 4 T n Þ F 2Cnþ1 F . S So, oðX 4 T 0 4 4 T N Þ F 2CN F ¼ ; and, since (X 4 T0 4 4 TN) is bounded, we obtain (X 4 T0 4 4 TN) = ;. Finally, since (AnB) (A 4 B), it is easy to prove that (Xn(T0 [ [ TN)) (X 4 T0 4 4 TN). Hence, Xn(T0 [ [ TN) = ;. h With this lemma, we avoid the hypothesis of orientability assumed in [18] when proving the following theorem. Theorem 9. Let T be a mesh topology on a finite node set N in X. Then ðN; TÞ is a mesh of X if and only if the simplices of T are non-degenerated and X jT j ¼ jXj. ð5Þ T 2T

That is, among all mesh topologies, valid meshes are those that satisﬁes the minimal volume criterion (5). P Proof. We just have to prove that if T 2T jT j ¼ jXj then ðN; TÞ is a mesh of X. Since condition (i) is already satisﬁed, we focusSon conditions (ii) and (iii) S of the Deﬁnition S 5. Firstly, suppose that T ¼ 6 X, since X T it comes j T 2T T 2T T 2T T j > jXj. But we always have S P T 2T jT j P j T 2T T j, which raises a contradiction. Hence, condition (ii) is true. P Secondly, Ssuppose that 9 T 1 6¼ T 2 2 T such that T1 \ T2 6¼ ;, then jT1j + jT2j > jT1 [ T2j. It comes T 2T jT j > j T 2T T j, which is the same contradiction. Thirdly, suppose that 9 T 1 6¼ T 2 2 T such that oT1 \ oT2 is not a k-face of T1 and T2 with k 2 {1, 0, . . . , d 1}. Let F1 be a face of T1 and let F2 be a face of T2 such that F1 \ F2 is not a k-face of F1 and F2 with k 2 {1, 0 , . . . , d 2}. If F 1 62 ofTg then exists T 3 2 T n fT 1 ; T 2 g whose F1 is a face. Since T3 is not degenerated, T3 overlap T2, it gives the same contradiction. So F 1 2 oT and similarly F 2 2 oT. But in this case, since F1 \ F2 is not a k-face of F1 and F2 with k 2 {1, 0 , . . . , d 2}, ðN; oTÞ cannot be a mesh of oX. h Finally, the only primitive used to construct a mesh topology is the star operation which consists in connecting one node S to a set of faces F T ðS; FÞ ¼ fconvexhull of fSg [ F where F 2 F and S 62 F g. ð6Þ

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4957

Lemma 10. Let S be a node and F a set of faces. If F is a mesh topology without boundary (i.e. oF ¼ ;), then T ðS; FÞ is a mesh topology and its boundary oðT ðS; FÞÞ ¼ F: This lemma validates the star operation as a mesh topology generator. Proof. Let F0 be a face of T ðS; FÞ that shares three distinct elements: T, T 0 and T 00 respectively convexhull of {S} union F, F 0 and F 00 . Then, we have F \ F0 = F 0 \ F0 = F 00 \ F0 = G. So, G shares three distinct elements of F, which is impossible because F is a mesh topology. Besides, by construction we have F oT ðS; FÞ. Let F be a face of T ðS; FÞ such that F 62 F. F is a face of the convexhull of {S} [ F 0 where F 0 2 F. We denote G = F \ F 0 , since F is without boundary, 9 F 00 6¼ F 0 2 F such that G shares F 0 and F00 . But in this case, F shares the convexhull of {S} [ F 0 and {S} [ F 00 , so F 62 oT ðS; FÞ. h Now, with Theorem 9 and the star operation (6), it is possible to generate a mesh of X. 2.3. Mesh generation The purpose of this section is to construct a mesh of X starting from a mesh of oX (given by some CAD tools for example). The generation technique proposed here does not relie on classical octree, frontal of Delaunay methods. Instead it proceeds by local cut and paste operations. We suppose that X has a connected boundary oX (this condition can be removed with Section 4). The topological mesh generator proceeds thus: Algorithm 1. Starting from ðN; oTÞ a mesh of oX, – build a simple mesh topology T ¼ T ðS 0 ; oTÞ of X where S0 is a node of the boundary (at this step, if X is not convex, ðN; TÞ is not a valid mesh, see Fig. 2 at the top) while the topology T has been changed do for each node or each edge of the mesh topology do – cut the local topology Tc around this node or this edge P – then, try to paste a better local topology Te so as to minimize T 2T jT j among all candidates obtained by the star operator: T ðS 00 ; oTc Þ done done Empirically, this algorithm always leads to the minimal volume jXj (computed with Formula (4)), and then, according to Theorem 9, ðN; TÞ is a valid mesh of X as in Fig. 2 at the end. However, the convergence of this algorithm has not been proved yet. The local topology Tc around a node S is composed by the elements of T whose vertices belong to S ¼ fSg [ NðSÞ where NðSÞ are the neighbors of S (Fig. 3(a)). The candidates for Te are T ðS 00 ; oTc Þ where S00 is S or one vertices of oTc or C the isobarycenter of the vertices of oTc . The local topology Tc around an edge [S, S 0 ] is composed by the elements of T whose vertices belong to S ¼ fS; S 0 g [ NðS; S 0 Þ where NðS; S 0 Þ ¼ NðSÞ \ NðS 0 Þ. In this case, the candidates for Te are T ðS 00 ; oTc Þ where S00 is S, S 0 or one vertices of oTc or C the middle of [S, S 0 ]. The cut and paste operation is licit in any dimension, because of the following result, partially proved in [18]. Proposition 11. With the notations of the Algorithm 1, we denote by Tf the elements of T whose all vertices, except exactly one, belong to S. If:

4958

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

Fig. 2. Mesh generation process by local optimization.

S5 S4

S6

S4 S3

S5

S1

S

S4

S6 S1

S3

S2

S5

S1

S3

S2

(a)

S6

C

(b)

S2 (c)

Fig. 3. The local topology Tc of S and two candidates Te .

(i) oTc is a mesh topology, (ii) the faces of oTc are also the faces of Tf whose all vertices belong to S, then T n Tc [ Te is also a mesh topology of X. That is, after the cut and paste operation, we still have a mesh topology (except when condition (ii) is not fulﬁlled).

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4959

Proof. Let F be a face of T n Tc [ Te that shares strictly more than two simplices. Firstly, since oTc is a mesh topology (condition (i)), T n Tc and Te are mesh topologies. So F is necessary a face of T n Tc and a face of Te . Secondly, let T be an element of T n Tc whose F is a face. Since T 62 Tc , all the vertices of T cannot belong to S. Since F is a face of Te , the vertices of F belong to S. Since F is a face of T, it comes T 2 Tf . Then, condition (ii) says that F belongs to oTc . Thirdly, since F 2 oTc , we have F 2 oT or F 2 oðT n Tc Þ, but in any case F cannot share more than one simplex in T n Tc . Besides, Lemma 10 says that oTe ¼ oTc , so F cannot share more than one simplex in Te , which raises a contradiction. Finally, it remains the proof of oðT n Tc [ Te Þ ¼ oT. • Let F 2 oT, we have 9 !T 2 T such that F 2 o{T}. – If T 62 Tc then F 62 oTc , so 8T 0 2 Te F 62 ofT 0 g. Hence, T is unique in T n Tc [ Te . – If T 2 Tc then F 2 oTc (because T is unique in Tc ) so, 9 !T 0 2 Te such that F 2 o{T 0 }. Since 8T 00 2 T n Tc F 62 ofT 00 g, T 0 is unique in T n Tc [ Te . At this point we have oT oðT n Tc [ Te Þ. • Let F 2 oðT n Tc [ Te Þ, we have 9 !T 2 T n Tc [ Te such that F 2 o{T}. – If T 62 Te then T 2 T n Tc and F 62 oTe . We have 8T 0 2 Tc F 62 fT 0 g, otherwise, T 0 is unique in Tc (because of the presence of T in T) so F 2 oTc , which is impossible because oTe ¼ oTc (Lemma 10). – If T 2 Te then F 2 oTe ¼ oTc (Lemma 10). So, 9 !T 0 2 Tc such that F 2 o{T 0 }. By construction, we have Te \ ðT n Tc Þ ¼ ;. Since T is unique in T n Tc [ Te , 8T 00 2 T n Tc F 62 ofT 00 g. Hence, T 0 is unique in T. It comes oðT n Tc [ Te Þ oT. h With the star operation, edges and faces can be swapped, the node S can be removed (only if S 62 oTc ) and the centroid C can be introduced. But, for a boundary node S 2 oX, the solution to get S 62 oTc is to connect all boundary faces F 2 oX to a ﬁctitious node (say node 0). With these (virtual) elements, the topology is without boundary (i.e. oT ¼ ;, Fig. 4). In this case, the faces of oX are the node 0 opposite faces. The introduction of node 0 allows the treatment of boundary nodes as internal ones. Thus, they do not have a particular treatment, which simpliﬁes the implementation of our dimension independent mesh generator. In particular, boundary nodes can be eliminated, moved or introduced if it does not induce volume changing (the coordinates of node 0 are locally the coordinates of node S).

0

Fig. 4. A topology without boundary and its virtual elements connected to node 0.

4960

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

Finally, this algorithm is more time consuming than a Delaunay mesh generator because it tries many more topologies. However, to improve a mesh by local optimizations turns out to be a good strategy in the following section. 2.4. Anisotropic adaptation Now that the minimal volume has been reached by topological optimization, we usually want to adapt the mesh to an isotropic or anisotropic metric ﬁeld, may it be Euclidian (i.e. constant over X) or not. Let ðN; TÞ be a mesh of X and let M be a continuous P1 metric ﬁeld on this mesh. We only use continuous P1 ﬁelds during mesh generation process, because they are easier to transport (than discontinuous P0 ﬁelds for example). To adapt the mesh to its metric ﬁeld, we proceed by the same local cut and paste technique. But in this case, we have to keep the minimal volume and to improve the quality of the elements. The anisotropic quality criterion [19] we use for an element T whose vertices are S0, . . . , Sd is c0

jT jM T hdM T

ð7Þ

;

where we compute the average metric

the Euclidian volume 2

the L average size

d 1 X MðS i Þ; d þ 1 i¼0 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ jT jM T ¼ jT j detðM T Þ;

MT ¼

hM T ¼

X 2 2 kS j S i kM T dðd þ 1Þ 06i Mx. This anisotropic quality criterion is Euclidian because a constant metric MT is considered over the element T. Note that c0 is such that (7) = 1 when T is equilateral (in the metric MT). Furthermore, in order to introduce or to remove enough nodes before optimizing the shape of the elements, a size criterion is also used d 1 min hM T ; . ð8Þ hM T The ﬁnal criterion c(T) is the worst between the quality criterion (7) and the size criterion (8). At this point, we can introduce an ordering relationship to compare local topologies: let T1 ¼ ðT i1 Þ16i6I 1 and T2 ¼ ðT i2 Þ16i6I 1 (such that 81 6 i1 < i2 6 I j cðT ij1 Þ 6 cðT ij2 Þ) be two local topologies of the same cavity (Fig. 3(b) and (c), for example) 8i < i0 cðT i1 Þ ¼ cðT i2 Þ; T1 < T2 () 9 i0 such that ð9Þ and cðT i10 Þ < cðT i20 Þ. Following this ordering relationship, the next algorithm ensure that the worst element of a maximal topology will be as good as possible. Instead, the use of a cost function [20], only provides an average result that enable degenerated or very poor elements. Moreover, compared to [20,21], which uses the same optimization technique, the only primitive needed here is the star operator.

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4961

Algorithm 2. The algorithm to make the mesh ðN; TÞ ﬁt the metric ﬁeld M is: while the topology T has been changed do for each node and each edge of the mesh do – cut the local topology around this node or this edge – among all minimal volume candidates (obtained by the star operator) paste one maximal topology in the ordering relationship between local topologies – transport the metric ﬁeld, if necessary done done The metric ﬁeld needs to be transported only if a new node C has beenP introduced. In this case, when C is n the isobarycenter of {S1, . . . , Sn}, the vertices of oTc , we take MðCÞ ¼ 1n i¼1 MðS i Þ and when C is the mid0 0 dle of [S, S ], we take M(C) = (M(S) + M(S ))/2. This choice is made for the sake of simplicity, but it remains inaccurate. The mesh we obtain is not directly perfect. By computing again the metric ﬁeld on this mesh and restarting the Algorithm 2, a good mesh can be iteratively produced by this adaptive process (as in Fig. 5 where the metric ﬁeld is deﬁned by [20]).

2.5. Conclusion The introduction of a metric M only changes volume and length calculation. This anisotropic mesh generation method is very robust and its implementation is less technical than for an hybrid Delaunay/frontal [22]. Furthermore, since the algorithm only uses local operations, this mesh generation and adaptation has been easily parallelized [23].

Fig. 5. Anisotropic adaptive process. (a) Initial mesh; ﬁrst step: (b) 1943 nodes in 73 s, (c) 25,066 nodes in 1301 s, and (d) 23,694 nodes in 352 s.

4962

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

The remainder of this paper focuses on the construction of a metric ﬁeld M in the a priori context of mesh generation.

3. Natural metric ﬁeld In this section, we want to address a diﬃcult problem in mesh generation: the meshing of thin and curved geometries (Fig. 6). For that purpose, we propose to take advantage of the anisotropic adaptivity driven by a metric ﬁeld M that places enough elements where it is necessary (i.e. through the thickness), while saving maximum nodes everywhere it is possible (i.e. in other directions), without much human intervention and only from geometrical informations. With several notations, the local anisotropy of a given geometry can be deﬁned (Section 3.1). Then, an approximation of the natural metric can be computed (Sections 3.2–3.4). Finally, the natural metric can be used in placing several element layers through the thickness (Sections 3.5 and 3.6). 3.1. Natural anisotropy deﬁnition Firstly, we need to introduce an ordering relationship between metrics: let M1 and M2 be two metrics ( 81 6 i < i0 k1i ¼ k2i ; ð10Þ M 1 < M 2 () 9 1 6 i0 6 d such that and k1i0 > k2i0 ; where kj1 6 6 kjd are the ordered eigenvalues of Mj. Secondly, the unit ball of a metric M is an ellipsoid EðMÞ ¼ fy 2 Rd such that kykM 6 1g. Thirdly, we denote by o12 X the median surface of X [7].

(a)

(b) Fig. 6. Initial mesh (a) and natural mesh (b) with 4 element layers (193 nodes) in 2D.

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4963

Fig. 7. Natural metric deﬁnition notations.

Finally, we deﬁne the natural metric on x 2 o12 X by the largest metric M such that x þ EðMÞ X (Fig. 7). We deﬁne the natural metric on x 2 X o12 X by the natural metric of x0 2 o12 X such that kx x0k2 is the distance of x to o12 X. If x0 is not unique, we take the averaged natural metric. This deﬁnition cannot be directly implemented because it would involve a inﬁnite number of operations. The proposed method to compute approximately the natural metric is only performed on the nodes S of the mesh. It can be divided in three steps: • detection of the local anisotropy by aggregation of elements around S (Section 3.2), • elliptic interpolation of this element patch (Section 3.3), • and interpolation error estimation to decide when to stop the aggregation (Section 3.4).

3.2. Successive neighborhoods To aggregate a collection of elements around a node S we proceed by successive neighborhoods. Therefore, we deﬁne Nk ðSÞ the kth order neighborhood around S (with k 2 N) by 8 < N0 ðSÞ ¼ fSg; S S ð11Þ Nl ðSÞ [ ðNl ðSÞ \ oXÞ : Nk ðSÞ ¼ NðNk1 ðSÞÞ n 06l kmax the elliptic interpolation error is too large, i.e. when 2 max kS 0 C k kMðSÞ 1 > 0.75. ð18Þ S0 2 N e k ðSÞ f k ðSÞ becomes really non-elliptic (75% is an In other words, we stop the neighborhood growth when N empirically good value for discriminating elliptic and non-elliptic neighborhoods). Then, the appropriate order is K = kmax (or K = kmin when kmin > kmax). Hence, the computed natural metric is M natural ðSÞ ¼ cK LðC K Þ1 with LðC K Þ ¼

Z

ð19Þ

ðy C K Þ ðy C K Þ dy

ð20Þ

CK ðSÞ

and

P cK ¼ P

S0 2 N e K ðSÞ S0 2 N e K ðSÞ

kS 0 C K k2LðCK Þ1 4

kS 0 C K kLðCK Þ1

.

ð21Þ

3.5. Several element layers through the thickness At this point, we know the local anisotropy of the geometry, especially the size and the direction of the thickness. Remember that we want to place enough elements where it is necessary (i.e. through the thickness) while saving maximum nodes everywhere it is possible (i.e. in other directions). 2 We denote by h2 1 P P hd the classiﬁed eigenvalues of Mnatural(S) and by R the rotation matrix such that 3 2 1 0 7 6 h2 7 6 1 7 > 6 .. 7R . M natural ðSÞ ¼ R6 ð22Þ 7 6 . 7 6 4 15 0 h2d

4966

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

In order to get n element layers through the thickness (n being provided by the end user), we have to divide the shortest size h1 by n. But, when the local geometry is isotropic, i.e. when the metric Mnatural(S) is isotropic, all the sizes h1, . . . , hd should be divided by n (in order to generate isotropic elements). In practice, the sizes h1, . . . , hp considered as the thickness are those between h1 and (1 + a)h1 where a > 0 (empirically, we use a = 1). Then, h1, . . . , hp are divided by n, while hp+1, . . . , hd remain unchanged. Finally, we obtain the partially divided natural metric 3 2 2 n 07 6 2 7 6 h1 7 6 7 6 . .. 7 6 7 6 7 6 2 n 7 6 7 6 2 7 > 6 hp 7R . M nnatural ðSÞ ¼ R6 ð23Þ 7 6 1 7 6 7 6 7 6 h2pþ1 7 6 7 6 .. 7 6 7 6 . 7 6 15 4 0 h2d 3.6. Numerical results The previous metric can perfectly treat three-dimensional curvatures (Fig. 10 with n = 6).

Fig. 10. Thin and curved geometry: initial mesh, natural mesh (5441 nodes) and a sectional view.

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4967

Fig. 11. Natural mesh (174,759 nodes) of a complex electrical device (courtesy of Schneider Electric) and a partial view.

Furthermore, the algorithm shows a great robustness, which enables the natural mesh generation of very complex 3D geometries (Fig. 11). In what concerns the computational cost of natural metric construction, the complexity of the algorithm is linear and the computation time on a 220 Mﬂops workstation is about twenty minutes for 105 nodes. Note that the anisotropy of these meshes does not induce any loss of accuracy in the solution given by our ﬁnite element solver [25]. Indeed, it respects the anisotropy of the solution which is often driven by the anisotropy of the geometry. 3.7. Conclusion Our method does not use skeletonization to extract the local thickness [26], and it does not need to extrude a mesh from its boundary or to determine opposite vertices [27,28], so as to get several element layers through the thickness. Instead, we automatically perform a metric ﬁeld that detects the local anisotropy of the geometry. The only a priori information needed to build this natural metric ﬁeld is: an initial mesh ðN; TÞ that deﬁnes the geometry of the part and n, the number of element layers the user wants through the thickness. To compute this natural metric, we gather an element collection around each node, which seems to be a good way to determine the local anisotropy. Then, among all the techniques to determine an elliptic interpolation of this element patch, the length tensor integration and inversion turn out to be particularly robust.

4968

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

Nevertheless, we proceed by successive neighborhoods to collect these element patches, which is certainly not the only possible method. However, we obtain satisfactory results, even on complex parts. 4. Multidomain metric ﬁeld At this point, we are able to generate a natural mesh of a single domain X. Now, let X be composed of several subdomains X = x1 [ [ xp. We want to replace the usual method, which consists in generating one mesh for each subdomain [3,29–31] by the generation of a single mesh which captures the interfaces between subdomains. Our main goal is then to represent internal surfaces in a 3D mesh, up to a desired precision and only from geometrical informations concerning the subdomains xi. Of course inside each subdomain, a natural mesh should be considered (Section 4.1). The ﬁrst idea to represent such internal surfaces, consists in imposing exactly these surfaces in the mesh, either by sub-mesh merging or by well known operations in Delaunay mesh generation technique [22]. However, keeping exactly this internal surfaces in the mesh represents a great constraint on all remeshing operations after. The alternative we develop here relies on a metric ﬁeld that tighten the mesh around internal surfaces (cf. Fig. 12). For that purpose, an interface tensor is added to the natural metric (Section 4.2). This tensor is computed from presence functions which is discussed in Section 4.3. 4.1. Subdomain natural metric The natural mesh of whole X is not very interesting. Instead, the natural metric must be computed on each subdomain: if S 2 xi, then the metric we compute on S is the natural metric of xi. Note that, if a node S belongs to several subdomains, then the natural metric in S does not really matter, because of interface reﬁnement in the following paragraph.

Fig. 12. A 2D multidomain mould mesh (1771 nodes) with a cavity, three prints and two regulation channels.

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4969

To extend the algorithm of the Section 3 onto xi X, it is suﬃcient to ignore the neighbors outside xi when we deﬁne Nk ðSÞ and also to ignore elements whose vertices are not all in xi when oTk ðSÞ is deﬁned. Then, we easily have Mnatural, a metric ﬁeld on X that is natural for each subdomain. Again, we can divide this natural metric on xi by ni, the number of element layers through xi. Note that, in the notation M nnatural , n may be diﬀerent on each subdomain. Hence, the mesh of X will be adapted diﬀerently on each subdomain x1, . . . , xp. But it is not enough to capture the interfaces between them. An additional term is needed and described in the following paragraph. 4.2. Interface tensor deﬁnition Let x X be a subdomain with non-empty interior. We denote by d(x, ox) the signed distance between x 2 Rd and ox (that is, d(x, ox) 6 0 when x 2 x and d(x, ox) > 0 when x 62 x). A strategy to tighten the mesh around the border of x X can be • to take a function g such that g = 1 inside x, g = 0 outside and 0 < g < 1 in the neighborhood of ox. For example 1 gex ðxÞ ¼ ð24Þ 1 þ eadðx;oxÞ with a 1 (Fig. 13) • and to triangulate its graph G ¼ x 2 Rdþ1 such that ðx1 ; . . . ; xd Þ 2 X and xdþ1 ¼ gðx1 ; . . . ; xd Þ ;

ð25Þ

so as to get a mesh that captures ox by orthogonal projection on X, • the triangulation of G being driven in Rdþ1 by the metric M nnatural ðx1 ; . . . ; xd Þ Mðx1 ; . . . ; xdþ1 Þ ¼ ; m2

ð26Þ

where m is the number of element layers the user wants around the interface ox.

1 0.8 0.6 0.4 0.2 0 2 2

1 1

0

0

–1

–1 –2 –2

Fig. 13. Function gex with dðx; oxÞ ¼ kx ck2 r.

4970

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

In fact, this strategy is not directly applied. Proposition 12. The induced metric by M on G (and thus on X by orthogonal projection) is M multidomain ðx1 ; . . . ; xd Þ ¼ M nnatural ðx1 ; . . . ; xd Þ þ m2 rgðx1 ; . . . ; xd Þ rgðx1 ; . . . ; xd Þ.

ð27Þ

That is, the previous strategy in Rdþ1 is equivalent to the use of metric (27) in X Rd . Proof. We denote by f the orthogonal projection on Rd : f(x1, . . . , xd+1) = (x1, . . . , xd). According to Riemannian geometry [15,16], the induced metric by M on X (with map f ) have components 1 1 M ijmultidomain ¼ ðofoxj Þ> M ofoxi . Since on G we have f1(x1, . . . , xd) = (x1, . . . , xd,g(x1, . . . , xd)) we get

of 1 oxi

og > ¼ ðd1i ; . . . ; ddi ; ox Þ . Hence, we obtain i > ! " og og og og ¼ d1j ; . . . ; ddj M nnatural ðd1i ; . . . ; ddi Þ> þ m2 M ijmultidomain ¼ d1j ; . . . ; ddj ; M d1i ; . . . ; ddi ; oxj oxi oxj oxi ! n "ij ij 2 ¼ M natural þ m ðrg rgÞ ;

which gives exactly the metric (27). h Note that the interface tensor $g $g is not a metric, because 0 is always an eigenvalue of x x with x 2 Rd when d P 2. But in addition with a metric, the result is also a metric. At this stage, the multidomain metric involves the computation of the costly function gex. In the next paragraph, gex is substituted by a cheapest function g such that g = 1 inside x, g = 0 outside and 0 < g < 1 in the neighborhood of ox. 4.3. The presence function and its gradient A subdomain x X can be represented by its presence function x ðxÞ

x

: X ! f0; 1g such that

¼ 1 () x 2 x.

ð28Þ 0

In practice, to compute eﬃciently the interface tensor, the function g taken into account is the P discontinuous interpolation of x . So, g: X ! [0, 1] is deﬁned (almost everywhere) by 8T 2 T

8x 2 T

gðxÞ ¼ gðT Þ ¼

jT \ xj . jT j

Proposition 13. The interpolation errors are !1=2 X kg x kL2 ðXÞ ¼ gðT Þð1 gðT ÞÞjT j ;

ð29Þ

ð30Þ

T 2T

kg

x kL1 ðXÞ

¼ 0.

So, considering g instead of small.

ð31Þ x

is conservative and precise when all elements T such that 0 < g(T) < 1 are

R R P Proof. We have kg x kL2 ðXÞ ¼ ð X ðg x Þ2 dxÞ1=2 ¼ ð T 2T T ðg x Þ2 dxÞ1=2 . Since g is constant on T R R R R we get T ðg x Þ2 dx ¼ gðT Þ2 jT j 2gðT Þ T x dx þ T 2x dx. As T x dx ¼ jT \ xj and 2x ¼ x , we R 2 2 obtain T ðg x Þ dx ¼ gðT Þ jT j 2gðT ÞjT \ xj þ jT \ xj. And using jT \ xj = g(T)jTj, it comes R 2 which proves Formula (30). Besides, we have T ðg x Þ dx ¼ gðT ÞjT jðgðT Þ 2gðT Þ þ 1Þ

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4971

T F n(F)

T(F)

T

ω

(a)

(b) 0

Fig. 14. Pixelization and P gradient notations.

kg

x kL1 ðXÞ

mula (31).

¼

R

X ðg

x Þ dx

¼

P

R

T 2T T ðg

x Þ dx

¼

P

jT \xj T 2T ð jT j jT j

jT \ xjÞ ¼ 0 which proves For-

h

To calculate g we pixelize or voxelize X up to a precision (the choice of this precision gives the thickness of the reﬁnement). Then, we light on the pixels or the voxels that belong to x (Fig. 14(a)) and we take 8T 2 T

jT \ xj lit voxels in T ¼ . jT j all voxels in T

ð32Þ

Hence, the only information we need about x is a function that says if a point is inside or not. The numerical gradient of g is a P0 discontinuous vector ﬁeld rg: X ! Rd deﬁned (almost everywhere) by 8T 2 T

8x 2 T rgðxÞ ¼ rgðT Þ ¼

X gðT ðF ÞÞ gðT Þ jF jnðF Þ jT ðF Þj þ jT j F 2oT

ð33Þ

with oT the set of T faces, T(F) the opposite element of T through F, and n(F) the outgoing normal of T on F (Fig. 14(b)). This P0 gradient needs boundary conditions. When F 2 oT, we take jT(F)j = 1, which corresponds to adiabatic boundary conditions. Since there are several presence functions g1, . . . , gp corresponding to diﬀerent subdomains, for each component 1 6 i 6 d we take the maximum

Fig. 15. The mesh that deﬁnes x, its voxelization and the initial mesh of X.

4972

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

! " rgðT Þi ¼ max rg1 ðT Þi ; . . . ; rgp ðT Þi ;

ð34Þ

which leads to a good treatment of triple points. Finally, we obtain a continuous P1 multidomain metric ﬁeld, deﬁned on each node S of the mesh by M multidomain ðSÞ ¼ M nnatural ðSÞ þ m2 rgðSÞ rgðSÞ with rgðSÞ ¼

X T 2TðSÞ

!, jT jrgðT Þ

X

ð35Þ

! jT j ;

T 2TðSÞ

where TðSÞ is the elements of T whose S is a vertex.

Fig. 16. Iterations between the metric generator and the mesh generator.

ð36Þ

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

Fig. 17. The ﬁnal mesh of X with internal surface ox captured (150,683 nodes), and a sectional view.

Fig. 18. A partial view of the previous sectional view.

4973

4974

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4.4. Numerical results We consider half of a symmetrical part x (Fig. 15(a)) in a rectangular domain X (Fig. 15(c)). This part is a biomedical application (courtesy of Sophysa), which is a benchmark for the multidomain ﬁnite element solver [1]. The voxelization of x (Fig. 15(b)) is accurate enough to keep all geometrical informations and is used to calculate eﬃciently the presence function of x. Starting from a very coarse mesh (Fig. 15(c)), we generate the multidomain metric and we adapt the mesh to this metric. This process is iterated several times (Fig. 16).

Fig. 19. The extracted meshes where g = 1, 0 < g < 1 and g = 0.

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

4975

At iteration 6, the boundary seems to be good (Fig. 17(a)) and we can see that the mesh is well adapted everywhere in X, not only on oX (Fig. 17(b)). The reﬁnement around the interface between x and X nx is really anisotropic (Fig. 18). Finally, we can extract the sub-meshes of x and Xnx and see the natural metric eﬀect (Fig. 19(a) and (c)). The elements T where 0 < g(T) < 1 are very thin (Fig. 19(b)). Hence, the interpolation error kg x kL2 ðXÞ is very low. 4.5. Conclusion To enrich a mesh around the interfaces between several subdomains, we use presence functions g to locate interfaces and their gradients $g to compute the interface tensor $g $g. Furthermore, this tensor gives a multidomain metric that works properly when it is used in addition to the natural metric we deﬁned previously. The only a priori information we need to build the multidomain metric ﬁeld is: for each subdomain, its natural metric and a function that returns if a point is inside or not. Like [32,33] we progressively detect and reﬁne the elements around interfaces. But, we can deal with complex geometries and the reﬁnement is anisotropic. Although this technique is still little node consuming, we obtain good results.

5. Conclusion The dimension independent mesh generator we introduce in this paper shows a great robustness in adapting an unstructured mesh to any anisotropic metric ﬁeld. Besides, we deﬁne a multidomain metric ﬁeld that captures internal surfaces between subdomains. This metric is computed only from geometry information and uses a natural metric ﬁeld inside each subdomain. Since this natural metric can treat thinness, curvature and changing thickness, we are now able to treat complex parts, whose meshes cannot be generated manually. The ﬁnite element solvers we use with the meshes discussed below, do not show obvious sensitivity to nearly degenerated elements (when they are stretched along the ﬂow direction). The metric ﬁeld method we propose here may be used with other anisotropic mesh generators. Although these algorithms were developed for material processing simulations, they are applicable to other numerical ﬁelds.

Acknowledgement The authors thank gratefully to the referee for helpful suggestions, and to the Rem3D Consortium partners for their support.

References [1] S. Batkam, T. Coupez, A new multi-domain approach for the thermal problem solution in 3D mold ﬁlling, in: ESAFORM Conference on Material Forming, Liege, 2001, pp. 31–34. [2] J. Bruchon, T. Coupez, A study of the 3D polymer foam formation based on the simulation of anisothermal bubble expansion, Mecanique Indust. 4 (2003) 331–338. [3] J. Barboza, L. Fourment, J.-L. Chenot, Contact algorithm for 3D multi-bodies problems: Application to forming of multimaterial parts and tool deﬂection, in: 5th World Congress on Computational Mechanics, Vienna, 2002.

4976

C. Gruau, T. Coupez / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4951–4976

[4] P.O. Bouchard, F. Bay, Y. Chastel, I. Tovena, Crack propagation modelling using an advanced remeshing technique, Comput. Meth. Appl. Mech. Engrg. 189 (2000) 723–742. [5] C. Gruau, T. Coupez, Anisotropic and multidomain mesh automatic generation for viscous ﬂow ﬁnite element method, in: N.-E. Wiberg, P. Dı`ez, Adaptive Modelling and Simulation, Go¨teborg, 2003. ISBN: 84-95999-30-7. [6] M. Bellet, C. Aliaga, O. Jaouen, Finite elements for a thermomechanical of analysis of solidiﬁcation processes, in: 9th International Conference on Modeling of Casting, Welding and Advanced Solidiﬁcation Processes, Aachen, 2000, pp. 10–17. [7] P.J. Frey, P.-L. George, Mesh generation—Application to ﬁnite elements, Hermes, 2000. [8] J. Dompierre, M.-G. Vallet, Y. Bourgault, M. Fortin, W.G. Habashi, Anisotropic mesh adaptation: Towards user-independent, mesh-independent and solver-Independent CFD. Part III: Unstructured meshes, Int. J. Numer. Meth. Fluids 39 (2002) 675–702. [9] J. Peraire, K. Morgan, Unstructured mesh generation including directional reﬁnement for aerodynamic ﬂow simulation, Finite Elements Anal. Des. 25 (1997) 343–355. [10] R. Lo¨hner, C. Yang, E. On˜ate, S. Idelssohn, An unstructured grid-based, parallel free surface solver, Appl. Numer. Math. 31 (1999) 271–293. [11] F. Alauzet, P.J. Frey, B. Mohammadi, Anisotropic mesh adaptation. Application to transient CFD problem, in: Second MIT Conference on Computational Fluid and Solid Mechanics, Boston, 2003, pp. 1851–1854. [12] M.J. Castro-Dı´az, F. Hecht, B. Mohammadi, O. Pironneau, Anisotropic unstructured mesh adaptation for ﬂow simulations, Int. J. Numer. Meth. Fluids 25 (1997) 475–491. [13] T. Coupez, Grandes transformations et remaillage automatique, PhD Thesis, France, 1991, pp. 123–210. [14] P.D. Zavattieri, E.A. Dari, G.C. Buscaglia, Optimization strategies in unstructured mesh generation, Int. J. Numer. Meth. Eng. 39 (1996) 2055–2071. [15] M.P. doCarmo, Riemannian Geometry, Birkha¨user, 1988. [16] J. Jost, Riemannian Geometry and Geometric Analysis, Springer-Verlag, 1998. [17] P.M. Morse, H. Feshbach, Stokes theorem, in: Methods of Theoretical Physics, Part I, New York, 1953, p. 43. [18] T. Coupez, Generation de maillage et adaptation de maillage par optimisation locale, Revue europeenne des elements ﬁnis 9 (2000) 403–423. [19] J. Dompierre, M.-G. Vallet, P. Labbe´, F. Guibault, On simplex shape measures with extension for anisotropic meshes, in: Workshop on Mesh Quality and Dynamic Meshing, Livermore, 2003, pp. 46–71. [20] C. Bottasso, On the choice of the optimization strategy and of the objective function for metric-driven mesh adaptation, in: International Conference on Adaptative Modelling and Simulation, Go¨teborg, 2003. [21] C.C. Pain, A.P. Umpleby, C.R.E. de Oliveira, A.J.H. Goddard, Tetrahedral mesh optimisation and adaptativity for steady-state and transient ﬁnite element calculations, Comput. Methods Appl. Mech. Engrg. 190 (2001) 3771–3796. [22] P.-L. George, H. Borouchaki, Delaunay triangulation and meshing. Application to ﬁnite elements, Hermes, 1998. [23] T. Coupez, H. Digonnet, R. Ducloux, Parallel meshing and remeshing, Appl. Math. Modell. 25 (2000) 153–175. [24] S.G. Advani, C.L. Tucker, The use of tensors to describe and predict ﬁber orientation in short ﬁber composites, J. Rheol. 31 (1987) 751–784. [25] T. Coupez, E. Bigot, 3D anisotropic mesh generation and adaptation with applications, in: European Congress on Computational Methods in Applied Sciences and Engineering, Barcelona, 2000. [26] K.-F. Tchon, M. Khachan, F. Guibault, R. Camarero, Constructing anisotropic geometric metrics using octrees and skeletons, in: 12th International Meshing Roundtable, Santa Fe, 2003, pp. 293–304. [27] R.V. Garimella, M.S. Shephard, Boundary layer meshing for viscous ﬂows in complex domains, in: 7th International Meshing Roundtable, Dearborn, 1998, pp. 107–118. [28] R.V. Garimella, Anisotropic tetrahedral mesh generation, PhD thesis, Rensselaer Polytechnic Institute, 1998. [29] S. Chakravarthy, A uniﬁed-grid ﬁnite volume formulation for computational ﬂuid dynamics, Int. J. Numer. Meth. Fluids 31 (1999) 309–323. [30] T. Tezduyar, Y. Osawa, The multi-domain method for computation of the aerodynamics of a parachute crossing the far wake of an aircraft, Comput. Methods Appl. Mech. Engrg. 191 (2001) 705–716. [31] M. Lesoinne, C. Farhat, Geometric conservation law for ﬂow problems with moving boundaries and deformable meshes, and their impact on aeroelastic computations, Comput. Methods Appl. Mech. Engrg. 134 (1996) 71–90. [32] T. Stalicky´, H.-G. Roos, Anisotropic mesh reﬁnement for problems with internal and boundary layers, Int. J. Numer. Meth. Engrg. 46 (1999) 1933–1953. [33] T. Tezduyar, S. Aliabadi, M. Behr, Enhanced-discretisation interface-capturing technique (EDICT) for computation of unsteady ﬂows with interfaces, Comput. Methods Appl. Mech. Engrg. 155 (1998) 235–248.